######################################################################
# TP 3 : Clustering avec le logiciel R 
#       (code R avec commentaires)
######################################################################


#===================================
# Exercice 1 : un petit exemple simple
#===================================
# Création du jeu de données
X <- matrix(c(5,4,4,5,1,-2,0,-3),4,2,byrow=TRUE)
colnames(X) <- c("X1","X2")
rownames(X) <- c("1", "2","3","4")
# Visualisation des 4 points
plot(X, pch=16, xlim=c(-3.5,5.5), ylim=c(-3.5,5.5))
text(X,labels=1:4,pos=4)
abline(h=0,lty=2)
abline(v=0,lty=2)

# kmeans
# On utilise la fonction kmeans en prenant les deux premiers points comme centres initiaux
init <-X[1:2,] 
km <- kmeans(X, centers=init) # tous les résultats sont stockés dans l'objet km

km$cluster    # permet d'avoir la partition en 2 classes
## 1 2 3 4 
## 2 2 1 1
              # ici les deux classes sont C1={3,4} et C2={1,2}

# On s'intéresse maintenant à la qualité de cette partition
Tot <- km$totss      #  cette sortie la somme de carrés totaux
W <- km$tot.withinss # cette sortie donne la somme de carrés intra-classe 
B <- km$betweenss # cette sortie donne la somme des carrés inter-classe
B/Tot*100 # donne le pourcentage d'inertie expliquée
## [1] 97.01493
# NB1 : ici les poids des individus sont tous w_i=1. Dans ce cas on parle 
#       de somme de carrés totaux plutôt que d'inertie
# NB2 : en anglais on parle de total, within, between sum of squares
# NB3 : on peut vérifier que l'inertie totale (T) est la somme de l'inertie intra (W)
#       et de l'inertie inter (B).
#       On en déduit que le pourcentage d'inertie expliqué par la partition est B/T ou 1-W/T
# NB4 : Ici, 97% de l'inertie des données est expliquée par la partition en deux classes.
# hclust avec lien maximum
# La fonction hclust prend en entrée une matrice de distance (ou de dissimilarité)
# qui peut être obtenue avec la fonction dist.

d <- dist(X) # permet de calculer les distances euclidiennes entre les points
tree <- hclust(d, method="complete") # tous les résultats sont stockés dans l'objet tree
tree$height # donne la hauteur des classes dans le dendrogramme
## [1] 1.414214 1.414214 8.944272
# Graphique du dendrogramme
plot(tree, main="CAH du lien maximum", xlab ="", sub="", hang=-1) 

cutree(tree,k = 2) # permet d'obtenir la partition en deux classes
## 1 2 3 4 
## 1 1 2 2
              # ici les deux classes sont C1={1,2} et C2={3,4}
# hclust avec Ward
# Pour retrouver vraiment le dendrogramme de Ward, il faut utiliser la méthode "ward.D" 
# avec le carré des distances euclidiennes divisées par 2

d <- dist(X) 
tree <- hclust(d^2/2, method="ward.D")
sum(tree$height)
## [1] 67
# NB1 : dans cet exercice, les poids des individus tous égaux à 1. On passe donc
#       en argument d^2/2 où d contient les distances euclidiennes entre les individus.
#       Une explication est fournie à la fin de l'énoncé du TP
# NB2 : avec ce paramétrage, on retrouve bien que la somme des hauteurs de l'arbre vaut 67
#       c'est à dire l'inertie totale (ici la somme des carrés totaux) calculée 
#       à la main et trouvée précédemment avec kmeans.

# Graphique du dendrogramme

plot(tree, main="CAH de Ward", xlab ="", sub="",hang=-1)

cutree(tree,k = 2) # On retrouve toujours les deux classes sont C1={1,2} et C2={3,4}
## 1 2 3 4 
## 1 1 2 2
#===================================
# Exercice 2 : Jeu de données fromages
#===================================
# Création du jeu de données
cheese <- read.table("../data/fromage.txt", header=TRUE, row.names=1)
# Pré-traitement des données
apply(cheese, 2, sd) # écart-type très différents
##    calories      sodium     calcium     lipides     retinol     folates 
##   91.914356  108.678923   72.528882    8.129642   24.163098   11.723339 
##   proteines cholesterol   magnesium 
##    6.959788   28.245755   11.318388
n <- nrow(cheese)
Z <- scale(cheese,center=TRUE,scale=TRUE)*sqrt(n/(n-1)) # on va donc centrer-réduire les données

