cours_clustering.Rmd

---
title: "Méthodes de Partionnement et d'apprentissage non supervisé"
subtitle: Classification Hiérarchique et Kmeans
author: "Anne Badel, Frédéric Guyon & Jacques van Helden"
date: "`r Sys.Date()`"
output:
  slidy_presentation:
    highlight: default
    incremental: no
    smart: no
    slide_level: 1
    self_contained: yes
    fig_caption: no
    fig_height: 5
    fig_width: 5
    keep_md: yes
    smaller: yes
    theme: cerulean
    toc: yes
    widescreen: yes
  beamer_presentation:
    colortheme: dolphin
    fig_caption: no
    fig_height: 5
    fig_width: 5
    fonttheme: structurebold
    highlight: default
    incremental: no
    keep_tex: no
    slide_level: 1
    theme: Montpellier
    toc: yes
  html_document:
    fig_caption: yes
    highlight: zenburn
    self_contained: no
    theme: cerulean
    toc: yes
    toc_depth: 3
    toc_float: yes
  ioslides_presentation:
    highlight: zenburn
    incremental: no
  pdf_document:
    fig_caption: yes
    highlight: zenburn
    toc: yes
    toc_depth: 3
  revealjs::revealjs_presentation:
    theme: night
    transition: none
    self_contained: true
    slide_level: 1
    css: ../slides.css
font-import: http://fonts.googleapis.com/css?family=Risque
font-family: Garamond
transition: linear
---


```{r include=FALSE, echo=FALSE, eval=TRUE}
library(knitr)
library(kableExtra)
library(FactoMineR)
library(clues)
library(RColorBrewer)
# library(formattable)

options(width = 300)

# options(encoding = 'UTF-8')

knitr::opts_chunk$set(
  fig.width = 5, fig.height = 5, 
  fig.path = 'figures/irisDeFisher_',
  fig.align = "center", 
  size = "tiny", 
  echo = TRUE, eval = TRUE, 
  warning = FALSE, message = FALSE, 
  results = TRUE, comment = "")

options(scipen = 12) ## Max number of digits for non-scientific notation
```

```{r  echo=F}
mes.iris <- iris[, 1:4]
mes.iris.scaled <- data.frame(scale(mes.iris, T, T))
nb.iris=nrow(mes.iris)
nb.var=ncol(mes.iris)
```


# Partitionnement et apprentissage

- On a une **représentation** des données

    + sous forme de valeurs réelles=vecteur de 
    + sous forme de catégories

- **Clustering**: on cherche a priori des groupes dans les données

- **Apprentissage**:

    + on connaît le partitionnement sur un jeu de données
    + on cherche le groupe (la classe) de nouvelles données

# Partionnement = Clustering

```{r are_there_clusters, echo=FALSE, out.width="50%", fig.cap="Y a-t-il des groupes ?"}
include_graphics(path = "img/figure1.png")
```


# Partionnement = Clustering

```{r yes_there_are_4_groups, echo=FALSE, out.width="50%", fig.cap="Oui, 4 groupes. "}
include_graphics(path = "img/figure2.png")
```

# Apprentissage

```{r 2_groups, echo=FALSE, out.width="50%", fig.cap="2 groupes. "}
include_graphics(path = "img/figure3.png")
```

# Apprentissage: Séparation linéaire

```{r 2_groups_bis, echo=FALSE, out.width="50%", fig.cap="2 groupes. "}
include_graphics(path = "img/figure4.png")
```


# Méthodes
Trois grands principes de méthodes basées sur:

- La géométrie
- Les probabilités (statistique)
- Les graphes

En fait, trois façons de voir les mêmes algorithmes

# Géométrie et distances

<div class="columns-2">

On considère les données comme des points de $R^n$ (*)


- géométrie donnée par distances
- distances = dissimilaritées imposées par le problème
- dissimilarités $\longrightarrow$ permettent visualisation de l'ensemble des points
- Détermination visuelle des groupes

(*) Espace Euclidien à $n$ dimensions, où 

- chaque dimension représente une des variables observées;
- un individu est décrit comme un vecteur à $n$ valeurs, qui correspond à un point dans cet espace. 

</div>


# Les données

