---
title: "Introduction to the evclust package"
author: "Thierry Denoeux"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
bibliography: tdenoeux.bib
vignette: >
  %\VignetteIndexEntry{Introduction to the evclust package}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
The package `evclust` contains methods for *evidential clustering*. In evidential clustering, cluster membership  uncertainty is represented by Dempster-Shafer mass functions. The user is invited to read the papers cited in this vignette to get familiar with the main concepts underlying evidential clustering. These papers can be downloaded from the author's web site, at <https://www.hds.utc.fr/~tdenoeux>. Here, we provide a guided tour of the main functions in the `evclust` package. It consists in six main functions, implementing six different evidential clustering algorithms:

* `ecm`: Evidential c-means algorithm [@masson08]
* `recm`: Relational Evidential c-means algorithm [@masson09a]
* `kevclus`: $k$-EVCLUS algorithm [@denoeux04b; @denoeux16a]
* `EkNNclus`: E$k$-NNclus algorithm [@denoeux15]
* `cecm`: Constrained Evidential c-means algorithm [@antoine12]
* `bootclus`: Generating credal partitions by bootstrapping Gaussian mixture models [@denoeux19]

You first need to install this package:

```{r, message=FALSE}
library("evclust")
```


## Evidential c-means (ECM) algorithm

The Evidential $c$-means (ECM) algorithm [@masson08] is a $c$-means-like algorithm that minimizes a cost function by searching alternatively  the space of prototypes and the space of credal partitions. Unlike  the hard and fuzzy $c$-means algorithms, ECM associates a prototype not only to clusters, but also to sets of clusters. The prototype associated to a set of clusters is defined as the barycenter of the prototypes of each single cluster in the set. The cost function to be minimized insures that objects close to a prototype have a large mass assigned to the corresponding set of clusters. 

Consider, for instance, the `fourclass` data. This dataset consists in four clusters in a two-dimensional space. 

```{r}
data(fourclass)
x<-fourclass[,1:2]
y<-fourclass[,3]
plot(x[,1],x[,2],pch=y,col=y,xlab=expression(x[1]),ylab=expression(x[2]))
```

We can run ECM with `c=4` clusters on this data as follows:

```{r}
clus<-ecm(x,c=4,type='full',alpha=1,beta=2,delta=sqrt(20),disp=FALSE)
```

The option `type='full'` is actually the option. It implies that mass functions in the credal partition will have $2^c$ focal sets. You can get basic information about the credal partition using the method `summary`:

```{r}
summary(clus)
```

We can restrict the focals set to pairs by changing the `type` argument.

```{r}
clus<-ecm(x,c=4,type='pairs',alpha=1,beta=2,delta=sqrt(20),disp=FALSE)
summary(clus)
```

A plot of the credal partition can be generated as follows:

```{r, fig.width=6, fig.height=6}
plot(clus,x,mfrow=c(2,2),ytrue=y,Approx=2)
```

In this plot, the lower and upper approximations of each cluster are plotted as solid and interrupted lines, respectively. Since we selected `Approx=2`, the lower and upper approximations are defined as follows. Let $A_i$ be the focal set of mass function $m_i$ with the largest mass. The lower approximation of cluster $\omega_k$  is  the set of objects such that $A_i=\{\omega_k\}$, while the upper approximation is  the set of objects such that $\omega_k \in A_i$. The outliers are the objects such that $A_i=\emptyset$. They are displayed as circles in the figure above.

## RECM

The ECM algorithm can only be used with attribute data. The Relational Evidential c-means algorithm (RECM) [@masson09a] is a version of ECM that can handle dissimilarity data, i.e., data that consist in a matrix on dissimilarities between $n$ objects. Consider, for instance, the Proteins dataset. It consists of a dissimilarity matrix derived from the structural comparison of 213 protein sequences. Each of these proteins is known to belong to one of four classes of globins: hemoglobin-$\alpha$ (HA), hemoglobin-$\beta$ (HB), myoglobin (M) and heterogeneous globins (G). We can display these data using Multidimensional Scaling:

