# A comparison of different clustering approaches for high-dimensional
presence-absence data

Gabriele d’Angella111Università di Bologna, Dipartimento di Scienze
Statistiche<EMAIL_ADDRESS>and Christian
Hennig222Università di Bologna, Dipartimento di Scienze Statistiche,
<EMAIL_ADDRESS>
###### Abstract

Presence-absence data is defined by vectors or matrices of zeroes and ones,
where the ones usually indicate a “presence” in a certain place. Presence-
absence data occur for example when investigating geographical species
distributions, genetic information, or the occurrence of certain terms in
texts. There are many applications for clustering such data; one example is to
find so-called biotic elements, i.e., groups of species that tend to occur
together geographically. Presence-absence data can be clustered in various
ways, namely using a latent class mixture approach with local independence,
distance-based hierarchical clustering with the Jaccard distance, $K$-modes, a
density-based approach, or also using clustering methods for continuous data
on a multidimensional scaling representation of the distances. These methods
are conceptually very different and can therefore not easily be compared
theoretically. We compare their performance with a comprehensive simulation
study based on models for species distributions.

This has been accepted for publication in Ferreira, J., Bekker, A., Arashi, M.
and Chen, D. (eds.) Innovations in multivariate statistical modelling:
navigating theoretical and multidisciplinary domains, Springer Emerging Topics
in Statistics and Biostatistics.

Keywords: multidimensional scaling; biogeography; cluster analysis; simulation
study; benchmarking; Jaccard’s distance.

##  1 Introduction

Presence-absence data comprises observations that are vectors or matrices of
zeroes and ones, where the ones usually indicate a “presence” in a certain
place. Presence-absence data are often high-dimensional (the number of
“places” may be large compared to the number of observations), and occur in a
large range of applications, for example when investigating geographical
species distributions, genetic information, or the occurrence of certain terms
in texts.

Often there is an interest in clustering presence-absence data. An example for
this is the search for “biotic elements”, which are groups of species showing
very similar distribution areas that can constitute evidence for the existence
of areas of endemism generated by the formation of barriers over geologic time
periods [16, 17].

Here we present a simulation study for comparing different approaches to
clustering such data. This is based on a model for generating artificial but
realistic presence-absence data with known clusters, the recovery of which can
be compared. Not only are we interested in ranking methods, we also
investigate characteristics of the data (such as degree of cluster overlap)
that drive the methods’ performances; the comparison of methods may strongly
depend on such characteristics. We will also explain in detail why some
methods do not perform that well. See [34] for guidelines on benchmark studies
to compare clustering methods.

There are various approaches to clustering presence-absence data. Two major
approaches are (1) the use of distance-based clustering methods such as
average linkage or partitioning around medoids [14] using distance measures
between the observations [9, 31], and (2) latent class mixture models [26].
Some of the most popular clustering methods such as K-Means [23] and model-
based clustering based on Gaussian mixtures [30] require Euclidean data and
cannot directly be applied to presence-absence data. As approach (3), Hausdorf
and Hennig [17, 18] have proposed to use methods for Euclidean data on the
output of a Multidimensional Scaling (MDS) [6] specifically with presence-
absence data for biotic element analysis. Furthermore, [3] generalised their
density-based clustering in this way to categorical data.

MDS is a set of techniques that generate a Euclidean representation so that
the resulting Euclidean distances approximate a given usually non-Euclidean
distance structure. Low dimensionality of the representation is often
desirable. Such a representation is often used for visualization (for which
low dimensionality is needed), but can also be used for making distance data
accessible to more elaborate statistical methodology that requires Euclidean
input, as is done here for cluster analysis.

A disadvantage of this approach is that the MDS generally will lose some
information in the original data (the lower the dimension, the larger the
loss). A key interest in the present study is whether this approach can
compete with the more direct approaches (1), and (2).

The relationship between clustering and scaling techniques has also been
discussed in de Leeuw and Heiser [11] and in Kruskal [25]. Simultaneous use of
MDS and clustering is treated by Desarbo [13] and Oh and Raftery [28]. Some
biogeographical work related to our study has been done by Vavrik [35]
(comparison of clustering methods for fossil data) and Ulrich and Gotelli [33]
(null models for simulating presence-absence data).

Section 2 has a basic description of the data, distance, and MDS. Section 3
introduces the involved clustering methods. Section 4 describes in detail how
the data for the simulations were generated. Section 5 discusses the results
of the study, and Section 6 concludes the paper.

##  2 Data and preprocessing

Data here are $m$-dimensional binary: $\mathbf{x}_{i}=(x_{i1},\ldots,x_{im}),\
i=1,{\dots},n$, for all $j:\ x_{ij}\in\\{0,1\\}$. The data are modelled with a
specific application area in mind, which is analysing presence (1) and absence
(0) of species in geographical regions, where $i$ indexes the species and $j$
indexes the region. Clusters are called “biotic elements” and are of interest
because they provide insight into natural history [16, 18]. Our results may
however also be informative for presence-absence data in other applications.

Methods are compared regarding their capability to retrieve true clusters in
the data as set in the simulation. The cluster membership of the observations
is denoted as $\mathbf{c}=(c_{1},\ldots,c_{n})$, where $c_{i}=k$ denotes that
observation $i$ is in cluster $k\in\\{1,\ldots,K\\}$.

Some of the compared clustering methods operate on distances. As distance we
use the popular Jaccard distance [22, 9, 12], for discussion and alternatives
see [9, 31, 14, 19]:

$d_{J}(\mathbf{x}_{1},\mathbf{x}_{2})=1-\frac{\sum_{j=1}^{m}\mathbbm{1}(x_{1j}=1\land
x_{2j}=1)}{\sum_{j=1}^{m}\mathbbm{1}(x_{1j}=1\lor x_{2j}=1)}$ (2.1)

