M1 Classification 2019/2020 - Correction du partiel

Initialisation du générateur de nombres aléatoires (permet d’obtenir les mêmes tirages au sort à chaque exécution du code) :

set.seed(10)

Exercice 1

1.

La dissimilarité de Dice est une mesure de dissimilarité entre des observations de variables binaires. En notant \(X_i^k\in\{0,1\}\), \(k\in\{1,\ldots,P\}\), les variables observées sur un individu \(e_i\) et \(X_j^k\in\{0,1\}\) les variables observées sur un autre individu \(e_j\), la dissimilarité de Dice est définie par \[\delta(e_i,e_j)=1-\frac{2a_{ij}}{2a_{ij}+b_{ij}+c_{ij}},\] où \(a_{ij}\) est le nombre d’indices \(k\in\{1,\ldots,P\}\) pour lesquels \(X_i^k=X_j^k=1\), \(b_{ij}\) le nombre d’indices \(k\) pour lesquels \(X_i^k=1\) et \(X_j^k=0\), et \(c_{ij}\) le nombre d’indices pour lesquels \(X_i^k=0\) et \(X_j^k=1\).

2.

• Pour \(e_1\) et \(e_2\) : on a \(a_{12}=2\), \(b_{12}=0\) et \(c_{12}=1\), donc \[\delta(e_1,e_2)=1-\frac{2\times2}{2\times2+0+1}=\frac15=0.2\]
• Pour \(e_1\) et \(e_3\) : on a \(a_{13}=0\), \(b_{13}=2\) et \(c_{13}=1\), donc \[\delta(e_1,e_3)=1-\frac{0}{0+2+1}=1\]
• Pour \(e_2\) et \(e_3\) : on a \(a_{23}=1\), \(b_{23}=2\) et \(c_{23}=0\), donc \[\delta(e_2,e_3)=1-\frac{2}{2+2+0}=\frac12=0.5\]

3.

tab = data.frame(X1=c(1,1,0),X2=c(0,1,1),X3=c(1,1,0),row.names=c("e1","e2","e3"))
tab

library(ade4)
delta = dist.binary(tab,method=5)
delta

          e1        e2
e2 0.4472136          
e3 1.0000000 0.7071068

On ne retrouve pas les mêmes valeurs (sauf pour \(\delta(e_1,e_3)\)) car dans la définition utilisée par dist.quant, il y a une racine carrée : autrement dit, les dissimilarités obtenues ici sont les carrées des dissimilarités telles que définies en cours. Il suffit donc de mettre au carré les valeurs pour retrouver les calculs de la question 2 :

delta^2

    e1  e2
e2 0.2    
e3 1.0 0.5

Exercice 2

library(ade4)
data("elec88")
tab = elec88$tab
tab

1.

Il s’agit de variables quantitatives correspondant à des pourcentages, par conséquent toutes les variables sont exprimées dans la même “unité” et donc directement comparables. On va donc choisir de calculer des distances euclidiennes simples.

remarque : on pourrait penser qu’il est plus logique d’utiliser une distance euclidienne normalisée car les ordres de grandeur des pourcentages entre certains candidats sont très différents et donc les variables ont des variances empiriques très différentes. Cependant utiliser une distance euclidienne normalisée donnerait une importance trop grande aux “petits candidats” dans l’évaluation des différences entre les observations, ce qui n’est a priori pas souhaitable pour analyser de telles données.

# calcul des distances euclidiennes
delta = dist.quant(tab,method = 1)

2.

# CAH avec distances euclidiennes et stratégie de Ward
cah = hclust(delta,method="ward.D2")
plot(cah)

3.

Le saut maximal du critère est obtenu pour le passage de 4 à 3 classes, ce qui veut dire que le nombre de classes à choisir est de 4.

Détermination de la partition pour 4 classes

cl_cah = cutree(cah,k=4)
cl_cah

 D1  D2  D3  D4  D5  D6  D7  D8  D9 D10 D11 D12 D13 D14 D15 D16 D17 D18 D19 D21 D22 D23 
  1   2   2   1   1   3   1   1   2   1   2   1   3   1   4   2   1   2   4   1   2   2 
D24 D25 D26 D27 D28 D29 D30 D31 D32 D33 D34 D35 D36 D37 D38 D39 D40 D41 D42 D43 D44 D45 
  2   1   1   1   1   1   3   2   2   2   3   1   2   1   1   1   2   1   1   1   1   1 
D46 D47 D48 D49 D50 D51 D52 D53 D54 D55 D56 D57 D58 D59 D60 D61 D62 D63 D64 D65 D66 D67 
  2   1   1   1   1   1   1   1   1   1   1   1   2   1   1   1   2   1   1   2   3   3 