# COMMENTAIRES : les écart-types des variables sont très différents. L'écart-type de protein
# est de 6.9 tandis que celui de sodium est de 108.6. La variable sodium aura donc
# un poids plus important dans le calcule des distances Euclidiennes entre les eaux
# et donc dans le clustering. On risque de mettre dans les mêmes classes des fromages
# qui se ressemblent sur les variables de forte variance. Pour éviter cela, on standardise
# les données avant d'appliquer le clustering. La variance de chaque variable est alors
# égale à 1.
# hclust avec Ward des données standardisées
d <- dist(Z) 
tree <- hclust(d^2/(2*n), method="ward.D") # ici les poids des individus sont 1/n
sum(tree$height) # on retrouve que l'inertie totale est égale à 9 (nombre de variables)
## [1] 9
# COMMENTAIRES : Afin de pouvoir retrouver exactement les hauteurs de la méthode de Ward
# avec la fonction hclust, il faut utiliser d^2/(2*n) (cf. détails à la fin de l'énoncé du TP).
# Avec ce paramètrage, on a bienla somme des hauteurs de l'arbre qui est égale
# à l'inertie totale des données. L'inertie totale est la somme de variances des variables.
# Ici les données sont standardisées, donc les variables ont une variance de 1
# et l'inertie des données (dans Z) vaut p=9. La somme des hauteurs vaut bien 9.
# Choix du nombre de classes
# Graphique du dendrogramme
plot(tree, hang=-1, main="CAH de Ward", sub="", xlab="") # 5 classes ?

# Graphique des hauteurs de l'arbre
barplot(sort(tree$height, decreasing = TRUE), main="Hauteurs de l'arbre")
abline(h = 0.4, col="red") # 4 ou 5 classes ?

# NB1 : la hauteur d'une classe dans le dendrogramme est la dissimilarité (mesure d'agrégation)
# entre les deux classes qui ont été agrégées pour l'obtenir. Si deux classes sont agrégées
# à une hauteur proche d'une des hauteurs des deux sous-classes, cela signifie que ces 2 classes 
# devraient être agrégées ensembles. On évite alors de couper l'arbre à ce niveau là. 
# C'est le cas pour la "grande classe" qui va de Pyrénées à Morbier par exemple. Rediviser
# cette classe en deux sous-classes ne serait pas forcément pertinent. Mais cela reste
# malgré tout une règle "au doigt mouillé". Ici les partitions en 2, 3, 4 et 5 classes
# semblent pertinentes. 
#
# NB2 : le graphique des hauteurs de l'arbre donne la même d'information (hauteur des classes)
# sous une forme différente) et aide également au choix du nombre de classes.
# La première barre donne la perte d'inertie expliquée lorsqu'on agrége les deux 
# dernières classes. On voit qu'on a intérêt à ne pas faire cette agrégation et à conserver deux classes.
# Idem  pour la seconde barre (on préfère garder 3 classes), la troisième (4 classes)
# et la quatrième barre (5 classes). A partir de 6 classes, les agrégation font
# perdre toujours à peu près la même quantité d'inertie (coude dans le barplot).
# Cela confirme l'idée d'une partition en 5 classes.
# Partition en 5 classes 
K <- 5
part <- cutree(tree,k=K)

# Interprétation via l'ACP
library(FactoMineR)
part <- as.factor(part)
levels(part) <- paste("cluster",1:K,sep="")
res <- PCA(data.frame(part, cheese), quali.sup=1,graph=FALSE) # ACP normée
res$eig # bonne représentation sur le plan factoriel 1-2 (76% d'information retrouvée)
##         eigenvalue percentage of variance cumulative percentage of variance
## comp 1 5.049113814            56.10126460                          56.10126
## comp 2 1.844418362            20.49353736                          76.59480
## comp 3 0.867800921             9.64223245                          86.23703
## comp 4 0.577384673             6.41538526                          92.65242
## comp 5 0.355210155             3.94677950                          96.59920
## comp 6 0.175311592             1.94790658                          98.54711
## comp 7 0.097434214             1.08260237                          99.62971
## comp 8 0.028431123             0.31590137                          99.94561
## comp 9 0.004895146             0.05439051                         100.00000
plot(res,habillage=1,invisible="quali",title="") # graphique des fromages sur le plan 1-2