where $\mathbbm{1}$ denotes the indicator function.

A topic of key interest here is the performance of clustering techniques for
Euclidean data that use an MDS representation of the distances as input. In
this study, two MDS methods were applied: classical scaling [32] as
implemented in the R-function cmdscale, and ratio MDS [6, 7] as implemented in
the R-package smacof, function mds. Two and three-dimensional MDS solutions
were used in the simulation study. Although this seems low, the results (as
provided in Section 5) were good enough that there is not much room for
improvement by higher dimensional solutions (in [28], lower-dimensional
proxies yielded clustering solutions almost as good as those exploiting
higher-dimensional configurations). An advantage of using a low dimensional
MDS is the ease of visualization of the data. However, in more general
applications, using a higher dimensional MDS solution may be helpful in case
that too much information is lost by a low dimensional solution.

However, there are no indications against the setting of $p>3$ when it comes
to the application of techniques like KM or GMM and this should be considered
in case the comparison with other clustering recovery methods shows that too
small a $p$ might lose too much information. In this project we deal with
multivariate binary data: not only does MDS constitute a powerful
visualization tool, it also enables us to resort to methods requiring
continuous data input for our clusters recovery purposes. These methods are
introduced in the next subsection.

##  3 Clustering methods

We used the following clustering methods.

###  3.1 Latent class analysis

Latent class analysis operates directly on the presence-absence data. The
method used here is the simplest of those that are referred to in the
literature as “latent class analysis” [26, chapter 3]. It models the data as a
mixture of locally independent Bernoulli distributions, i.e., the different
variables/regions are assumed independent within clusters, although they can
depend over the whole dataset. Data are modelled as i.i.d. according to the
density [5, chapter 6]:

$\displaystyle
f_{\mathbf{\eta}}(\mathbf{x}_{i})=\sum_{k=1}^{K}\pi_{k}\prod_{j=1}^{m}\theta_{jk}^{x_{ij}}(1-\theta_{jk})^{1-x_{ij}},$
(3.1)

where $\theta_{jk}$ is the probability of positive response for variable $j$
in cluster $k$, namely the probability that a species in the $k^{th}$ group
inhabits the $j^{th}$ cell. The parameter vector
${\mathbf{\eta}}=((\pi_{1},\theta_{11},\ldots,\theta_{1m}),\ldots,(\pi_{K},\theta_{1K},\ldots,\theta_{mK}))$
stores mixture proportions and probabilities for presence per variable and
cluster. The parameters can be fitted by Maximum Likelihood, using the EM-
algorithm [36]. Observations can then be assigned to clusters by maximizing
the estimated posterior probability of the observations to belong to the
clusters [26, 5]. The Bayesian Information Criterion (BIC) can be used to
estimate the number of clusters, but in our study the number of clusters was
taken as fixed and known. In our study we used the R-function poLCA in the
R-package with the same name for computing this kind of latent class analysis.

###  3.2 Methods operating on distances

Some methods were used that take the Jaccard distance matrix as input, namely
standard Single, Average, and Complete Linkage clustering [14, Chapter 4] with
the dendrogram cut in such a way that the required number of clusters is
produced. These were computed using the R-function hclust. Furthermore we used
Partitioning Around Medoids (PAM; [24]) as computed by the R-function pam, and
the K-modes algorithm [20], which tries to optimise the PAM objective function
using the simple matching distance as computed by R-function kmodes in package
klaR with parameter fast set to FALSE to avoid error messages.

###  3.3 Methods operating on Euclidean data

For clustering the Euclidean MDS-output, we used two of the most popular
clustering methods, namely K-means [23] and the Gaussian mixture model [14,
chapter 6]. K-means was computed by the R-function kmeans (using parameters
nstart = 100, iter.max = 100 in order to allow for a more stable performance
than granted by the default values), the Gaussian mixture model was fitted by
the function Mclust in package mclust [30]. More precisely, data is assumed to
be generated by a model with density

$\displaystyle
f_{\mathbf{\eta}}(\mathbf{x}^{\texttt{mds}}_{i})=\sum_{k=1}^{K}{\pi_{k}}\phi_{\mathbf{a}_{k},\mathbf{\Sigma}_{k}}(\mathbf{x}^{\texttt{mds}}_{i}),$
(3.2)

where $K$ is the number of mixture components,
$\phi_{\mathbf{a},\mathbf{\Sigma}}$ is the density of the multivariate
Gaussian distribution with mean vector $\mathbf{a}$ and covariance matrix
$\mathbf{\Sigma}$, $\pi_{1},\ldots,\pi_{K}$ are the mixture proportions, and
$\eta$ is a vector collecting all parameters. These parameters can be
estimated by Maximum Likelihood as implemented in the EM-algorithm. mclust
provides several constrained covariance matrix models (e.g., all covariance
matrices equal) besides a fully flexible model, and the mclust\- software uses
the Bayesian Information Criterion (BIC) to select the best one [30].
Observations are assigned to clusters by

$\displaystyle\hat{P}(c_{i}=k|\mathbf{x}^{\texttt{mds}}_{i})=\frac{{\hat{\pi}_{k}}\phi_{\mathbf{\hat{a}}_{k},\mathbf{\hat{\Sigma}}_{k}}(\mathbf{x}^{\texttt{mds}}_{i})}{f_{\mathbf{\hat{\eta}}}(\mathbf{x}^{\texttt{mds}}_{i})}$
(3.3)

and estimating the cluster membership labels as follows:

$\displaystyle\hat{c_{i}}=\operatorname*{arg\,max}_{k\in\\{1,\dots,K\\}}\hat{P}(c_{i}=k|\mathbf{x}^{\texttt{mds}}_{i}).$
(3.4)