```{r}
data(protein)
z<- cmdscale(protein$D,k=2)
plot(z[,1],z[,2],xlab='axis 1',ylab='axis 2')
```

We can see that the data points are grouped in four clusters. We can run RECM on this dataset at follows,

```{r, fig.width=6, fig.height=6}
clus <- recm(D=protein$D,c=4,type='pairs',alpha=0.2,beta=1.1,delta2=20,disp=FALSE)
summary(clus)
plot(clus,X=z,mfrow=c(2,2),ytrue=protein$y)
```

## k-EVCLUS

The RECM algorithm requires to store the whole dissimilarity matrix. Consequently, it is not suitable to handle large dissimilarity datasets. The k-EVCLUS algorithm, recently introduced in [@denoeux16a], overcomes this limitation. It is an improved version of the EVCLUS algorithm introduced in [@denoeux04b]. The EVCLUS algorithm applies  ideas from Multidimensional Scaling to clustering: given a dissimilarity matrix, it finds a credal partition such that the degrees of conflict between mass functions match the dissimilarities, dissimilar objects being represented by highly conflicting mass functions; this is achieved by iteratively minimizing a stress function. The k-EVCLUS algorithm uses a more efficient optimization procedure, and uses only the dissimilarities between each object and $k$ randomly chosen other objects. As a result, the complexity is reduced from quadratic to linear. 

As an example, let us generate a dataset from the same distribution as the `fourclass` dataset, but with $n=1000$ objects.

```{r}
n<-1000
x<-matrix(rt(2*n,df=5),n,2)
y<-c(rep(1,n/4),rep(2,n/4),rep(3,n/4),rep(4,n/4))
x[(y==2)|(y==4),2]<- 5 + x[(y==2)|(y==4),2]
x[(y==3)|(y==4),1]<- 5 + x[(y==3)|(y==4),1]

plot(x[,1],x[,2],pch=y,col=y,xlab=expression(x[1]),ylab=expression(x[2]))
```

There are $1000\times 999/2=499,500$ pairwises distances for these data. However, by choosing $k=50$, we can use only $1000\times 50=50,000$ distances, i.e., only 10 percent of the original distances to generate a credal partition.

```{r k-EVCLUS, eval=TRUE, include=TRUE}
clus<-kevclus(x=x,c=4,k=50,disp=FALSE)
summary(clus)
```

For a partition created by `kevclus`, the `plot.credpart` functions generates two plots: the usual plot of the credal partition in the attribute space, and the *Shepard diagram*. This is a plot of the degrees of conflict vs. the transformed dissimilarities. If the argument `X` is not supplied,  only the Shepard diagram is plotted.

```{r, eval=TRUE}
plot(clus,X=x,ytrue=y)
```

## EKNN-clus

The clustering methods described above require the number of cluters to be fixed in advance. In contrast, EK-NNclus [@denoeux15] is another evidential clustering algorithm that automatically determines the number of clusters. Starting from an initial partition, E$K$-NNclus iteratively reassigns objects to clusters using the E$K$-NN rule, until a stable partition is obtained. After convergence, the cluster membership of each object is described by a Dempster-Shafer mass function assigning a mass to each cluster and to the whole set of clusters. The mass assigned to the set of clusters can be used to identify outliers. The method can be implemented in a competitive Hopfield neural network, whose energy function is  related to the plausibility of the partition. The procedure can thus be seen as searching for the most plausible partition of the data. The  E$K$-NNclus algorithm  depends on two parameters, the number $K$ of neighbors and a scale parameter $q$, which can be fixed using simple heuristics.

As an example, let us apply E$K$-NNclus to the `s2` dataset. This dataset contains 5000 two-dimensional vectors grouped in 15 Gaussian clusters.

