######################################################################
# 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