The BIC is also often used to estimate the number of clusters, but this is
treated as fixed and known here, which also means that there is no need to use
one of the many available methods (that often lead to contradictory results)
to decide the number of clusters for the other methods [15].

Furthermore, [3] have proposed to run the density-based clustering method
pdfCluster (R-function and package name) [2] on MDS-output in order to handle
categorical or binary data. The method estimates the density by kernel
methods. Clustering is then performed by finding density level sets using
neighbourhood graphs. As opposed to the other methods in the study, pdfCluster
will implicitly determine the number of clusters rather than allowing for it
to be fixed. In order to avoid errors, we used the parameter settings
`graphtype="pairs"` and `n.grid` equal to the number of observations.

K-means, the Gaussian mixture model, and pdfCluster are combined with the
different MDS techniques and two dimensionalities of the MDS output. The
combinations are called kc2, ks2, kc3, ks3, gc2, gs2, gc3, gs3, pdfc2, pdfs2,
pdfc3, pdfs3 in the results section. “k”, “g”, “pdf” refer to K-means, the
Gaussian mixture, and pdfCluster, respectively. “c” and “s” denote cmsscale,
and smacof’s mds, respectively. The number in the end is the number of MDS
dimensions used.

##  4 The simulation

In order to compare the methods, we simulate datasets of $n$ species, each of
which is represented by a presence-absence vector of length $m$, and cluster
them.

Simulations enable us to investigate the features of the data that drive the
performance of the clustering methods. Real biogeographical data sets can make
this assessment difficult, as it is hard to separate results due to
singularities in data from those due to methods properties [8]. Among the
plethora of aspects that can be considered, in this project we examined the
level of clusters’ overlap, the cluster size, the width of the cluster
specific areas and the number of clusters: in particular, we compared a
situation with three clusters endowed with their specific geographic area with
a situation where on top of these three groups a cluster of universal
spreaders is added, as often exist in reality [10]. Such a cluster consists of
species that are widespread across the whole map. These species are not
informative regarding biotic elements, however they qualify as clusters in a
data analytic sense, namely as a set of observations that behaves in a certain
way. In the following, we will call the clusters made up of species that are
not universal spreaders “proper”, because these species inhabit a specific
group of cells and can signal the existence of a biotic element.

The 24 parameter combinations for the simulated data are described below,
after the presentation of the algorithm that was used to generate this kind of
data.

###  4.1 Data generation

The algorithm implemented to simulate each species resorts repeatedly to a
random draw from a distribution that chooses one category of the $m$ possible
ones, each of which has its own probability of being picked. At each step, one
category has to be picked, and therefore these probabilities sum to one. This
distribution is called here “categorical distribution”. Let $j=1,{\dots},m$ be
the categories. Each category has a probability $p_{j}$ of being picked and
$\sum_{j=1}^{m}p_{j}=1$. The probability mass function of our categorical
distribution is:

$f_{\mathbf{p}}(y)=\prod_{j=1}^{m}{p_{j}}^{[y=j]}$ (4.1)

where $[y=j]$ is equal to 1 if $y=j$ and 0 otherwise.

The species in our data set inhabit the geographic units (cells) of a map. To
each of the $m$ categories corresponds a cell: at every drawing step, a
category is picked and its associated binary indicator is set to 1, meaning
that the species is present in this geographic unit of the map. A value of 0
means that the cell has not been drawn and the species is absent.

Each cluster will have its specific area on the map, consisting of a set of
cells whose probability to be inhabited by the species of the cluster will be
larger than elsewhere. To simplify, there will be a unique probability inside
these areas ($p_{in}$) and a unique probability outside of them ($p_{out}$),
with $p_{in}>p_{out}$. Once the cluster specific areas have been determined,
the probability vectors specific to each cluster can be produced. Species
belonging to a given cluster will mainly inhabit the cluster specific area,
while their occurrence outside of it will be rarer. The rareness of species
outside of their specific cluster area will depend on the level of overlap
$\omega$, see below.

We define the following notation:

  1. (a)

observations $\mathbf{x}_{i}=(x_{i1},\ldots,x_{im})$ taking values in the
product set $\\{0,1\\}^{m}$;

  2. (b)

cluster labels $k=1,{\dots},K$;

  3. (c)

species range $r_{i}=\sum_{j=1}^{m}x_{ij}$, for $i=1,\dots,n$;

  4. (d)

vector of species cluster labels $\mathbf{c}=(c_{1},\ldots,c_{n})$;

  5. (e)

vector of cell cluster labels $\mathbf{d}=(d_{1},\ldots,d_{m})$. A cell is
given the label $k$ when it belongs to the specific area of a cluster, i.e.,
the set of cells with elevated probability for cluster $k$. The probability to
be inhabited by a species belonging to that cluster is $p_{in}$. Each cell can
belong to at most one (exactly one in this study) cluster specific area.

The following parameters of the simulation are fixed in advance:

  * •

the number of species $n$ and the number of cells $m$. $n$ was either 300 or
400, see below. $m$ was 60;

  * •

the level of overlap $\omega=\frac{p_{out}}{p_{in}}$. Three levels of overlap
are used, i.e. low ($\omega=0.05$), medium ($\omega=0.2$) and high
($\omega=0.4$);

  * •

the number of clusters $K$. This was either 3 (with $n=300$), or 4 (with
$n=400$). There were always 3 proper clusters, and a cluster of universal
spreaders was added or not;

  * •

the number of species belonging to each cluster $k$, i.e., $|\\{i:c_{i}=k\\}|$
for $k=1,\dots,K$. Two possibilities are investigated: one with all clusters
having 100 observations and one with clusters of various sizes (a small
clusters with 50 species, one with 100 and one with 150 species). If a group
of universal spreaders is present, it consists of 100 observations;

  * •