```{r, fig.width=6, fig.height=6}
data(s2)
plot(s2[,1],s2[,2],xlab=expression(x[1]),ylab=expression(x[2]),pch='.')
```

To apply E$K$-NNclus, we need to fix $k$ and $q$. Results presented [@denoeux15] suggest that a value of $K$ of the order of two or there times $\sqrt{n}$ and $q$ greater than 0.5 are suitable. For instance, let use set $K=200$ and $q=0.9$. The algorithm is initialized with a random partition in 500 clusters.  It is run `ntrials=5` times, and the result with the lowest value of the energy function is kept.

```{r, fig.width=6, fig.height=6}
y0<-sample(500,size=5000,replace=TRUE)
clus<- EkNNclus(s2,K=200,y0=y0,ntrials=5,q=0.9,disp=FALSE)
## Summary of the partition
summary(clus)
## Plot of the partition
plot(clus,X=s2,Outliers=FALSE)
```

We can see that E$K$-NNclus generally finds 15 or 14 clusters for these data. In constrast with the previous algorithms, E$K$-NNclus generates a credal partition with normalized mass functions. Outliers tend to have a large mass $m(\Omega)$ on the frame of discernment, but the credal partition is otherwise usually very close to a hard partition.

## Generating a credal partition with a large number of clusters

If no restriction is imposed on the focal sets, the number of parameters to be estimated in evidential clustering grows exponentially with the number $c$ of clusters, which makes it intractable unless $c$ is small. If we allow masses to be assigned to pairs of clusters, as suggested in [\@denoeux04b] and [\@masson08], the number of focal sets becomes proportional to $c^2$, which is manageable for moderate values of $c$ (say, until 10), but still impractical for larger $n$. It is clear, however, that only a few  pairs of contiguous clusters will be assigned some mass during the learning process. To determine which pairs of clusters can potentially become focal sets, we can use the following three-step approach [@denoeux16a]:

1. In the first step, a clustering algorithm (`ecm`, `recm` or `kevclus`) is run in the basic configuration, with focal sets of cardinalities 0, 1 and (optionally) $c$. A credal partition $M_0$ is obtained. 
2. The similarity between each pair of clusters $(\omega_j,\omega_\ell)$ is analyzed, and we determine the set $P_K$ of pairs $\{\omega_j,\omega_\ell\}$ that are mutual $K$ nearest neighbors.
3. The clustering algorithm is run again, starting from the previous credal partition $M_0$, and  adding as focal sets the pairs in $P_K$.

As an example, consider again the `s2` dataset with the `ecm` algorithm. We first construct a credal partition with $c=15$ clusters and $f=16$ focal sets (the empty set and the sigletons):

```{r, fig.width=6, fig.height=6}
data(s2)
clus<-ecm(x=s2,c=15,type='simple',Omega=FALSE,delta=1,disp=FALSE)
plot(x=clus,X=s2,Outliers = TRUE)
```

We then search for pairs of clusters that are $K=2$ mutual nearest neighbors:

```{r}
P<-createPairs(clus,k=2)
P$pairs
```
Finally, we run again `ecm`, with the pairs of clusters found in the previous step:

```{r, fig.width=6, fig.height=6}
clus1<-ecm(x=s2,c=15,type='pairs',Omega=FALSE,pairs=P$pairs,g0=clus$g,delta=1,disp=FALSE)
plot(x=clus1,X=s2,Outliers = TRUE,Approx=2)
```

This time, we can see that the outer approximations of neighboring clusters overlap, reflecting uncertainty in the partition.

## Evidential clustering with pairwise constraints

The clustering task can be made easier if we have some prior knowledge of the underlying partition. Such knowledge can take the form of *must-link* (ML) constraints (pairs of objects that are known to belong to the same clusters) and *cannot-link* (CL) constraints (pairs of objects that are known to belong to different clusters). The CECM algorithm [@antoine12] is a variant of ECM that makes it possible to handle such constraints. 