[Ces données sont un classique des méthodes d'apprentissage](https://onlinelibrary.wiley.com/doi/epdf/10.1111/j.1469-1809.1936.tb02137.x)

Dans un premier temps, regardons les données.

```{r iris_data}
dim(mes.iris)
head(mes.iris)
```

# Les variables

```{r iris_summary}
summary(mes.iris)
```


# Visualisation des données

On peut ensuite essayer de visualiser les données
  
```{r plot_iris, fig.width=7, fig.height=6, out.width="60%"}
plot(mes.iris)
```

*****

# Cas d'étude : TCGA Breast Invasive Cancer (BIC)

- Présentation du cas d'étude (Jacques van Helden A COMPLETER)

# TP : analyse de données d'expression

- TP clustering : 
[[html](TP_clustering.html)]
[[pdf](TP_clustering.pdf)]
[[Rmd](https://raw.githubusercontent.com/DU-Bii/module-3-Stat-R/master/seance_4/TP_clustering.Rmd)]

- Première partie : chargement des données

# Géométrie et distances

Sur la base d'une distance (souvent euclidienne)

- Partionnement:

    + Moyennes mobiles ou K-means : séparation optimale des groupes connaissant le nombre de groupes
    + Méthode agglomérative ouhierarchical clustering

- Classification:

    + attribution K plus proches voisins (K Nearest Neighbor)
    + séparation linéaire ou non linéaire

# Distances

Définition d'une distance : fonction positive de deux variables

1. $d(x,y) \ge 0$
2. $d(x,y) = d(y,x)$
3. $d(x,y) = 0 \Longleftrightarrow x = y$
4. **Inégalité triangulaire :** $d(x,z) \le$ d(x,y)+d(y,z)

Si 1,2,3 : dissimilarité

# Distances utilisées dans R

- distance euclidienne ou distance $L_2$:
  $d(x,y)=\sqrt{\sum_i (x_i-y_i)^2}$
  
- distance de manahattan ou distance $L_1$:
  $d(x,y)=\sum_i |x_i-y_i|$
    
- distance du maximum ou L-infinis, $L_\infty$:
  $d(x,y)=\max_i |x_i-y_i|$
    
```{r distances, echo=FALSE, out.width="70%"}
include_graphics(path = "img/distance.png")
```


# Distances utilisées dans R

- distance de Minkowski $l_p$:
$$d(x,y)=\sqrt[p]{\sum_i (|x_i-y_i|^p}$$

- distance de Canberra (x et y valeurs positives):
$$d(x,y)=\sum_i \frac{x_i-y_i}{x_i+y_i}$$

- distance binaire ou distance de Jaccard ou Tanimoto: proportion de propriétés communes

# Autres distances non géométriques (pour information)

Utilisées en bio-informatique:

- Distance de **Hamming**: nombre de remplacements de caractères (substitutions)

- Distance de **Levenshtein**: nombre de substitutions, insertions, deletions entre deux chaînes de caractères

$$d("BONJOUR", "BONSOIR")=2$$

- Distance d'**alignements**: distances de Levenshtein avec poids (par ex. matrices BLOSSUM)

- Distances d'**arbre** (Neighbor Joining)

- Distances **ultra-métriques** (phylogénie UPGMA)

# Distances plus classiques en génomique

Comme vu lors de la séance 3, il existe d'autres mesures de distances :

- **Jaccard** (comparaison d'ensembles): $J_D = \frac{A \cap B}{A \cup B}$

- Distance du $\chi^2$ (comparaison de tableau d'effectifs)

Ne sont pas des distances, mais indices de dissimilarité :

- **Bray-Curtis** (en écologie, comparaison d'abondance d'espèces)
- **Jensen-Shannon** (comparaison de distributions)
  
**Remarque** : lors du TP, sur les données d'expression RNA-seq, nous utiliserons le **coefficient de corrélation de Spearman** et la distance dérivée, $d_c = 1-r$

# Distances entre groupes

- **Single linkage** : élements les plus proches des 2 groupes

$$D(C_1,C_2) = \min_{i \in C_1, j \in C_2} D(x_i, x_j)$$
  
- **Complete linkage** : éléments les plus éloignés des 2 groupes

$$D(C_1,C_2) = \max_{i \in C_1, j \in C_2} D(x_i, x_j)$$
  
- **Group average** : distance moyenne

$$D(C_1,C_2) = \frac{1}{N_1 N_2} \sum_{i \in C_1, j \in C_2} D(x_i, x_j)$$
  
- **Ward**

$d^2(C_i,C_j) = I_{intra}(C_i \cup C_j)-I_{intra}(C_i)-I_{intra}(C_j)$

$D(C_1,C_2) = \sqrt{\frac{N_1N_2}{N_1 + N_2}} \| m_1 -m_2 \|$

# Distances entre groupes

```{r group_distances, echo=FALSE, out.width="90%"}
include_graphics(path = "img/groupes.png")
```


# Les données

[Ces données sont un classique des méthodes d'apprentissage](https://onlinelibrary.wiley.com/doi/epdf/10.1111/j.1469-1809.1936.tb02137.x)

Dans un premier temps, regardons les données

```{r}
dim(mes.iris)
head(mes.iris)
```

---

```{r}
str(mes.iris)
summary(mes.iris)
```


# Visualisation des données

On peut ensuite essayer de visualiser les données

  - par un `plot`
  
```{r plot_4variables, out.width="60%"}
plot(mes.iris)
```

# Visualisation des données - coloration par espèces


```{r plot_4variables_variety, out.width="60%"}
species.colors <- c(setosa = "#BB44DD", virginica = "#AA0044", versicolor = "#4400FF")
plot(mes.iris, col = species.colors[iris$Species], cex = 0.7)
```


# Visualisation des données

- par la fonction `image()`
  
```{r var_image, out.width="60%"}
image(1:nb.var, 1:nb.iris ,t(as.matrix(mes.iris)), xlab = "variables", ylab = "Observation nb", las = 1)
```

# Nettoyage des données (1): données manquantes

Avant de commencer à travailler, il est nécessaire de commencer par vérifier que :

- il n'y a pas de données manquantes

```{r}
sum(is.na(mes.iris))
```

# Nettoyage des données (2) : variables constantes

- aucune  variable n'est constante (aucune variable n'a une variance nulle)
  
```{r}
iris.var <- apply(mes.iris, 2, var)
kable(iris.var, digits = 3, col.names = "Variance")
sum(apply(mes.iris, 2, var) == 0)
```


# Normalisation

Afin de pouvoir considérer que toutes les variables sont à la même échelle, il est parfois nécessaire de normaliser les données.

- soit 

    - en centrant (ramener la moyenne de chaque variable à $0$)
  
```{r}
mes.iris.centre <- scale(mes.iris, center=TRUE, scale=FALSE)
```

- soit 

    - en centrant (ramener la moyenne de chaque variable $0$)
    - et mettant à l'échelle (ramener la variance de chaque variable à $1$)

```{r}
mes.iris.scaled <- scale(mes.iris, center=TRUE, scale=TRUE)
```


# On peut visuellement regarder l'effet de la normalisation :

## par un plot des données 

```{r plot_raw_vs_scaled, out.width="60%", fig.width=12, fig.height=6}
plot(mes.iris, main = "Raw variables")
```

**!** ne pas faire si "grosses" données


----

## ... par une boîte à moustaches (boxplot)

```{r boxplots_raw_vs_norm, out.width="70%", fig.width=8, fig.height=6}
par(mfrow = c(1,2))
par(mar = c(7, 4.1, 4.1, 1.1)) # adapt margin sizes for the labels
boxplot(mes.iris, main = "Raw data", las = 2)
boxplot(mes.iris.scaled, main = "scaled", las = 2)
par(mar = c(5.1, 4.1, 4.1, 2.1)) # Restore original margin sizes
```

----

## ... par une image

```{r image_raw_vs_norm, out.width="80%", fig.width=12, fig.height=6}
par(mfrow=c(1,2))
image(1:nb.var, 1:nb.iris, t(as.matrix(mes.iris)), main="Raw data")
image(1:nb.var, 1:nb.iris, t(as.matrix(mes.iris.scaled)), main="Scaled data")
```


----

## ... par une projection sur une ACP

```{r,  out.width="90%", fig.width=12, fig.height=6}
par(mfrow = c(1,2))
biplot(prcomp(mes.iris), las = 1, cex = 0.7,
       main = "Données non normalisées")
biplot(prcomp(mes.iris, scale = TRUE), las = 1, cex = 0.7,
       main = "Données normalisées")
```

# La matrice de distances

Nous utilisons ici la distance euclidienne. 

```{r echo=F}
mes.iris.scaled <- scale(mes.iris, center=T, scale=F)
iris.acp <- FactoMineR::PCA(mes.iris, graph=F, scale.unit=F)
iris.scaled.acp <- FactoMineR::PCA(mes.iris.scaled, graph=F)
```

```{r}
iris.euc <- dist(mes.iris)
iris.scale.euc <- dist(mes.iris.scaled)
```

```{r heatmap_euclidian_dist, out.width="80%", fig.width=12, fig.height=6}
par(mfrow = c(1,2))
image(t(as.matrix(iris.euc)), main = "Données brutes", las = 1)
image(t(as.matrix(iris.scale.euc)), main = "Données normalisées", las = 1)

par(mfrow = c(1,1))
```

# La classification hiérarchique

## Principe

- **classification hiérarchique** : mettre en évidence des liens hiérachiques entre les individus

    - classification hiérarchique **ascendante** : partir des individus pour arriver à des classes / cluster
    
    - classification hiérarchique **descendante** : partir d'un groupe qu'on subdivise en sous-groupes /clusters jusqu'à arriver à des individus.

---

## Notion importante, cf distances

- ressemblance entre individus = distance
- ressemblance entre groupes d'invidus = critère d'aggrégation

    - lien simple
    - lien complet
    - lien moyen
    - critère de Ward

---

## L'algorithme

### étape 1 :

- départ : n individus = n clusters distincts
- calcul des distances entre tous les individus

    - choix de la métrique à utiliser en fonction du type de données

- regroupement des 2 individus les plus proches => (n-1) clusters

---

## au départ

```{r hclust_initial, echo=FALSE, out.width="90%"}
include_graphics(path = "img/hclust1.png")
```


---

## identification des individus les plus proches

```{r hclust_closest, echo=FALSE, out.width="90%"}
include_graphics(path = "img/hclust2.png")
```


---

## construction du dendrogramme

```{r hclust_dendrogram, echo=FALSE, out.width="90%"}
include_graphics(path = "img/hclust3.png")
```

---

## étape j :

- calcul des dissemblances entre chaque groupe obtenu à l'étape $(j-1)$

- regroupement des deux groupes les plus proches => $(n-j)$ clusters
  
---

## calcul des nouveaux représentants 'BE' et 'CD'

```{r hclust_nvx_representants, echo=FALSE, out.width="90%"}
include_graphics(path = "img/hclust4.png")
```

---

## calcul des distances de l'individu restant 'A' aux points moyens

```{r hclust_update_dist, echo=FALSE, out.width="90%"}
include_graphics(path = "img/hclust5.png")
```


---

## A est plus proche de ...

```{r hclust_closest_a, echo=FALSE, out.width="90%"}
include_graphics(path = "img/hclust6.png")
```

---

## dendrogramme

```{r hclust_dendrogram2, echo=FALSE, out.width="90%"}
include_graphics(path = "img/hclust7.png")
```

---

## pour finir

```{r hclust_finalize, echo=FALSE, out.width="90%"}
include_graphics(path = "img/hclust8.png")
```

---

- à l'étape $(n-1)$, tous les individus sont regroupés dans un même cluster

## dendrogramme final

```{r hclust_dendrogram_final, echo=FALSE, out.width="90%"}
include_graphics(path = "img/hclust9.png")
```

---

## Je ne fais pas attention à ce que je fais ...

```{r dont_care, out.width="70%"}
iris.hclust <- hclust(iris.euc)
plot(iris.hclust, hang = -1, cex = 0.5)
```

... c'est à dire aux options des fonctions `dist()` et `hclust()`

---

## Sur données normalisées

```{r dont_care_norm, out.width="70%"}
iris.scale.hclust <- hclust(iris.scale.euc)
plot(iris.scale.hclust, hang = -1, cex = 0.5)
```

---

```{r dont_care_raw_vs_norm, out.width="80%", fig.width=10, fig.height=7}
par(mfrow = c(2, 1))
plot(iris.hclust, hang = -1, cex = 0.5, main = "Données brutes")
plot(iris.scale.hclust, hang = -1, cex = 0.5, main = "Normalisées")
```

---

## En utilisant une autre métrique

```{r hclust_euclidian_vs_manhattan, out.width="80%", fig.width=10, fig.height=7}
iris.scale.max <- dist(mes.iris.scaled, method = "manhattan")
iris.scale.hclust.max <- hclust(iris.scale.max)
par(mfrow=c(2,1))
plot(iris.scale.hclust, hang=-1, cex=0.5, main = "Euclidian dist")
plot(iris.scale.hclust.max, hang=-1, cex=0.5, main = "Manhattan dist")
```

---

## En utilisant un autre critère d'aggrégation

```{r linkage_rule, out.width="80%", fig.width=10, fig.height=7}
iris.scale.hclust.single <- hclust(iris.scale.euc, method="single")
iris.scale.hclust.ward <- hclust(iris.scale.euc, method="ward.D2")
par(mfrow=c(2,1))
plot(iris.scale.hclust.single, hang=-1, cex=0.5, main = "Single linkage")
plot(iris.scale.hclust.ward, hang=-1, cex=0.5, main = "Ward linkage")
par(mfrow=c(1,1))
```

# Les k-means

Les individus dans le plan

```{r kmeans_initial, echo=FALSE, out.width="80%"}
include_graphics(path = "img/kmeans0.png")
```


# L'algorithme

## étape 1 :

- $k$ centres provisoires tirés au hasard
- $k$ clusters créés à partir des centres en regroupant les individus les plus proches de chaque centre
- obtention de la partition $P_0$ 

---

## Choix des centres provisoires 

```{r kmeans_init_centers, echo=FALSE, out.width="80%"}
include_graphics(path = "img/kmeans1.png")
```


---

## Calcul des distances aux centres provisoires 

```{r kmeans_dist_to_centers, echo=FALSE, out.width="80%"}
include_graphics(path = "img/kmeans2.png")
```

---

## Affectation à un cluster

```{r kmeans_cluster_assignment, echo=FALSE, out.width="80%"}
include_graphics(path = "img/kmeans3.png")
```

---

## Calcul des nouveaux centres de classes

## Etape j :

- construction des centres de gravité des k clusters construits à l’étape $(j-1)$

- $k$ nouveaux clusters créés à partir des nouveaux centres suivant la même règle qu’à l’étape $0$

- obtention de la partition $P_j$

```{r kmeans_update_centers, echo=FALSE, out.width="80%"}
include_graphics(path = "img/kmeans4.png")
```

---

## Fin : 

- l’algorithme converge vers une partition stable

## Arrêt : 

- lorsque la partition reste la même, ou lorsque la variance intra-cluster ne décroit plus, ou lorsque le nombre maximal d’itérations est atteint.

```{r kmeans_stop, echo=FALSE, out.width="80%"}
include_graphics(path = "img/kmeans5.png")
```

---

## Un premier k-means en 5 groupes

```{r calc_kmeans, results = TRUE}
iris.scale.kmeans5 <- kmeans(mes.iris.scaled, center=5)
iris.scale.kmeans5
```

---

```{r cluster_sizes}
iris.scale.kmeans5$cluster
table(iris.scale.kmeans5$cluster)
```

---

## Visualisation des clusters

```{r cluster_viz}
plot(iris.scaled.acp, col.ind = iris.scale.kmeans5$cluster, choix = "ind", cex = 0.5)

```

---

## Combien de clusters ?

Quand une partition est-elle bonne ?
  
- si les individus d’un même cluster sont proches

    - homogénéité maximale à l’intérieur de chaque cluster

- si les individus de 2 clusters différents sont éloignés

    - hétérogénéité maximale entre chaque cluster

---

## Classification hiérarchique

La coupure de l’arbre à un niveau donné construit une partition. la coupure doit se faire :

- après les agrégations correspondant à des valeurs peu élevées de l’indice

- avant les agrégations correspondant à des niveaux élevés de l’indice, qui dissocient les groupes bien distincts dans la population.

#

```{r plot_iris_ward, out.width="95%", fig.width=10, fig.height=5}
plot(iris.scale.hclust.ward, hang=-1, cex=0.5)
```

---

## K-means
```{r}
I.intra = numeric(length=10)
I.intra[1] = kmeans(mes.iris.scaled, centers=2)$totss
for (i in 2:10) {
  kmi <- kmeans(mes.iris.scaled, centers=i)
  I.intra[i] <- kmi$tot.withinss
}
```

# 

```{r}
plot(1:10, I.intra, type="l")
```

#

```{r}
iris.scale.kmeans3 <- kmeans(mes.iris.scaled, center=3)
plot(iris.scaled.acp, col.ind=iris.scale.kmeans3$cluster, choix="ind")
```

# Heatmap

```{r headmapt, out.width="65%", fig.width=8, fig.height=8}
heatmap(mes.iris.scaled, margins = c(7,4), cexCol = 1.4, cexRow = 0.5)
```

---

```{r echo=F}
variete <- iris[,5]
```

```{r heatmap_rowsidecol, out.width="65%", fig.height=8, fig.width=8}
my_group <- as.numeric(as.factor(substr(variete, 1 , 2)))
my_col <- brewer.pal(3, "Set1")[my_group]
heatmap(mes.iris.scaled, RowSideColors = my_col, 
        margins = c(7,4), cexCol = 1.4, cexRow = 0.5)
```

# Comparaison de clustering: Rand Index

Mesure de similarité entre deux clustering

à partir du nombre de fois que les classifications sont d'accord

$$R=\frac{m+s}{t}$$

- m=nombre de paires dans la même classe dans les deux classifications
- s=nombre de paires séparées dans les deux classifications
- t=nombre de paires totales

# Comparaison de clustering: Adjusted Rand Index

$$ ARI=\frac{RI-Expected RI}{Max RI -Expected RI}$$

- ARI=RI normalisé
- Prend en compte la taille des classes
- ARI=1 pour classification identique
- ARI $\simeq$ 0 pour classification aléatoire (peut être <0)
- Adapté pour nombre de classe différent entre les deux classifications
et taille de classe différente

# Comparaison des résultats des deux classifications

- par une table de confusion

```{r results = TRUE}
cluster.kmeans3 <- iris.scale.kmeans3$cluster
cluster.hclust5 <- cutree(iris.scale.hclust.ward, k=5)
table(cluster.hclust5, cluster.kmeans3)
```

---

- par une visualisation

```{r out.width="95%", fig.width=12, fig.height=6}
par(mfrow=c(1,2))
plot(iris.scaled.acp, col.ind=cluster.kmeans3, choix="ind", title="kmeans en 3 groupes", cex=0.6)
plot(iris.scaled.acp, col.ind=cluster.hclust5, choix="ind", title="hclust en 5 groupes", cex=0.6)
par(mfrow=c(1,1))
```


# Comparaison avec la réalité

## La réalité

```{r}
variete <- iris[,5]
table(variete)
```

#

```{r out.width="65%", fig.width=6, fig.height=6}
plot(iris.scaled.acp, col.ind=variete, choix="ind", cex=0.8)
```

---

# Comparer k-means avec la réalité

```{r confusion_table, results = TRUE}
conf.kmeans <- table(variete, cluster.kmeans3)
kable(conf.kmeans, caption = "Confusion table: 3-clusters k-means versus actual class")
```

---

## Setosa vs others

### Visualisation

```{r  results=TRUE}
variete2 <- rep("notSetosa", 150)
variete2[variete=="setosa"] <- "setosa"
variete2 = factor(variete2)
table(variete2)
```

---

```{r out.width="95%", fig.width=12, fig.height=6}
par(mfrow=c(1,2))
plot(iris.scaled.acp, col.ind=variete2, title="Actual species", cex=0.6)
cluster.kmeans2 <- kmeans(mes.iris.scaled, center=2)$cluster
plot(iris.scaled.acp, col.ind=cluster.kmeans2, title="2-group k-means", cex=0.6)
par(mfrow=c(1,1))
```


---

### Table de confusion et calcul de performances

```{r}
conf.kmeans <- table(variete2, cluster.kmeans2)
kable(conf.kmeans)
```

---

  - table de confusion, taux de bien prédits, spécificité, sensibilité, ...
```{r}
TP <- conf.kmeans[1,1]
FP <- conf.kmeans[1,2]
FN <- conf.kmeans[2,1]
TN <- conf.kmeans[2,2]
P <- TP + FN          # nb positif dans la réalité
N <- TN + FP          # nb négatif dans la réalité
FPrate <- FP / N      # = false alarm rate
Spe <- TN / N      # = spécificité 
Sens <- recall <- TPrate <- TP / P      # = hit rate ou recall ou sensibilité ou coverage
PPV <- precision <- TP / (TP + FP)
accuracy <- (TP + TN) / (P + N)
F.measure <- 2 / (1/precision + 1/recall)
performance <- c(FPrate, TPrate, precision, recall, accuracy, F.measure, Spe, PPV)
names(performance) <- c("FPrate", "TPrate", "precision", "recall", "accuracy", "F.measure", "Spe", "PPV")
```

---

```{r}
kable(performance, digits=3)
```

---

  - rand index et adjusted rand index
```{r}
clues::adjustedRand(as.numeric(variete2), cluster.kmeans2)
```

---

## Versicolor vs !Versicolor

### Visualisation

```{r results = TRUE}
variete2 <- rep("notVersicolor", 150)
variete2[variete=="versicolor"] <- "versicolor"
variete2 = factor(variete2)
table(variete2)
```

---


```{r out.width="90%", fig.width=12, fig.height=6}
par(mfrow=c(1,2))
plot(iris.scaled.acp, col.ind = variete2, cex  =0.7)
cluster.kmeans2 <- kmeans(mes.iris.scaled, center=2)$cluster
plot(iris.scaled.acp, col.ind = cluster.kmeans2, cex = 0.7)
par(mfrow=c(1,1))
```

---

## Table de confusion et calcul de performances

```{r}
conf.kmeans <- table(variete2, cluster.kmeans2)
kable(conf.kmeans)
TP <- conf.kmeans[1,1]
FP <- conf.kmeans[1,2]
FN <- conf.kmeans[2,1]
TN <- conf.kmeans[2,2]
P <- TP + FN          # nb positif dans la réalité
N <- TN + FP          # nb négatif dans la réalité
FPrate <- FP / N      # = false alarm rate
Spe <- TN / N      # = spécificité 
Sens <- recall <- TPrate <- TP / P      # = hit rate ou recall ou sensibilité ou coverage
PPV <- precision <- TP / (TP + FP)
accuracy <- (TP + TN) / (P + N)
F.measure <- 2 / (1/precision + 1/recall)
performance <- c(FPrate, TPrate, precision, recall, accuracy, F.measure, Spe, PPV)
names(performance) <- c("FPrate", "TPrate", "precision", "recall", "accuracy", "F.measure", "Spe", "PPV")
```

---

```{r}
kable(performance, digits=3)
clues::adjustedRand(as.numeric(variete2), cluster.kmeans2)
```

# Pros et cons des différents algorithmes

| Algorithme | Pros | Cons |
|-------------|------------------------------|------------------------|
| **Hiérarchique** | L'arbre reflète la nature imbriquée de tous les sous-clusters | Complexité quadratique (mémoire et temps de calcul) $\rightarrow$ quadruple chaque fois qu'on double le nombre d'individus  |
| | Permet une visualisation couplée dendrogramme (groupes) + heatmap (profils individuels) | |
| | Choix a posteriori du nombre de clusters (élagage) | |
|  |  |  |
|  |  |  |
| **K-means** | Rapide (linéaire en temps), peut traiter des jeux de données énormes (centaines de milliers de pics ChIP-seq) | Positions initiales des centres est aléatoire $\rightarrow$ résultats changent d'une exécution à l'autre |
| | | Distance euclidienne (pas appropriée pour transcriptome par exemple) |


#

Contact: <anne.badel@univ-paris-diderot.fr>