the number of cells constituting each cluster specific area on the map, i.e.,
$|\\{j:d_{j}=k\\}|$ for $k=1,\dots,K$. Two possibilities are investigated, one
with equal-sized specific areas (all areas 20 cells) and one with
heterogeneous specific areas (10, 20, and 30 cells). With unequal area sizes
and unequal cluster sizes, cluster sizes were chosen in line with the area
sizes, i.e., the larger the cluster specific area, the larger the cluster.

The range $r_{i}$ of the $i^{th}$ species is instead randomly drawn. If the
species belongs to a proper cluster, it is uniformly drawn between 1 and the
size of the cluster specific area. If it is a universal spreader ($u.s.$), its
range is at least as large as the widest cluster specific area on the map,
namely

$\displaystyle r_{u.s.}=\max_{k}|\\{j:d_{j}=k\\}|,$

and the maximum is $m$. As a consequence, a universal spreader will always
cover a region of the map that involves more cells than the habitat of any
proper species. These choices have been made inspired by real datasets of this
kind that we have seen; particularly, species ranges are highly variable.

Here is the algorithm for the $i^{th}$ species draw, provided that it is not a
universal spreader:

  1. 1.

pick $r_{i}$ according to a discrete uniform distribution with domain
$\\{1,\dots,|\\{j:c_{i}=d_{j}\\}|\\}$;

  2. 2.

set $\mathbf{x}^{0}_{i}=\mathbf{0}_{m}$, namely the observation is initialized
as a vector of zeros of length $m$;

  3. 3.

for $t=0,\dots,r_{i}-1$:

    1. (a)

compute $p^{t}_{in,c_{i}}$ and $p^{t}_{out,c_{i}}$ such that:

