Importation et inspection des données
#chargement rapide
library(data.table)
X <- data.table::fread("wave100k4.txt",header=TRUE)
#premières lignes
head(X)
#structure
str(X)
## Classes 'data.table' and 'data.frame': 100000 obs. of 4 variables:
## $ X07: num -0.14 -0.06 2.4 5.02 1.73 -0.1 2.97 6.79 4.86 0.76 ...
## $ X11: num 0.53 3.99 3.85 1.44 0.75 1.52 2.43 2.44 0.77 2.94 ...
## $ X15: num 7.2 2.98 4.59 2.93 5.22 4.97 2.37 0.41 0.44 3.82 ...
## $ X10: num 1 2.07 1.7 2.9 -0.2 0.7 4.36 1.62 2.92 1.08 ...
## - attr(*, ".internal.selfref")=<externalptr>
#stat. descriptives
summary(X)
## X07 X11 X15 X10
## Min. :-3.600 Min. :-1.980 Min. :-4.330 Min. :-2.330
## 1st Qu.: 1.090 1st Qu.: 2.070 1st Qu.: 1.090 1st Qu.: 1.890
## Median : 2.510 Median : 3.170 Median : 2.510 Median : 2.990
## Mean : 2.671 Mean : 3.328 Mean : 2.665 Mean : 2.995
## 3rd Qu.: 4.230 3rd Qu.: 4.530 3rd Qu.: 4.210 3rd Qu.: 4.100
## Max. : 9.890 Max. : 9.770 Max. : 9.310 Max. : 8.510
Carte de Kohonen
#package
library(kohonen)
#carte
set.seed(0)
carte <- kohonen::som(X=as.matrix(X),grid=somgrid(7,7,"rectangular"))
print(carte$grid)
## $pts
## x y
## [1,] 1 1
## [2,] 2 1
## [3,] 3 1
## [4,] 4 1
## [5,] 5 1
## [6,] 6 1
## [7,] 7 1
## [8,] 1 2
## [9,] 2 2
## [10,] 3 2
## [11,] 4 2
## [12,] 5 2
## [13,] 6 2
## [14,] 7 2
## [15,] 1 3
## [16,] 2 3
## [17,] 3 3
## [18,] 4 3
## [19,] 5 3
## [20,] 6 3
## [21,] 7 3
## [22,] 1 4
## [23,] 2 4
## [24,] 3 4
## [25,] 4 4
## [26,] 5 4
## [27,] 6 4
## [28,] 7 4
## [29,] 1 5
## [30,] 2 5
## [31,] 3 5
## [32,] 4 5
## [33,] 5 5
## [34,] 6 5
## [35,] 7 5
## [36,] 1 6
## [37,] 2 6
## [38,] 3 6
## [39,] 4 6
## [40,] 5 6
## [41,] 6 6
## [42,] 7 6
## [43,] 1 7
## [44,] 2 7
## [45,] 3 7
## [46,] 4 7
## [47,] 5 7
## [48,] 6 7
## [49,] 7 7
##
## $xdim
## [1] 7
##
## $ydim
## [1] 7
##
## $topo
## [1] "rectangular"
##
## $neighbourhood.fct
## [1] bubble
## Levels: bubble gaussian
##
## $toroidal
## [1] FALSE
##
## attr(,"class")
## [1] "somgrid"
#noeud d'appartenance des individus - 20 premiers
print(head(carte$unit.classif,20))
## [1] 7 3 26 39 21 20 33 48 46 19 15 32 26 12 26 40 13 35 46 26
#effectifs par noeud
print(nk <- table(carte$unit.classif))
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 2213 2099 1766 2474 1764 2118 1730 2036 2249 1991 2223 2054 1815 2311 2031 1894
## 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
## 1883 2080 2788 1636 1603 2008 2474 2036 2313 2177 1995 2019 1717 2089 2254 2330
## 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
## 1870 1592 1753 2088 2175 2120 2018 1977 1947 1856 1793 2096 2199 2191 1725 2005
## 49
## 2425
#le noeud le plus "peuplé"
print(max(nk))
## [1] 2788
print(which.max(nk))
## 19
## 19
#fonction pour code couleur - degradé de bleu
degrade.bleu <- function(n){
return(rgb(0,0.5,1,alpha=seq(0,1,1/n)))
}
#effectifs dans les noeuds
plot(carte,type="count",palette.name = degrade.bleu)
#codebooks
#moyennes conditionnellement aux noeuds
print(carte$codes)
## [[1]]
## X07 X11 X15 X10
## V1 1.7970928 5.0957129 2.13818646 5.7217036
## V2 1.3743381 6.0086096 3.49156137 4.9570448
## V3 -0.5784362 3.9278500 4.16040314 2.7050619
## V4 1.1087058 4.4569350 4.60889777 3.7459237
## V5 -0.7659054 2.7488395 6.02395963 2.0813969
## V6 0.8198278 1.4690022 6.16778609 1.9984366
## V7 0.2925451 2.1421029 6.80180134 0.4264744
## V8 1.1913965 6.4651841 1.47362969 4.5947988
## V9 2.5012925 5.4757050 3.26444954 3.2792503
## V10 2.6256589 4.0374639 3.24225302 4.5613086
## V11 0.6336316 2.5190351 4.81004172 3.3081917
## V12 2.0531233 2.3975348 5.44040739 0.9179408
## V13 -0.1007406 3.1804598 5.09848764 0.2293138
## V14 0.8880825 3.5574868 5.99640353 1.8782007
## V15 3.2763116 5.9581072 1.10002910 3.6061481
## V16 0.4180288 4.3040731 2.86194640 4.4127729
## V17 0.7025150 5.4730675 2.79685067 2.8560540
## V18 1.0131049 4.3518780 4.14605475 1.6331038
## V19 0.2701448 2.2676279 4.14565088 1.4923676
## V20 -0.6948966 1.1701382 5.24057317 1.0521238
## V21 1.2163149 1.0917278 4.71782106 0.3555973
## V22 2.9631642 6.6388704 2.44197311 5.0447366
## V23 1.7229178 4.6675830 1.33393198 4.0265451
## V24 3.8017393 4.9659939 1.84237186 5.3585663
## V25 1.2917868 3.0980478 2.85952690 2.7727695
## V26 2.5281657 3.0666244 4.51718785 2.4773081
## V27 1.8549596 1.0131558 4.34814222 2.3270976
## V28 2.1422922 2.2983429 3.32103681 0.7953525
## V29 2.6904892 5.6221442 0.42926306 5.3563950
## V30 3.1491304 3.5827995 0.94421224 4.8645394
## V31 4.0001153 4.3714102 2.01253242 3.5497861
## V32 3.6104214 2.3294897 1.02695838 3.0808161
## V33 3.4220861 1.8108754 2.91042917 3.3181222
## V34 2.2584597 1.0823619 2.20627452 1.7854108
## V35 3.9102526 1.3158928 3.61053101 1.4853538
## V36 3.3523993 4.0422063 0.08111549 3.2917325
## V37 5.0313174 4.6395052 0.46584314 4.2889474
## V38 5.0419996 2.9617318 1.63331476 4.3071782
## V39 5.3975231 2.4217549 2.30519763 2.4595071
## V40 4.4833898 1.8385922 -0.44491785 2.5350817
## V41 4.1498049 1.7612525 1.39116525 1.3025677
## V42 3.3299732 3.1842295 2.45101161 1.8351995
## V43 5.9918170 3.2363289 -0.81413637 3.1446332
## V44 6.5108637 2.6750116 0.84113668 3.5831011
## V45 4.7774005 2.4692039 -0.01674306 4.5583925
## V46 5.0743216 0.9089923 1.27165874 3.3028399
## V47 6.2445521 1.3972535 -0.44530434 3.4852358
## V48 6.2523541 1.8481793 0.48161652 1.8616431
## V49 5.0755377 3.5662321 0.59247717 2.2846420
#attention, $codes est une liste
#pour obtenir le codebook du 1er noeud
print(carte$codes[[1]][1,])
## X07 X11 X15 X10
## 1.797093 5.095713 2.138186 5.721704
#affichage des codebooks - interprétation des noeuds
#possible parce que 4 variables seulement
plot(carte,type="codes")
## Warning in par(opar): argument 1 does not name a graphical parameter
#caractérisation des noeuds (bis)
#couleurs plus ou moins intenses en fonction de la valeur
#des variables
par(mfrow=c(2,2))
for (j in 1:ncol(X)){
plot(carte,type="property",property=carte$codes[[1]][,j],palette.name=degrade.bleu,main=colnames(X)[j])
}
par(mfrow=c(1,1))
CAH à partir de la carte de Kohonen
#distance entre les noeuds
D <- dist(carte$codes[[1]])
#affichage de contrôle
print(round(as.matrix(D)[1:10,1:10],2))
## V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
## V1 0.00 1.85 4.49 3.30 6.36 6.65 7.80 1.99 2.81 2.09
## V2 1.85 0.00 3.70 2.28 5.47 6.07 6.90 2.11 2.10 2.38
## V3 4.49 3.70 0.00 2.10 2.30 3.54 4.01 4.51 3.61 3.82
## V4 3.30 2.28 2.10 0.00 3.35 3.81 4.67 3.82 2.24 2.24
## V5 6.36 5.47 2.30 3.35 0.00 2.04 2.20 6.68 5.21 5.20
## V6 6.65 6.07 3.54 3.81 2.04 0.00 1.90 7.34 5.38 5.00
## V7 7.80 6.90 4.01 4.67 2.20 1.90 0.00 8.08 6.05 6.23
## V8 1.99 2.11 4.51 3.82 6.68 7.34 8.08 0.00 2.76 3.33
## V9 2.81 2.10 3.61 2.24 5.21 5.38 6.05 2.76 0.00 1.93
## V10 2.09 2.38 3.82 2.24 5.20 5.00 6.23 3.33 1.93 0.00
#cah
cah <- hclust(D,method="ward.D2",members=nk)
plot(cah,hang=-1)
#découpage en 3 groupes
groupes_cah <- cutree(cah,k=3)
print(groupes_cah)
## V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 V17 V18 V19 V20
## 1 1 2 1 2 2 2 1 1 1 2 2 2 2 1 1 1 2 2 2
## V21 V22 V23 V24 V25 V26 V27 V28 V29 V30 V31 V32 V33 V34 V35 V36 V37 V38 V39 V40
## 2 1 1 1 2 2 2 3 1 1 1 3 3 3 3 1 1 3 3 3
## V41 V42 V43 V44 V45 V46 V47 V48 V49
## 3 3 3 3 3 3 3 3 3
#effectifs par groupe
print(table(groupes_cah))
## groupes_cah
## 1 2 3
## 17 15 17
#code couleur pour groupes
couleur.groupes <- function(k){
return(c('orchid','chartreuse','gold'))
}
#carte avec groupes finaux
plot(carte,type="property",property=groupes_cah,palette.name=couleur.groupes)
Associer les individus à leur groupe définitif
#associer les individus au groupe final
finalGroup <- groupes_cah[carte$unit.classif]
print(head(finalGroup,20))
## V7 V3 V26 V39 V21 V20 V33 V48 V46 V19 V15 V32 V26 V12 V26 V40 V13 V35 V46 V26
## 2 2 2 3 2 2 3 3 3 2 1 3 2 2 2 3 2 3 3 2
#effectifs par groupe final
print(table(finalGroup))
## finalGroup
## 1 2 3
## 35711 30373 33916
#vérification pour le groupe 1
#à partir des effectifs des noeuds associés
print(sum(nk[groupes_cah == 1]))
## [1] 35711
Interprétation des groupes
#moyenne marginale de X15
print(mean(X$X15))
## [1] 2.665058
#moyenne conditionnelle des individus du groupe 2 (chartreuse)
print(mean(X$X15[finalGroup==2]))
## [1] 4.916023
#moyenne marginale de X07
print(mean(X$X07))
## [1] 2.671258
#moyenne conditionnelle des individus du groupe 2 (chartreuse)
print(mean(X$X07[finalGroup==2]))
## [1] 0.7681306