D68 D69 D70 D71 D72 D73 D74 D75 D76 D77 D78 D79 D80 D81 D82 D83 D84 D85 D86 D87 D88 D89 
  3   1   1   1   2   1   1   4   2   1   1   1   1   1   1   3   3   1   2   2   1   1 
D90 D91 D92 D93 D94 D95 
  1   1   1   3   2   1 

4.

Partition obtenue par k-moyennes pour k=4

km = kmeans(tab,4)
cl_km = km$cluster

Comparaison des deux classifications

table(cl_cah,cl_km)

      cl_km
cl_cah  1  2  3  4
     1  0  9 39 10
     2  1  1  0 21
     3  0 10  0  0
     4  3  0  0  0

Les deux partitions sont assez différentes, même si on peut noter certaines similarités faibles : - la classe 1 des k-moyennes correspond approximativement à la classe 3 de la CAH, - la classe 3 de la CAH est incluse dans la classe 2 des k-moyennes.

5.

On pourrait chercher à affecter un poids \(p_i\) à chaque département (donc chaque observation) dépendant du nombre d’habitants \(N_i\). On pourrait prendre simplement par exemple \(p_i=N_i\). Voici comment alors ces poids pourrait intevenir dans les méthodes : - pour la CAH, ces poids pourraient modifier la formule de calcul dand la stratégie d’aggrégation de Ward. En effet, cette stratégie consiste à calculer des inerties intra et inter-classes. L’inertie d’un nuage de points \(X_i\in\R^P\), \(i\in\{1,\ldots,n\}\) est définie par la formule \[I(X) = \frac1n\sum_{i=1}^n\|X_i-\bar X\|^2,\quad\text{ avec }\bar X=\frac1n\sum_{i=1}nX_i.\] Cette formule pourrait être modifiée ainsi pour tenir compte des pondérations : \[I(X) = \frac1n\sum_{i=1}^np_i\|X_i-\bar X\|^2,\quad\text{ avec }\bar X=\frac1n\sum_{i=1}^np_iX_i.\] - pour les k-moyennes, on pourrait de même faire intervenir les poids dans le calcul des barycentres de chaque groupe, suivant la même formule que celle définie ci-dessus pour \(\bar X\).

6.

?kmeans
km$totss

[1] 6588.886

km$withinss

[1]  244.1800  788.2300 1046.1272  943.1948

km$tot.withinss

[1] 3021.732

km$betweenss

[1] 3567.154

• totss siginifie “total sum of squares” ; il s’agit en fait de l’intertie totale, au coefficient de normalisation près.
• withinss signifie “within-sum of squares” ; c’est le vecteur des inerties (au coefficient de normalisation près) de chaque groupe dans la partition
• tot.withinss signifie “total within-sum of squares” ; c’est la somme du vecteur withinss, et ça correspond donc à l’inertie intra (au coefficient de normalisation près)
• enfin betweenss signifie “between sum of squares” et correspond donc à l’inertie inter (toujours à coefficient de normalisation près)

Ces quantités sont données dans la sortie de kmeans car c’est précisément la quantité tot.withinss que l’algorithme des k-moyennes cherche à minimiser par une stratégie de minimisation alternée, comme on l’avait vu dans le cours. Or la stratégie de Ward cherche également à minimiser l’inertie intra (ou maximiser l’inertie inter, ce qui revient au même), même si l’algorithme qui en résulte est très différent puisqu’il s’agit d’une classification ascendante hiérarchique.

On peut donc directement comparer la performance des deux méthodes en comparant les inerties intra. Il faut juste faire attention au coefficient de normalisation poru que les valeurs soient comparables.

La valeur de l’inertie intra pour la CAH est précisément la valeur du critère obtenue à l’étape aboutissant à une partition en 4 classes. Sur le dendrogramme, cette valeur semble être autour de 33 ; mais on peut récupèrer sa valeur exacte en affichant cah$height

cah$height

 [1]  1.256981  1.442221  1.449138  1.489966  1.509967  1.667333  1.673320  1.679286
 [9]  1.702939  1.732051  1.783255  1.881489  1.959592  1.971463  1.974842  1.994994
[17]  2.034699  2.144761  2.181742  2.281812  2.370654  2.391652  2.509980  2.537716
[25]  2.580698  2.600000  2.638181  2.649528  2.764055  2.884441  2.954093  3.018830
[33]  3.028751  3.039737  3.243455  3.272614  3.419357  3.457359  3.525148  3.644173
[41]  3.743439  4.176123  4.248529  4.279408  4.331282  4.338202  4.433960  4.532108
[49]  4.882622  4.947053  5.037195  5.088035  5.100000  5.113195  5.199679  5.242082
[57]  5.312250  5.704969  5.810500  5.980524  6.344551  6.373121  6.467869  6.634255
[65]  6.923487  6.974812  7.311778  7.343932  7.464136  7.729315  7.841283  7.861679
[73]  8.718945  9.096739  9.123961  9.546727  9.793942  9.933469 11.854624 13.001641
[81] 14.467331 15.553179 15.707641 15.780368 17.073733 18.118157 20.662510 23.920008
[89] 26.943182 33.868641 43.736449 45.334620 51.010111