$\left\\{\begin{array}[]{ll}\frac{p^{t}_{out,c_{i}}}{p^{t}_{in,c_{i}}}=\omega\\\
|\\{j:d_{j}=c_{i}\land x^{t}_{ij}=0\\}|\cdot{p^{t}_{in,c_{i}}}+|\\{j:d_{j}\neq
c_{i}\land x^{t}_{ij}=0\\}|\cdot{p^{t}_{out,c_{i}}}=1\end{array}\right.$

    2. (b)

specify the vector $\mathbf{p}^{t}_{i}$ according to the following rule:

$p^{t}_{ij}=\left\\{\begin{array}[]{ll}p^{t}_{in,c_{i}},&\text{ if
}c_{i}=d_{j}\land x^{t}_{ij}=0\\\ p^{t}_{out,c_{i}},&\text{ if }c_{i}\neq
d_{j}\land x^{t}_{ij}=0\\\ 0&\text{ otherwise }\end{array}\right.$

    3. (c)

draw j from a categorical distribution (eq. 4.1) with probability vector
$\mathbf{p}^{t}_{i}$;

    4. (d)

set $x^{t+1}_{ij}=1$.

  4. 4.

set the $i^{th}$ observation in the data set equal to
$\mathbf{x}^{r_{i}}_{i}$.

Note that the probability vector $\mathbf{p}^{t}_{i}$ is updated after each
draw in such a way that every category can only be drawn once. Step 3 makes
sure that $\omega$ (i.e., the chosen ratio between the probability to draw a
cell without and within the cluster specific area) is respected by the
probability vector $\mathbf{p}^{t}_{i}$. Species generated in this way may
well have presences outside the cluster specific area, and absences within it,
although the probability to be present within it is larger. This reflects
realistically observed patterns.

The algorithm for the draw of the $i^{th}$ observation being a universal
spreader is slightly different:

  1. 1.

set $\mathbf{x}^{0}_{i}=\mathbf{0}_{m}$, namely the observation is initialized
as a vector of zeros of length $m$;

  2. 2.

pick $r_{i}$ according to a discrete uniform distribution with domain
$\\{r_{u.s.},\dots,m\\}$;

  3. 3.

for $t=0,\dots,r_{i}-1$:

    1. (a)

pick $j$ according to a discrete uniform distribution with domain
$\\{1,\dots,m\\}\setminus\\{j:x^{t}_{ij}=1\\}$;

    2. (b)

set $x^{t+1}_{ij}=1$.

  4. 4.

set the $i^{th}$ observation in the data set equal to
$\mathbf{x}^{r_{i}}_{i}$.

###  4.2 Scenarios

24 simulation scenarios have been used, which are listed in Table 1.

Table 1: The 24 simulation scenarios resulting from the combination of the investigated data parameters. scenario | $\omega$ | $K$ | sizes | areas | case | $\omega$ | $K$ | sizes | areas
---|---|---|---|---|---|---|---|---|---
1 | 0.05 | 3 | $=$ | $=$ | 13 | 0.05 | 3 | $=$ | $\neq$
2 | 0.20 | 3 | $=$ | $=$ | 14 | 0.20 | 3 | $=$ | $\neq$
3 | 0.40 | 3 | $=$ | $=$ | 15 | 0.40 | 3 | $=$ | $\neq$
4 | 0.05 | 3+$u.s.$ | $=$ | $=$ | 16 | 0.05 | 3+$u.s.$ | $=$ | $\neq$
5 | 0.20 | 3+$u.s.$ | $=$ | $=$ | 17 | 0.20 | 3+$u.s.$ | $=$ | $\neq$
6 | 0.40 | 3+$u.s.$ | $=$ | $=$ | 18 | 0.40 | 3+$u.s.$ | $=$ | $\neq$
7 | 0.05 | 3 | $\neq$ | $=$ | 19 | 0.05 | 3 | $\neq$ | $\neq$
8 | 0.20 | 3 | $\neq$ | $=$ | 20 | 0.20 | 3 | $\neq$ | $\neq$
9 | 0.40 | 3 | $\neq$ | $=$ | 21 | 0.40 | 3 | $\neq$ | $\neq$
10 | 0.05 | 3+$u.s.$ | $\neq$ | $=$ | 22 | 0.05 | 3+$u.s.$ | $\neq$ | $\neq$
11 | 0.20 | 3+$u.s.$ | $\neq$ | $=$ | 23 | 0.20 | 3+$u.s.$ | $\neq$ | $\neq$
12 | 0.40 | 3+$u.s.$ | $\neq$ | $=$ | 24 | 0.40 | 3+$u.s.$ | $\neq$ | $\neq$

Equal ($=$) group sizes imply that there are 100 observations per cluster (be
it a proper one or a group of universal spreaders). When the group sizes are
different ($\neq$), there is a cluster with 50 species, one with 100, and one
with 150 (there are always 100 $u.s.$ species). When the areas are equal
($=$), the 60 cells are evenly divided across the clusters; when the areas
differ ($\neq$), the smallest cluster specific area is made up of 10
cells/variables and the biggest one consists of 30 cells/variables.

(a) Heatmap for case 1: no $u.s.$ group, equal sizes, equal areas, low
overlap.

(b) Heatmap for case 3: no $u.s.$ group, equal sizes, equal areas, high
overlap

(c) Heatmap for case 4: with $u.s.$ group, equal sizes, equal areas, low
overlap.

(d) Heatmap for case 13: no $u.s.$ group, equal sizes, different areas, low
overlap.

Figure 1: Presence-absence heatmaps for cases 1, 3, 4 and 13. Rows correspond
to species, columns correspond to cells. Colored bars represent species
clusters (rows) and cluster specific areas (columns).

We visualize the spread of the species in four exemplary cases using heatmaps.
They represent the species in the rows and the cells in the columns: in each
row, grey units represent cells where the species is present. In Figure 1(a)
the overlap is low and all groups have the same number of members and cluster
specific cells. In Figure 1(b) the high level of overlap makes it rather tough
to distinguish the cluster specific areas. Case 4 (Figure 1(c)) differs from
Case 1 in that it includes a group of universal spreaders. In Case 13 (Figure
1(d)), clusters have specific areas of differing sizes, therefore species
belonging to the green cluster will be able to inhabit up to thirty cells.

##  5 Results

###  5.1 General results

The simulations were evaluated using the Adjusted Rand Index [21] in order to
compare clusterings generated by the methods to the true clusterings, a
measure of similarity between data clusterings that ranges from -1 to 1: two
independent random partitions are expected to return an ARI of zero, whereas
an ARI equal to 1 implies identical clusterings. Not only did we try to
investigate what methods led to a satisfactory clustering recovery, but we
also checked what data features had a significant impact on the methods’
performance.

Each of the 24 scenarios reported in Table 1 was simulated ten times. This
configuration generated a data set with 3120 rows: thirteen methods applied on
24 scenarios, each simulated ten times.

(a) Interaction plot showing the average ARI values (vertical axis) for the
three levels of overlap ($\omega$).

(b) Interaction plot showing the average ARI values (vertical axis) depending
on the presence (yes) or absence (no) of universal spreaders in the simulation
setup.

(c) Interaction plot showing the average ARI values (vertical axis) with equal
or unequal cluster sizes.

(d) Interaction plot showing the average ARI values (vertical axis) with equal
or unequal cluster specific area widths.

Figure 2: Interaction plots of the average ARI values in the 24 cases (each
simulated ten times). The lines referring to each of the twelve clustering
recovery methods are indicated in the legends.

Figure 2 shows interaction plots of the results. They show the ARI means on
the vertical axis and the levels of the various data features on the
horizontal axis. The different lines in the plots refer to the clustering
methods. Figure 3 shows boxplots of the ARI values achieved by the methods by
overlap and presence of universal spreaders.

Unsurprisingly, methods are generally better with lower overlap, while some
methods are helped whereas others are harmed by universal spreaders. Overall
patterns regarding cluster sizes and cluster specific area widths are not that
striking.

The performance of hierarchical clustering is generally rather bad, and
clearly dominated by the other methods. Average linkage performs somewhat
better than the other two methods for low overlap and without universal
spreaders, but it is still worse than all the non-hierarchical methods; in
fact the interaction plots could be interpreted as showing two clusters of
methods, namely the hierarchical methods achieving low ARIs, and most other
methods achieving better ARIs with K-modes and pdfs2 in between. We explore
the reasons for this in Section 5.2.

K-modes performs a bit better than the hierarchical methods, but worse than
the other methods. It becomes mostly better with universal spreaders. In
particular, the use of the simple matching distance seems to be a worse choice
for this kind of data than Jaccard for PAM.

The latent class analysis belongs to the cluster of better methods, but it
cannot compete with K-means and the Gaussian mixture after MDS, as can be seen
in particular in Figure 2(d). It does absolutely and relatively better with
low overlap. PAM is among the best methods with low overlap, but deteriorates
even stronger than the latent class analysis with more overlap.

Out of the different variants of K-means and the Gaussian mixture model,
regarding the overall average ARI those based on a 3-dimensional classical MDS
(kc3, gc3) perform best, as can be seen in Figures 2(c), 2(d). There is
considerable variation over the different experimental factor levels though.
Figure 3 shows that gc2 is excellent with low overlap and without universal
spreaders. kc3 performs much better with universal spreaders and does not
stand out without them. Using the ratio (smacof) MDS is mostly worse than
using classical MDS. Particularly it does not work well with the Gaussian
mixture model and 3-dimensional MDS output with low and medium overlap and no
universal spreaders, although universal spreaders help that method,
particularly with medium overlap. Unequal sizes of cluster specific areas
generally favour the classical MDS. The combination of medium overlap and
universal spreaders generally favours the 3-dimensional MDS strongly, which
otherwise has a small but not so substantial advantage. The K-means methods
are overall competitive, but relatively weaker compared to the Gaussian
mixture model for unequal cluster sizes. High overlap favours the K-means
methods. pdfCluster is overall not quite as good as K-means and the Gaussian
mixture. Overall it has a somewhat similar performance to the latent class
analysis, with pdfc3 doing mostly better and the other three versions doing
mostly worse. The impossibility to fix the number of clusters for pdfCluster
puts it at a slight disadvantage, as the correct number of clusters was fixed
for the other methods. pdfCluster mostly found the correct number of clusters,
but occasionally (particularly with the smacof-MDS) it would put everything
together in a single cluster.

Figure 3 gives some information about the variability of the results. pam and
gc3 have most stable results for low overlap (without universal spreaders also
gc2), whereas other methods, gs3 in particular, have substantial variation.
Otherwise stability can mostly be found among bad results (e.g., the linkage
methods).

Figure 3: Boxplots of the ARIs achieved by the clustering methods, grouped by
level of overlap and presence of a cluster of universal spreaders (sea-green
boxplots: no $u.s.$ group; golden boxplots: with $u.s.$ group).

###  5.2 More detailed insight

Figure 4: Heatmap from a simulation of scenario 24: high overlap, with $u.s.$
group, various sizes and areas. ARI$=0.3538$. Rows sorted by LCA clustering.
True cluster membership on the left (blue=1, red=2, green=3, gold=$u.s.$),
cluster specific areas on top (light-blue=1, orange=2, dark red=3), LCA
cluster membership on the right (numbers).

In order to explore behavior of some methods in more detail and to understand
their weaknesses, we use some exemplary visualizations.

For examining the performance of latent class analysis, we consider a
simulation from scenario 24 of Table 1 with high overlap and different cluster
sizes. The heatmap in Figure 4 shows the four cluster latent class solution
indirectly by the order of rows; colors on the side denote the true clusters,
colors on top denote the cluster specific areas. A typical behavior of latent
class analysis with large overlap was to build a cluster with sparse species,
i.e., with small ranges, regardless to which original cluster they belonged.
The smaller true clusters were hard to identify. Latent class analysis treats
presences and absences symmetrically, which leads to a tendency to group
sparse species together; also species with large range are easily grouped
together, which is a good thing with a universal spreader cluster. Therefore
its existence helps the latent class analysis to often achieve a larger ARI.
It can also happen, however, that the universal spreaders cluster is split up
into two clusters (relatively larger and smaller ranges), as happened in
Figure 4. The symmetric treatment of presences and absences is not good for
identifying the proper clusters, at least not with large overlap. These rely
on presences much more than on absences, and therefore the focus of the
Jaccard distance on joint presences rather than joint absences will help.

(a) Ratio smacof 2D mapping colored according to true cluster membership. The
u’s are universal spreaders.

(b) Ratio smacof 2D mapping colored according LCA clustering.

(c) Classical scaling 2D configuration colored according to true cluster
membership. The u’s are universal spreaders.

(d) Classical scaling 2D configuration colored according to LCA clustering.

Figure 5: Ratio smacof MDS and classical scaling, based on the Jaccard’s
distance matrix computed on the data set simulated under scenario 24 (heatmap
in Figure 4), colored by true clustering (left) and LCA clustering (right).

Figure 5 shows two-dimensional MDS maps of these data, highlighting the true
cluster memberships ((a) and (c)) and the latent class clusters ((b) and (d)),
for classical and ratio (smacof) MDS. Although not based on the MDS input, the
latent class solution looks geometrically sensible in terms of both MDS
visualizations, actually more so than the true clusters. However, latent class
clusters 1 and 2 are very heterogeneous and involve large distances, merging
observations from two or more true clusters together. This is in line with the
fact that for large overlap with universal spreaders, starting from the MDS
output, K-means, which prioritizes within-cluster homogeneity over separation
and geometrically visible cluster shapes, does better than the Gaussian
mixture model (Figure 3). The true clusters 1 and 2 are apparently hard to
separate for any method, but they are still more homogeneous than latent class
cluster 2, which holds the majority of both of them.

Figure 5 also serves for understanding why classical scaling performs better
in our experiments. Geometrically, the ratio MDS is not worse at arranging the
true clusters into a sensible shape, but both K-means and the Gaussian mixture
model are connected to multivariate Gaussian distributions, which are
characterised by linear relations between the variables, whereas the ratio MDS
produces nonlinear cluster shapes and boundaries between true clusters. The
ratio MDS in itself does not fail, but it does not serve as well as input for
the Euclidean cluster methods. The classical MDS also separates the universal
spreaders on one side of the plot, whereas in the ratio MDS they are, although
also nicely separated, surrounded by observations from the proper clusters.
Presence of universal spreaders helped some methods to achieve better results
simply because they are often more easily “visible” as clusters to the
methods.

(a) Single linkage clustering.
ARI$=-0.0039$.

(b) Complete linkage clustering.
ARI$=-0.0051$.

(c) Average linkage clustering. ARI$=0.0004$.

Figure 6: Heatmaps of a simulation of scenario 22 ($\omega=0.05$, with $u.s.$
cluster, various cluster sizes and areas), sorted according to the three
hierarchical clustering solutions. Colours on the left side indicate true
cluster membership of the species, while those on top refer to the cluster
specific area of the cells.

Regarding the weak performance of the hierarchical methods, Figure 6 shows
three heatmaps from data from scenario 22, with low overlap, and universal
spreaders, and clusters have different sizes and different widths of their
specific areas. These are shown together with the dendrograms form the three
hierarchical clusterings.

Single linkage (Figure 6(a)) clearly suffers from its well known chaining
issue [14, Chapter 4]. Fixing the number of clusters at 4, all observations
are merged, and only a few very sparse species with a large Jaccard distance
to all the other species are isolated as one-point clusters. Looking at the
colors indicating the true clusters, the dendrogram as a whole is not perfect,
but to some extent in line with the true clusters. A single clustering
obtained by cutting the dendrogram at a certain height would however require a
very large number of clusters to have at least the bigger ones of them in some
agreement with the true clusters. The same phenomenon occurred with most other
datasets, causing ARIs generally close to zero.

Complete linkage is not affected by chaining, but Figure 6(b) shows a
different problem. The Jaccard distance assigns its maximum value 1 to any two
species that do not have joint presences. As occurs often in real data, the
simulated data have several species with a small range that easily have no
overlap. Complete linkage can’t therefore join them in the same cluster, and
they are still separated at the top level of the dendrogram. The dendrogram as
a whole can be seen as even better in line with the true clustering, but in
order to reflect this in a clustering with a fixed number of clusters, the
method would need to integrate some sparse species with some clusters that
include species to which they have the maximum Jaccard distance, and complete
linkage cannot do that.

Unfortunately, average linkage (Figure 6(c)) cannot solve both of these
problems at the same time, but is rather also affected by chaining, isolating
some sparse species when cutting for lower numbers of clusters. Once more,
some useful agreement with the true clustering can be found in the dendrogram
as a whole, but with a fixed number of clusters 4, the ARI is about zero, and
a too large number of clusters would be needed to achieve a substantially
better value.

The best ARIs that could have been achieved by optimal cutting are 0.427 for
single linkage (181 clusters), 0.469 for complete linkage (35 clusters), and
0.522 for average linkage (81 clusters). These are still lower than those for
the MDS-based methods, latent class analysis, and PAM, largely due to the many
small clusters involved. Still it shows that the results of the hierarchical
methods are ultimately also connected to the true clusterings, though this
requires looking deeper into the dendrogram.

A major lesson to learn from these insights regards the role of the Jaccard
distance. On one hand, it is advantageous by treating presences and absences
in an asymmetric manner. The latent class analysis suffers from treating them
symmetrically, compared to the MDS-based methods. On the other hand, it also
creates problems, particularly for the hierarchical methods, by assigning a
maximum distance to pairs of sparse species that have no cell in common. [17,
18] have preferred the Kulczynski dissimilarity, which is not a metric, but
this also assigns maximum dissimilarity if there is no presence in common,
although it may help by assessing sparse species as closer to species with a
larger range of which they are a subset. [19] propose to involve geographical
distances into the distance computation, and this may indeed improve matters,
if only in a way that depends strongly on the geography of the specific real
data, which will make setting up general models for simulation a more
difficult task.

##  6 Conclusions

We have run a simulation study in order to evaluate the clustering
performance, on multivariate presence-absence data, of a combination of
multidimensional scaling and clustering methods for quantitative data
(K-means, Gaussian mixture models, and pdfCluster), compared with methods that
either operate directly on the data (such as the latent class analysis), or on
a Jaccard or simple matching distance matrix computed from the data
(hierarchical methods, PAM, and K-modes). Two different MDS methods and MDS
outputs of dimensionality 2 and 3 were involved.

The results suggest that the MDS-based techniques can be a valuable tool to
cluster such data, returning by and large better ARI values than the
competitors. The hierarchical methods did particularly badly, mainly due to
the fact that cutting the dendrogram at the correct number of clusters proved
inadequate; however it would be hard to repair this problem by any automatic
method for deciding the number of clusters, because the required number of
clusters for achieving a better ARI would be very large, and therefore
potentially undesirable by researchers who want a simple and interpretable
clustering solution. K-modes worked slightly better than the hierarchical
methods.

Latent class analysis proved to be competitive when the data features were not
demanding, but its performance deteriorated faster than that of MDS-based
methods in more complicated setups. Similarly, PAM performed very well for low
overlap but markedly worse, also in relative terms, with higher overlap.

For the MDS-based methods, classical MDS did a better job than smacof’s ratio
MDS, due to the fact that classical MDS arranged the data in such a way that
clusters could be separated more linearly, which makes the job of K-means and
the Gaussian mixture model easier. However, without any linearity requirement,
pdfCluster shared this behaviour. Cluster methods such as spectral clustering
[27] that do not involve linearity could also be tried with ratio MD. 3
dimensions worked overall slightly better than 2, but there is no reason to
believe that increasing the dimensionality above 3 will improve results
substantially. Comparing Gaussian mixture models and K-means, there is no
clear winner. Both of these approaches worked well in some situations, and
both worked mostly better than pdfCluster; the fact that the latter estimates
the number of clusters automatically may have affected its results though.

Regarding characteristics of the data, unsurprisingly, increasing levels of
overlap among clusters turned out to be the most important factor, making
separation of the true clusters more difficult for all methods.

Interesting extensions for future research could be involving other distance
measures and clustering methods, with possible alternatives suggested above,
and in particular involving methods to estimate the number of clusters,
although this is a hard task due to the large number of possible combinations
between clustering methods for fixed $K$ and methods to choose $K$.

## References

  * [1] Anderlucci, L., Hennig, C. (2014). The clustering of categorical data: A comparison of a model-based and a distance-based approach. Commun. Stat. A Theor. 43, 704–721.
  * [2] Azzalini A, Torelli N (2007) Clustering Via Nonparametric Density Estimation. Stat. Comput. 17, 71–80.
  * [3] Azzalini, A., Menardi, G (2016) Density-based clustering with non-continuous data. Comput. Stat. 31, 771–798.
  * [4] Bartholomew, D. J. (2008). Analysis of Multivariate Social Science Data. CRC Press, Boca Raton.
  * [5] Bartholomew, D. J., Moustaki, I., Knott, M. (2011). Latent Variable Models and Factor Analysis: A Unified Approach. Wiley, New York.
  * [6] Borg, I., Groenen, P. (2005). Modern Multidimensional Scaling: Theory and Applications. Springer, Berlin.
  * [7] Borg, I., Groenen, P., Mair, P. (2013). Applied Multidimensional Scaling. Springer, Berlin.
  * [8] Casagranda, M. D., Taher, L., Szumik, C. A. (2012). Endemicity analysis, parsimony and biotic elements: A formal comparison using hypothetical distributions. Cladistics, 28, 645–654.
  * [9] Cheetham, A., Hazel, J. (1969). Binary (presence-absence) similarity coefficients. J. Paleontol. 43, 1130–1136.
  * [10] Chen, Y. (2013). Biotic element analysis of reptiles of China: A test of vicariance model. Curr. Zool. 59, 449–457.
  * [11] de Leeuw, J., Heiser, W. (1982). Theory of Multidimensional Scaling. In Handbook of Statistics Vol. 2, ed. by P. R. Krishnaiah and I. N. Kanal. North Holland, Amsterdam, pp. 285-316.
  * [12] Dennis, R. L. H., Williams, W. R., Shreeve, T. G. (1998). Faunal structures among European butterflies: Evolutionary implications of bias for geography, endemism and taxonomic affiliation. Ecography, 21, 181–203.
  * [13] Desarbo, W., Howard, D., Jedidi, K. (1991). Multiclus: A new method for simultaneously performing multidimensional scaling and cluster analysis. Psychometrika, 56, 121–136.
  * [14] Everitt, B., Landau, S., Leese, M., Stahl, D. (2011). Cluster analysis. Wiley, New York.
  * [15] Halkidi, M., Vazirgiannis, M., Hennig, C. (2015). Method-independent indices for cluster validation and estimating the number of clusters, in Handbook of Cluster Analysis ed. by C. Hennig, M. Meila, F. Murtagh, R. Rocci, R., CRC Press, Boca Raton, pp. 595-618.
  * [16] Hausdorf, B. (2002). Units in biogeography. Syst. Biol., 51, 648–652.
  * [17] Hausdorf, B., Hennig, C. (2003). Biotic element analysis in biogeography. Syst. Biol., 52, 717–723.
  * [18] Hennig, C., Hausdorf, B. (2004). Distance-based parametric bootstrap tests for clustering of species ranges. Comp. Stat. Data An., 45, 875–895.
  * [19] Hennig, C. Hausdorf, B. (2006). Design of dissimilarity measures: A new dissimilarity between species distribution areas. In Data Science and Classification, ed. by V. Batagelj, H. Bock, A. Ferligoj, A. Ziberna, Springer, Berlin, pp. 29-37.
  * [20] Huang, Z. (1997). A Fast Clustering Algorithm to Cluster Very Large Categorical Data Sets in Data Mining. In KDD: Techniques and Applications, ed. by H. Lu, H. Motoda and H. Luu, World Scientific, pp. 21-34.
  * [21] Hubert, L., Arabie, P. (1985). Comparing partitions. J. Classif. 2, 193–218.
  * [22] Jaccard, P. (1901): Distribution de la florine alpine dans la Bassin de Dranses et dans quelques regiones voisines. B. Soc. Vaudoise.Sci. Nat., 37, 241–272.
  * [23] Jain, A. K. (2010). Data clustering: 50 years beyond k-means. Pattern Recogn. Lett. 31, 651–666.
  * [24] Kaufman, L., Rousseeuw, P.J. (1990) Finding Groups in Data: An Introduction to Cluster Analysis. Wiley, New York.
  * [25] Kruskal, J. (1977). The relationship between multidimensional scaling and clustering, in Classification and Clustering, ed. by J. Van Ryzin, Academic Press, New York, pp. 17-44.
  * [26] McCutcheon, A. L., Hagenaars, J. A. (2002). Applied Latent Class Analysis. Cambridge University Press, Cambridge.
  * [27] Ng, A. Y., Jordan, M. I., Weiss, Y. (2002). On spectral clustering: analysis and an algorithm, in Advances in Neural Information Processing Systems 14 ed. by T. G. Dietterich, S. Becker, Z. Ghahramani, MIT Press, Cambridge, MA, pp. 849-856.
  * [28] Oh, M.-S., & Raftery, A. E. (2007). Model-based clustering with dissimilarities: A Bayesian approach. J. Comput. Graph. Stat. 16, 559–585.
  * [29] Schwarz, G. (1978). Estimating the Dimension of a Model. The Ann. Stat. 6, 461–464.
  * [30] Scrucca L., Fop M., Murphy T. B., Raftery A. E. (2016). mclust 5: clustering, classification and density estimation using Gaussian finite mixture models. R J., 8, 289–317.
  * [31] Shi, G. R. (1993). Multivariate data analysis in palaeoecology and palaeobiogeography - a review. Palaeogeogr. Palaeocl. 105, 199–234.
  * [32] Torgerson, W. S. (1958). Theory and Methods of Scaling. Wiley, New York.
  * [33] Ulrich, W., Gotelli, N.J. (2012). A null model algorithm for presence–absence matrices based on proportional resampling. Ecol. Model. 244, 20–27.
  * [34] Van Mechelen, I., Boulesteix, A.-L., Dangl, R., Dean, N., Guyon, I., Hennig, C., Leisch, F., Steinley, D.: Benchmarking in cluster analysis: A white paper. arXiv (2018), https://arxiv.org/abs/1809.10496
  * [35] Vavrek, M. J. (2016). A comparison of clustering methods for biogeography with fossil datasets. PeerJ, 4, e1720.
  * [36] Vermunt, J. K., Magidson, J. (2002). Latent class cluster analysis, in Applied Latent Class Analysis ed. by A. L. McCutcheon, J. A. Hagenaars. Cambridge University Press, Cambridge, pp.89-106.