Consider, for instance, the Four-class dataset. It is actually composed of four clusters, but let us assume that the partition of interest has only two groups, with clusters 1-2 and 3-4 in the first group. We can display the corresponding binary partition:

```{r}
data(fourclass)
x<-as.matrix(fourclass[,1:2])
y<-as.vector(fourclass[,3])
## We group classes 1-2 and 3-4
y<-2-(y<=2)
plot(x[,1],x[,2],pch=y,col=y)
```

Depending on the initial conditions, `ecm` may converge to this partition, or to a partition close to the following one:

```{r}
y1<-((fourclass[,3]==1) | (fourclass[,3]==3))+1
plot(x[,1],x[,2],pch=y1,col=y1)
```

To guide the clustering algorithm to the desired partition, we can provide ML and CL constraints. They can be randomly generated using the `create_MLCL` function:

```{r}
const<-create_MLCL(y,nbConst=30)
```

Here, we generated 30 constraints, of ML or VL type. These constraints can be passed to  function `cecm`. Other important arguments of `cecm` are the coefficient `bal`, which determines the tradeoff between the stress and penalty terms in the cost function, and `distance`, which determines the distance function used. The default (`distance=0`) is the Euclidean function. By setting `distance=1`, we get generalized Euclidean (Mahalanobis) distances for each cluster.

```{r, message=FALSE}
clus<-cecm(x,c=2,type='full',ntrials=5,ML=const$ML,CL=const$CL,delta=10,bal=0.5,distance=1,disp=FALSE)
plot(clus,X=x,ytrue=y,Outliers=TRUE,Approx=1)
```

More sophisticated active learning stategies for selecting the constraints are described in [@antoine12].

## BootClus

This is the most recent algorithm for generating a credal partition. It works by bootstrapping mixture models [@denoeux19]. The current implementation uses the `mclust` package and is thus restricted to Gaussian mixtures, but the approach can be applied to  any mixture model. The idea is to compute bootstrap percentile intervals on each pairwise probability $P_{ij}$ that objects $i$ and $j$ belong to the same cluster. We then search for a credal partition such that the belief and plausibilities that each pair of objects $i$ and $j$ belong to the same cluster approximate the bounds of the confidence intervals in the least squares sense.

Let us illustrate this approach with the Faithful geyser data that are part of package `mclust`:

```{r}
data("faithful")
X<-faithful
```

We run function `bootclus` with `B=100` bootstrap samples to shorten execution time (better results are obtained with at least 500 bootstrap samples). The argument `alpha` specifies the levels of the quantiles used to compute the confidence intervals; here, by choosing 5\% and 95\% quantiles, we obtain equal-tail 90\% confidence inetrvals. The argument `G=3` is passed to function `Mclust`called inside `bootclus` though the list `param`, to obtain a partition with three clusters:

```{r}
param=list(G=3)
res.faithful<-bootclus(X,alpha=c(0.05,0.95),B=100,param=param)
```

Function `bootclus` returns the credal partition `Clus`, an object of class `credpart`:

```{r}
plot(res.faithful$Clus,X)
```

The confidence inetrvals as well as the approximating belief-plausibility intervals are also returned by function `bootclus` as arrays `CI` and `BelPl`, respectively. the following graphs show the pairwise belief and and plausibility degrees vs. the lower and upper bounds of the confidence intervals, which allows us to assess the quality of the approximation:

```{r,fig.width=5, fig.height=5}
plot(as.dist(res.faithful$CI[1,,]),as.dist(res.faithful$BelPl[1,,]),pch=".",xlab="",ylab="")
title(xlab="Lower bound of 90% CI",ylab="Belief",line=2.2,cex.lab=1.7)
abline(0,1)
plot(as.dist(res.faithful$CI[2,,]),as.dist(res.faithful$BelPl[2,,]),pch=".",xlab="",ylab="")
title(xlab="Upper bound of 90% CI",ylab="Plausibility",line=2.2,cex.lab=1.7)
abline(0,1)
```



## References