Il s’agit de la valeur d’indice 90, égale à 33.9.

Cette valeur est donc à comparer à tot.withinss poru les k-moyennes, en renormalisant la valeur par le nombre d’observations :

km$tot.withinss / nrow(tab)

[1] 32.14609

On a \(32.1 < 33.9\) donc sur ce critère la méthode des k-moyennes aboutit à une meilleure partition, mais les deux valeurs sont en fait assez proches, et donc la différence n’est pas vraiment significative.

7.

La méthode des k-moyennes donne des résultats différents à chaque essai lorsqu’on entre simplement le nombre de classes comme paramètre (et pas les positions des centres initiaux), car alors les \(k\) centres initiaux sont choisis au hasard parmi les observations.

Voici quelques essais et les tables de contingence des partitions obtenues :

km1 = kmeans(tab,4)
cl_km1 = km1$cluster
km2 = kmeans(tab,4)
cl_km2 = km2$cluster
km3 = kmeans(tab,4)
cl_km3 = km3$cluster
table(cl_km1,cl_km2)

      cl_km2
cl_km1  1  2  3  4
     1  0  0  4  0
     2 31  0  0  0
     3  0  0  0 39
     4  0 20  0  0

table(cl_km1,cl_km3)

      cl_km3
cl_km1  1  2  3  4
     1  3  1  0  0
     2 13  0  0 18
     3  0 13  0 26
     4  0 10 10  0

table(cl_km2,cl_km3)

      cl_km3
cl_km2  1  2  3  4
     1 13  0  0 18
     2  0 10 10  0
     3  3  1  0  0
     4  0 13  0 26

Ici les essais 1 et 2 donnent des partitions identitiques, et l’essai 3 a donné une partition différente.

km1$tot.withinss

[1] 3021.732

km2$tot.withinss

[1] 3021.732

km3$tot.withinss

[1] 3699.058

Les deux premiers essais ont donc abouti à une meilleure partition.

8.

Lorsque nstart est fixé à une valeur plus grande que 1, la partition des k-moyennes est répétée nstart fois, avec des choix aléatoires des centres initiaux, et la partition donnant la plus petite valeur de tot.withinss est conservée. Ceci permet de réduire l’aléa car les valeurs de tot.withinss obtenues seront plus proches de celle de la partition idéale, et donc les partitions elles-mêmes seront probablement plus semblables.

9.

On effectue 100 essais avec nstart=1 (qui est la valeur par défaut) ; on enregistre les valeurs de tot.withinss correspondantes dans un vecteur et on calcule moyenne et variance empiriques du vecteur obtenu :

N = 100
vec = rep(0,N)
for(l in 1:N)
{
  vec[l] = kmeans(tab,4)$tot.withinss
}
mean(vec)

[1] 3101.888

sd(vec)

[1] 168.2828

A présent on effectue la même chose avec nstart=10 :

N = 100
vec = rep(0,N)
for(l in 1:N)
{
  vec[l] = kmeans(tab,4,nstart=10)$tot.withinss
}
mean(vec)

[1] 3021.732

sd(vec)

[1] 0

On remarque que l’écart-type est nul, ce qui signifie que les 100 partitions obtenues sont en fait toutes identiques, alors que ce n’était pas le cas avec nstart=1. De plus la valeur de tot.withinss est plus faible que la moyenne des tot.withinss pour nstart=1, ce qui était attendu, et donc la partition est meilleure.

10.

On fait à présent varier le nombre de classes k. Les valeurs possibles pour k vont de 1 à n=94 le nombre d’observations. Cependant la valeur \(k=n\) ne semble pas autorisée pour la fonction kmeans, donc on va aller jusqu’à \(n-1\) seulement. On peut décider d’autre part de fixer nstart à 10 pour ces essais, puisqu’on a vu que ça stabilisait les résultats.

n = nrow(tab)
vec = rep(0,n-1)
for(k in 1:(n-1))
  vec[k] = kmeans(tab,k,iter.max=10,nstart=10)$tot.withinss
plot(vec)

Le saut maximal est cette fois obtenu pour le passage de \(k=1\) à \(k=2\), ce qui amène à retenir la partition à deux classes. Ce n’est donc pas le même nombre de classes que celui retenu avec la CAH.