Practicum 8
Analysons le jeu de données iris, composé de trois espèces et tentons de regrouper les données. Ci-dessus figurent trois exemples de modèles, avec une technique différente. Dans le modèle M\(_1\), on note un groupe composé de deux individus. Idem pour M\(_4\), qui présente un groupe composé d’un seul iris. Des méthodes hiérarchiques robustes, basées sur la médiane ou une distance pondérées, offrent un meilleur résultat ici. Un graphe en coordonnées parallèles est utile ici pour voir ce que la valeur 42 du dernier dendrogramme représente de particulier. Dans ce graphique, les variables sont standardisées et représentées côte à côte en abscisse. L’ordre de la représentation graphique a une importance. La fonction cutree permet de séparer pour une distance \(h\) donnée les observations en groupe.
library(MASS); library(cluster)
data(iris)
#Tracé parallèle, avec couleurs pour chaque espèce
colnames(iris) <- c("Longueur sépales","Largeur sépales", "Longueur pétales",
"Largeur pétales", "Espèces")
groups <- c(rep(2,50),rep(3,50),rep(4,50))
style <- rep(1,150); style[42] <- 3
groups[42] <- 1 #Isoler le point 42
parcoord(iris[,c(2,1,3,4)],col=groups, lwd=style, lty=style)
Pour les modèles hiérarchiques, on utilise la fonction hclust. Cetter dernière exige comme paramètre d’entrée une matrice de distance que l’on obtient avec dist). Pour représenter un dendrogramme, on utilise la fonction plot ou plclust sur l’objet obtenu par hclust.
La première question qui se pose lors de l’utilisation d’une méthode hiérarchique est la standardisation des données. Recentrer les données, les mettre sur une échelle commune (en divisant par l’erreur standard ou l’étendue) va avoir un impact. Comme lors de l’ACP, les variables dont les mesures sont les plus grandes ont plus d’impact
#Calculer un objet distance et fournir à hclust
m1 <- hclust(d=dist(iris[,1:4]), method="single")
m2 <- hclust(d=dist(iris[,1:4],method="canberra"), method="single")
#Données normalisées, probablement une mauvaise idée ici!
m3 <- hclust(d=dist(apply(iris[,1:4],2,scale)), method="median")
m4 <- hclust(d=dist(apply(iris[,1:4],2,scale)), method="ward.D")
#Tracer l'arborescence et inclure la répartition en k groupes
plot(m1, main="Dendrogramme",sub="Iris",xlab="distance unique, Euclidienne")
rect.hclust(m1, k=3)
plot(m2, main="Dendrogramme",sub="Iris",xlab="distance unique, Canberra")
rect.hclust(m2, k=3)
plot(m3, main="Dendrogramme",sub="Iris normalisés",xlab="distance médiane, Euclidienne")
rect.hclust(m3, k=3)
plot(m4, main="Dendrogramme",sub="Iris normalisés",xlab="distance au barycentre, Euclidienne")
rect.hclust(m4, k=3)
Choix du nombre de groupes par le biais des silhouettes: on veut maximiser la valeur \(s_{i, K}\) moyenne, tout en obtenant des profils équilibrés. La fonction silhouette du paquet cluster permet la représentation graphique. La fonction agnes est l’équivalent de hclust.
par(mfrow=c(1,2))
#Même résultat avec hclust
m1 <- agnes(x=iris[,1:4],metric="canberra", method="single")
s1 <- silhouette(cutree(m1, k = 3), dist(iris[,1:4]))
plot(s1,main="Silhouette de m1",col=2:4); abline(v=mean(s1[,3]))
m5 <- hclust(d=dist(iris[,1:4]), method="complete")
s5 <- silhouette(cutree(m5, k = 3), dist(iris[,1:4]))
plot(s5,main="Silhouette de m3",col=2:4); abline(v=mean(s5[,3]))
par(mfrow=c(1,1))
plot(princomp(iris[,1:4])$scores[,1:2], col=s5[,'cluster']+1, pch=1,cex=1.5,
xlab="Composante principale 1",ylab="Composante principale 2",
main="Partitionnement des données")
points(princomp(iris[,1:4])$scores[,1:2], col=as.numeric(as.factor(iris[,5]))+1,, pch=20)
points(princomp(iris[,1:4])$scores[,1:2], col=s1[,'cluster']+1, pch=20)
Toutes les méthodes hiérarchiques donnent de mauvais résultats ici, avec des groupes qui ne sont pas équilibrés ou trop petits. Nous verrons la classification plus tard avec l’analyse discriminante pour ces même données.
On peut aussi utiliser une méthode \(K\)-moyenne pour séparer les individus. Quelquefois, l’algorithme va diviser l’espèce setosa en deux sous-espèces. On spécifie à la fonction kmeans le nombre de groupes désirées en fournissant soit les coordonnées du centre, soit en donnant un entier spécifiant le nombre de groupes souhaités.
kmoy <- kmeans(iris[,1:4],centers=3)
plot(1:150, kmoy$cluster, col=groups, pch=4)
dE2 <- daisy(iris[,1:4])^2 #Calcule distance
sk2 <- silhouette(kmoy$cl, dE2)Comme vu en classe, pam, pour partition autour de médoïdes, est une version robuste de l’algorithme de \(k\)-médianes pour \(k\) arbitraire. Plus robuste que les \(k\)-moyennes, cet algorithme d’apprentissage non supervisé sélectionne à partir de réprésentants (les médoïdes) et d’une matrice de dissimilarité (dont on minimise la somme, au contraire de la somme des distances euclidiennes entre observations). Puisque pam est un algorithme gourmand (\(k\)-médoïdes est en général NP complet), l’algorithme clara est recommandé pour les jeux de données de plus de 2000 observations.
par(mfrow=c(1,2))
kpam <- pam(iris[,1:4], k=3)
plot(sk2,col=c(2:4),main="Silhouette, K-moyennes")
abline(v=mean(sk2[,3]))
plot(s.pam <- silhouette(kpam),main="Silhouette avec PAM",col=c(2:4))
abline(v=mean(s.pam[,3]))
Le résultat est plus probant, les setosa étant clairement séparés des autres (normalement, le partitionnement des données sert pour les cas où on ne connaît pas les étiquettes des groupes; le cas actuel est donc plutôt atypique). On profite de l’occasion pour illustrer les diagrammes d’Andrews sur les iris.
par(mfrow=c(1,1))
library(andrews)
andrews(iris,type=4,clr=5,ymax=2)
Encore une fois, les versicolor et virginica sont très semblables et les Iris setosa ressortent du lot. Les faces de Chernoff permettent de distinguer visuellement les aberrances et les caractéristiques.
aplpack::faces(iris[sort(sample.int(150,20)),1:4])
## effect of variables:
## modified item Var
## "height of face " "Longueur sépales"
## "width of face " "Largeur sépales"
## "structure of face" "Longueur pétales"
## "height of mouth " "Largeur pétales"
## "width of mouth " "Longueur sépales"
## "smiling " "Largeur sépales"
## "height of eyes " "Longueur pétales"
## "width of eyes " "Largeur pétales"
## "height of hair " "Longueur sépales"
## "width of hair " "Largeur sépales"
## "style of hair " "Longueur pétales"
## "height of nose " "Largeur pétales"
## "width of nose " "Longueur sépales"
## "width of ear " "Largeur sépales"
## "height of ear " "Longueur pétales"Méthode basée sur un modèle: mélange de lois normales
Voir https://cran.r-project.org/web/packages/mclust/vignettes/mclust.html pour des exemples d’utilisation du paquet mclust. Chaque classe est modélisée comme une variable multinormale et le graphique de base représente le critère de Schwartz (BIC) pour les différents modèles de matrices de covariance en fonction du nombre de groupes. La fonction choisit le modèle qui donne la plus grande valeur de \(-\)BIC.
suppressPackageStartupMessages(library(mclust))
irisBIC <- mclustBIC(iris[,-5])
plot(irisBIC, xlab="Nombre de composantes")
irisSummary3 <- summary(irisBIC, data = iris[,-5], G = 3)
irisclust2 <- Mclust(iris[,-5], G=2,modelNames = "VEV")
irisclust3 <- Mclust(iris[,-5], G=3,modelNames = "VEV")
coordProj(data = iris[,-5], dimens = c(2,4), what="classification",
parameters = irisSummary3$parameters, z= irisSummary3$z)
title("Classification")
coordProj(data = iris[,-5], dimens = c(2,4), what="uncertainty",
parameters = irisSummary3$parameters, z= irisSummary3$z)
title("Incertitude")
coordProj(data = iris[,-5], dimens = c(2,4), what="errors",
parameters = irisSummary3$parameters, z= irisSummary3$z,truth = iris[,5])
title("Erreurs")
par(mfrow=c(1,2))
plot(sclust2 <- silhouette(irisclust2$classification, dist(iris[,-5])),main="Silhouette pour mclust",col=c(2:3))
abline(v=mean(sclust2[,3]))
plot(sclust3 <- silhouette(irisclust3$classification, dist(iris[,-5])),main="Silhouette pour mclust",col=c(2:4))
abline(v=mean(sclust3[,3]))
Si nous nous sommes concentrés en classe sur les méthodes de partionnement des données basées sur les modèles gaussiens, il existe de nombreuses extensions plus récentes qui font appel à d’autre distributions, principalement des lois elliptiques (de Student-\(t\) (teigen), normale-inverse-gaussienne ou hyperboliques ("MixGHD) et leurs extensions asymmétriques). Ces dernières peuvent capturer les queues lourdes et l’asymmétrie souvent présente dans les groupes. D’autres paquets s’attardent à la fusion de classes gaussiennes (lorsque la représentation elliptique n’est pas appropriée, il est possible que mclust retourne un nombre de groupes plus élevé que nécessaire.
Autres techniques et jeux de données pour pratique
D’autres technique que nous n’avons pas couvert en classe (tirées de la vaste littérature en apprentissage automatique) incluent la partionnement spectral. Cette dernière est très utile pour les groupes connexes, mais dont les frontières ne sont pas convexes. La méthode permet une transformation des données, qui dans un nouvel espace peuvent être plus facilement partitionnées.
On note aussi le partionnement par délocalisation (mean shift). Ces méthodes fonctionnent bien quand les groupes ne sont pas nécessairement ellipsoidaux.
Il est également possible de démarrer les \(K\)-moyennes à des centres initiaux à l’aide de l’algorithme \(K\)-moyennes ++, qui utilise la matrice de dissimilarité pour obtenir de meilleurs estimés.
library(kernlab)
sclust <- specc(as.matrix(iris[,-5]), centers=3) #aussi via clusterSim::speccl
plot(iris[,-5], col=sclust, pch=4) 
library(LPCM)
msclust <- ms(iris[,-5], h=0.12) #h bande passante
detach("package:kernlab", unload=TRUE)
library(flexclust, quietly=TRUE)
kmeanppclust <- flexclust::kcca(x=as.matrix(iris[,-5]),k=3) #K-moyennes++
kmedppclust <- flexclust::kcca(x=as.matrix(iris[,-5]),k=3, family=kccaFamily("kmedians"))
plot(iris[,-5], col=clusters(kmeanppclust), pch=4) 
plot(iris[,-5], col=clusters(kmedppclust), pch=4) 
Cette fonction fait aussi le boulot pour \(K\)-moyennes++.
kmpp <- function(X, k) {
n <- nrow(X)
C <- numeric(k)
C[1] <- sample(1:n, 1)
for (i in 2:k) {
dm <- pracma::distmat(X, X[C, ])
pr <- apply(dm, 1, min); pr[C] <- 0
C[i] <- sample(1:n, 1, prob = pr)
}
kmeans(X, X[C, ])
}
kmpp(as.matrix(iris[,-5]),3)Vous pouvez apréhender les diverses méthodes sur quelques exemples tirées du paquet cluster.datasets ou générer vos propres exemples par le biais de clusterSim et mlbench.
library(clusterSim)
dat1 <- shapes.two.moon(numObjects=200)
plot(dat1$data, col=dat1$clusters, xlab="x",ylab="y")
dat2 <- shapes.worms(180)
plot(dat2$data,col=dat2$clusters, xlab="x",ylab="y")
dat3 <- shapes.blocks3d(300,1,3,3,3)
plot(dat3$data,col=dat3$clusters, xlab="x",ylab="y")
dat4 <- shapes.circles3(numObjects=c(120,180,240))
plot(dat4$data,col=dat4$clusters, xlab="x",ylab="y")
library(mlbench)
dat5 <- mlbench.smiley(n=500, sd1 = 0.2, sd2 = 0.2)
plot(dat5$x,col=dat5$classes, xlab="x",ylab="y")
dat6 <- mlbench.spirals(100,1,0.025)
plot(dat6$x,col=dat6$classes, xlab="x",ylab="y")
#dat7 <- shapes.bulls.eye(numObjects=c(100,100))
#plot(dat7$data,col=dat7$clusters, xlab="x",ylab="y")