Eléments de correction du TP d’analyse de données

Master 2 Nutrition

Présentation du problème
La composition chimique des six éléments les plus importants (bicarbonates HCO3, sulfates SO4, chlorure
Cl, calcium Ca, magnésium Mg, sodium Na) a été relevée sur 20 bouteilles d’eau minérale. L’objectif est de

a. décrire les différentes variables et leurs relations,

b. dresser une typologie des eaux.

Questions
1. Charger les données Eaux.txt puis les packages FactoMineR et factoextra avec les commandes
library(FactoMineR) et library(FactoMineR).

eaux <- read.table("Eaux.txt",sep="\t",header=T)

head(eaux)

## HCO3 SO4 CL CA MG NA. ident

## 1 341 27 3 84 23 2 Aix
## 2 1300 40 40 190 85 150 Badoit
## 3 263 23 9 91 5 3 Beckeris
## 4 287 3 5 44 24 23 Cayranne
## 5 298 9 23 96 6 11 Chambon
## 6 200 15 8 70 2 4 Cristalline

row.names(eaux) <- eaux$ident

2. Calculer quelques statistiques simples sur les variables (moyennes, écart-types, corrélations).

apply(eaux[,1:6],2,mean)

## HCO3 SO4 CL CA MG NA.

## 464.84 154.64 35.64 118.96 21.48 130.00

apply(eaux[,1:6],2,sd)

## HCO3 SO4 CL CA MG NA.

## 845.72324 365.20301 74.95325 130.59239 29.02172 401.92910

cor(eaux[,1:6])

1
## HCO3 SO4 CL CA MG
## HCO3 1.00000000 0.0530771356 0.75991803 0.06016449 0.10805261
## SO4 0.05307714 1.0000000000 -0.01112288 0.95508574 0.81424692
## CL 0.75991803 -0.0111228812 1.00000000 -0.04532338 -0.09274089
## CA 0.06016449 0.9550857402 -0.04532338 1.00000000 0.88686404
## MG 0.10805261 0.8142469206 -0.09274089 0.88686404 1.00000000
## NA. 0.77751167 -0.0004516221 0.99297333 -0.04576128 -0.07878847
## NA.
## HCO3 0.7775116724
## SO4 -0.0004516221
## CL 0.9929733306
## CA -0.0457612773
## MG -0.0787884696
## NA. 1.0000000000

3. Pour réaliser une ACP sur les données :

a. Quelles variables utiliser?

b. Faut-il normer ou non les variables?

Les variables actives en ACP sont quantitatives. Pour réaliser une ACP sur les eaux, nous utiliserons les
6 teneurs mesurées sur les différents échantillons d’eau (bicarbonates HCO3, sulfates SO4, chlorure Cl,
calcium Ca, magnésium Mg, sodium Na). Les statistiques simples (moyennes et écart-types) calculées
à la question 2 montrent que les ordres de grandeur des différentes teneurs sont très différents avec
une variabilité différente sur chaque variable. En effet, la teneur en magnésium est la variable qui a le
plus faible écart-type (environ 29) alors que la teneur en bicarbonate présente l’écart-type le plus élevé
(environ 845 soit près de 30 fois plus que pour le magnésium). Si les variables n’étaient pas normées,
l’ACP ne tiendrait quasiment pas compte des deux teneurs les moins variables (Mg et CL). En normant
les variables, un même poids (une même importance) sera accordé aux 6 teneurs.

library(FactoMineR)
library(factoextra)

## Loading required package: ggplot2

## Welcome! Related Books: `Practical Guide To Cluster Analysis in R` at https://goo.gl/13EFCZ

acp.res<-PCA(eaux[,1:6],graph=F) # donnees normalisees par defaut

acp.res$eig # valeurs propres pour le choix du nbr d'axes

## eigenvalue percentage of variance

## comp 1
2.789615328 46.49358880
## comp 2
2.686215357 44.77025595
## comp 3
0.335429722 5.59049537
## comp 4
0.149646400 2.49410667
## comp 5
0.034194889 0.56991482
## comp 6
0.004898304 0.08163839
## cumulative percentage of variance
## comp 1 46.49359
## comp 2 91.26384
## comp 3 96.85434

2
## comp 4 99.34845
## comp 5 99.91836
## comp 6 100.00000

fviz_eig(acp.res)

Scree plot

40
Percentage of explained variances

30

20

10

1 2 3 4 5 6
Dimensions
D’après les valeurs des inerties cumulées, les deux premières composantes principales portent à elles seules
91.26% de la variabilité totale des données. Conserver les deux premières composantes principales suffit donc
largement pour l’interprétation. La représentation graphique des inerties portées par les 6 composantes de
l’ACP mène au même choix. On représente donc les variables et les eaux sur le premier plan factoriel.

fviz_pca_var(acp.res,axes=c(1,2))

3
 Variables − PCA
1.0 NA.
CL
HCO3

0.5

SO4
CA
MG
Dim2 (44.8%)

0.0

−0.5

−1.0

−1.0 −0.5 0.0 0.5 1.0

Dim1 (46.5%)

fviz_pca_ind(acp.res,axes=c(1,2))

4
 Individuals − PCA
St_Yorre

5.0
Dim2 (44.8%)

Vichy_Celestins
2.5

Badoit Hepar
Contrex

Vittel
0.0 Perrier Ondine
Ogeu
Ferita
Chambon
St_Hyppo
Thonon
St_Cyr Aix
Evian
Cayranne
Beckeris
Laurier
Veri
Cristalline
Ribes
Volvic
Spa
Viladrea

−2 0 2 4 6
Dim1 (46.5%)

4. Pourquoi serait-il judicieux de faire une ACP sur le logarithme des différentes variables? Réaliser cette
ACP.

Le graphe des individus est peu lisible. Ceci est dû à la présence dans l’échantillon de quelques eaux (Hépar
et Contrex en particulier) très différentes des autres du point de vue de leur composition. Ceci a pour effet
d’“écraser” le reste du nuage. La transformation logarithmique permettra de représenter de façon plus dilatée
les autres eaux et de mettre en évidence des différences plus fines entre leurs compositions.

# Corrélations entre les teneurs log-transformées

cor(log(eaux[,1:6]))

## HCO3 SO4 CL CA MG NA.

## HCO3 1.0000000 0.3911386 0.44616482 0.5865177 0.53313976 0.4598051
## SO4 0.3911386 1.0000000 0.27831193 0.5736158 0.53933746 0.1916941
## CL 0.4461648 0.2783119 1.00000000 0.2545694 0.04859288 0.8591840
## CA 0.5865177 0.5736158 0.25456937 1.0000000 0.72956091 0.1649545
## MG 0.5331398 0.5393375 0.04859288 0.7295609 1.00000000 0.1579806
## NA. 0.4598051 0.1916941 0.85918400 0.1649545 0.15798059 1.0000000

# ACP normée
acp.res.log<-PCA(log(eaux[,1:6]),graph=F)
# Valeurs propres de l'ACP
acp.res.log$eig

## eigenvalue percentage of variance cumulative percentage of variance

5
## comp 1 3.09276595 51.546099 51.54610
## comp 2 1.57902121 26.317020 77.86312
## comp 3 0.59150386 9.858398 87.72152
## comp 4 0.35035819 5.839303 93.56082
## comp 5 0.30606497 5.101083 98.66190
## comp 6 0.08028581 1.338097 100.00000

fviz_eig(acp.res.log)

Scree plot

50
Percentage of explained variances

40

30

20

10

1 2 3 4 5 6
Dimensions

# Qualités de représentation des eaux

rowSums(acp.res.log$ind$cos2[,1:2])

## Aix Badoit Beckeris Cayranne

## 0.8776607 0.8214127 0.5677568 0.1434473
## Chambon Cristalline Contrex Evian
## 0.2597497 0.6475559 0.9112307 0.5973527
## Ferita Hepar Laurier Ogeu
## 0.0222564 0.9039158 0.7700569 0.6235305
## Ondine Perrier Ribes St_Cyr
## 0.8978507 0.2494251 0.8655147 0.4784917
## St_Hyppo St_Yorre Spa Thonon
## 0.5633011 0.9018862 0.6520556 0.5642958
## Veri Vichy_Celestins Viladrea Vittel
## 0.7547235 0.7188471 0.9352480 0.8379097
## Volvic
## 0.8145571

6
# Projection des variables et des individus sur le premier plan factoriel.
# Ne sont représentées ici que les eaux dont le cos2 (qui mesure la qualité de
# représentation sur le plan factoriel) est supérieur à 0.5.
fviz_pca_var(acp.res.log,axes=c(1,2),repel=T)

Variables − PCA
1.0

NA.
CL
0.5
Dim2 (26.3%)

0.0
HCO3

SO4
CA
−0.5
MG

−1.0

−1.0 −0.5 0.0 0.5 1.0

Dim1 (51.5%)

fviz_pca_ind(acp.res.log,axes=c(1,2),select.ind = list(cos2 = 0.5))

7
 Individuals − PCA
St_Yorre
3
Vichy_Celestins

Ogeu
St_Hyppo
Dim2 (26.3%)

1
Viladrea
Spa Volvic Laurier Badoit
Cristalline
Ribes
Veri
0
Beckeris
Thonon
−1 Ondine
Evian
Vittel
Aix
−2 Hepar
Contrex

−2 0 2 4
Dim1 (51.5%)

5. Combien de composantes principales retiendriez-vous? Justifiez.

D’après le critère de Kaiser qui nous indiquerait de conserver les composantes portant plus de 1/6=16% de la
variabilité totale des données, nous choisirions deux composantes. D’après le critère du coude appliqué à la
représentation graphique des inerties portées par chaque axe, nous pourrions choisir trois composantes. Dans
la mesure où les deux premières composantes principales portent une part très importante de la variabilité
totale des données (77.86%), conserver deux axes semble raisonnable.

6. En quoi la représentation des variables sur le cercle des corrélations illustre-t-elle la matrice des cor-
rélations?

Toutes les variables sont bien représentées sur le cercle des corrélations ce qui autorise à utiliser les positions
des variables les unes par rapport aux autres pour décrire les liens qui existent entre elles. On distingue
trois groupes de variables : i) log(NA) et log(CL), ii) log(HCO3) et iii) les trois autres log-teneurs. Les
groupes i) et iii) sont peu liés. En effet, les coefficients de corrélations entre log(NA) et log(SO4), log(CA)
et log(MG) sont très faibles (moins de 0.2 dans les trois cas). De même, les coefficients de corrélations
entre log(CL) et log(SO4), log(CA) et log(MG) sont très faibles variables (moins de 0.3). log(HCO3) est
modérément corrélé aux variables des groupes i) et iii) (coefficients de corrélations de l’ordre de 0.4 ou 0.5).
Dans i), log(NA) et log(CL) sont très fortement et positivement corrélées entre elles(0.86). De même dans
iii) les trois log-teneurs sont corrélées positivement.

7. Comment interprétez-vous les composantes principales retenues?

On utilise la position des variables (bien représentées) par rapport aux axes. L’axe 1 oppose des eaux riches
en MG, CA, HCO3 à des eaux plus pauvres en ces trois éléments. L’axe 2 oppose des eau riches en sel (NA
et CL) à des eaux moins riches en sel.

8
 8. Interprétez la position des eaux sur le premier plan factoriel.

On prend soin de n’interpréter que les positions des eaux les mieux représentées sur le plan (cf calcul des
cos2). On utilise la signification de chaque composante principale pour interpréter la position des eaux. St
Yorre, Vichy-Célestins, Badoit sont des eaux très riches en bicarbonates de sodium et fortement minéralisées
d’où leur position sur le quart supérieur droit du premier plan factoriel. Ondine, Vittel, Contrex et Hépar
sont très fortement minéralisées et sont en particulier riches en calcium et magnésium mais sont pauvres
en sodium ce qui explique leur représentation dans le quart inférieur droit du premier plan factoriel. Très à
gauche du plan sont représentées trois eaux dont Volvic qui sont toutes les trois très faiblement minéralisées.
Les autres eaux, plus proches de l’origine du repère, ont un profil plutôt intermédiaire, ni très riches, ni très
pauvres en les 6 éléments mesurés.

9. Proposez une classification des eaux par la méthode de classification ascendante hiérarchique (CAH).

a. Combien de classes retenez-vous?

b. Représentez les classes obtenues sur le premier plan factoriel de l’ACP et donnez-leur une inter-
prétation.

#centrer-reduire
tab<-scale(eaux[,1:6], center=TRUE,scale=TRUE)
# Calcul des distances entre individus (au carré)
data.dist <- dist(tab,method = "euclidean")
#---------------------------------------------
# CAH avec distance de Ward
data.ward <- hclust(data.dist,method="ward.D2")
#---------------------------------------------
# Representation graphique
plot(data.ward,hang=-1)
# hang pour aligner les identifiants
# on choisit de distinguer 3 groupes:
rect.hclust(data.ward,k=3)

9
 Cluster Dendrogram
8
4
Height

Contrex
Hepar
Ferita
Thonon
Cayranne
Aix
Evian
Spa
Viladrea
Volvic
Perrier
Ogeu
Beckeris
Chambon
St_Hyppo
Laurier
St_Cyr
Cristalline
Ribes
Veri
Badoit
Ondine
Vittel
St_Yorre
Vichy_Celestins
data.dist
hclust (*, "ward.D2")

#---------------------------------------------
# Representation des individus avec classif cah
# assignation classes sur vecteur individus :
cluster <- cutree(data.ward, k = 3)
eaux$ClassesCAH <- cluster

L’allure du dendrogramme suggère le choix de 3 groupes. Une classification en 3 groupes correspond à un

grand saut d’inertie sur le dendrogramme.

fviz_pca_ind(acp.res.log,axes=c(1,2),select.ind = list(cos2 = 0.6),habillage=as.factor(cluster))

10
 Individuals − PCA
St_Yorre
3
Vichy_Celestins

Ogeu
Dim2 (26.3%)

1 Groups
Viladrea
Spa Volvic Laurier Badoit a 1

Cristalline
Ribes a 2
0 Veri
a 3

−1 Ondine

Vittel
Aix
−2 Hepar
Contrex

−2 0 2 4
Dim1 (51.5%)

by(eaux$HCO3,eaux$ClassesCAH,mean)

## eaux$ClassesCAH: 1
## [1] 300.4286
## --------------------------------------------------------
## eaux$ClassesCAH: 2
## [1] 403
## --------------------------------------------------------
## eaux$ClassesCAH: 3
## [1] 2253

by(eaux$SO4,eaux$ClassesCAH,mean)

## eaux$ClassesCAH: 1
## [1] 42.28571
## --------------------------------------------------------
## eaux$ClassesCAH: 2
## [1] 1333
## --------------------------------------------------------
## eaux$ClassesCAH: 3
## [1] 156

by(eaux$CL,eaux$ClassesCAH,mean)

## eaux$ClassesCAH: 1

11
## [1] 14.90476
## --------------------------------------------------------
## eaux$ClassesCAH: 2
## [1] 10.5
## --------------------------------------------------------
## eaux$ClassesCAH: 3
## [1] 278.5

by(eaux$CA,eaux$ClassesCAH,mean)

## eaux$ClassesCAH: 1
## [1] 82.85714
## --------------------------------------------------------
## eaux$ClassesCAH: 2
## [1] 520.5
## --------------------------------------------------------
## eaux$ClassesCAH: 3
## [1] 96.5

by(eaux$MG,eaux$ClassesCAH,mean)

## eaux$ClassesCAH: 1
## [1] 15.33333
## --------------------------------------------------------
## eaux$ClassesCAH: 2
## [1] 97
## --------------------------------------------------------
## eaux$ClassesCAH: 3
## [1] 10.5

by(eaux$NA.,eaux$ClassesCAH,mean)

## eaux$ClassesCAH: 1
## [1] 16.7619
## --------------------------------------------------------
## eaux$ClassesCAH: 2
## [1] 11.5
## --------------------------------------------------------
## eaux$ClassesCAH: 3
## [1] 1437.5

Pour interpréter les classes obtenues par classification ascendante hiérarchique, on peut représenter les classes
sur le premier plan factoriel de l’ACP et utiliser les positions des eaux dans chaque classe. On peut également
calculer les teneurs moyennes en chacun des 6 éléments dans chaque classe. Ici le groupe 3 comprend les
eaux les plus riches en bicarbonate et sodium. Le groupe 2 correspond aux eaux les plus minéralisées, ie
riches en calcium et magnésium. Le groupe 1 correspond aux autres eaux.

10. Proposez une classification des eaux par la méthode des k-means en utilisant le même nombre de classes
qu’avec la CAH.

a. La nouvelle classification est-elle identique à la précédente?

b. Interprétez les classes obtenues par la méthode des k-means.

12
res.k.means <- kmeans(eaux[,1:6],centers=3)
fviz_pca_ind(acp.res.log,axes=c(1,2),select.ind = list(cos2 = 0.6),
habillage=as.factor(res.k.means$cluster))

Individuals − PCA
St_Yorre
3
Vichy_Celestins

Ogeu
Dim2 (26.3%)

1 Groups
Viladrea
Spa Volvic Laurier Badoit a 1

Cristalline
Ribes a 2
0 Veri
a 3

−1 Ondine

Vittel
Aix
−2 Hepar
Contrex

−2 0 2 4
Dim1 (51.5%)

## Comparaison des deux classifications

table(cluster,res.k.means$cluster)

##
## cluster 1 2 3
## 1 0 0 21
## 2 0 2 0
## 3 1 0 1

Les k-means définissent ici trois classes très proches de celles obtenues par classification ascendante hiérar-
chique. La seule différence porte sur le classement de l’eau Vichy-Célestins. On interprète donc les trois
classes de la même façon que précédement.

# Caractérisation de la classification
res.k.means$totss

## [1] 24808384

res.k.means$withinss

## [1] 0.0 45363.5 2805049.9

13
res.k.means$tot.withinss

## [1] 2850413

res.k.means$betweenss

## [1] 21957971

On peut récupérer (cf script ci-dessus) la répartition de l’inertie correspondant à cette classification. Avec
une inertie inter-classes de 21957971 pour une inertie totale de 24808384, cette classification explique 89%
de la variabilité des données. D’autre part, nous aurions pu nous servir du graphe suivant pour choisir un
nombre de classes en recherchant un coude dans la représentation du rapport inertie inter-classes/inertie
totale en fonction du nombre de classes. Ce critère nous aurait amené à choisir 3 classes.

inertie.expl <- rep(0,times=10)

for (k in 2:10){
clus <- kmeans(tab,centers=k,nstart=5)
inertie.expl[k] <- clus$betweenss/clus$totss }
plot(1:10,inertie.expl,type="b",xlab="Nb. de groupes",ylab="% inertie expliquée")
1.0
0.8
% inertie expliquée

0.6
0.4
0.2
0.0

2 4 6 8 10

Nb. de groupes
Remarque : il aurait été possible d’utiliser le logarithme des 6 teneurs pour construire les deux classifications.
En utilisant le logarithme des 6 teneurs, nous aurions eu plus de facilité à distinguer les eaux faiblement
minéralisées des autres eaux (pour les mêmes raisons que pour l’ACP). Nous aurions obtenu les résultats
ci-dessous (seuls les résultats de la CAH sont fournis).

tab<-scale(log(eaux[,1:6]), center=TRUE,scale=TRUE)
data.dist <- dist(tab,method = "euclidean")
data.ward <- hclust(data.dist,method="ward.D2")
plot(data.ward,hang=-1)
rect.hclust(data.ward,k=4)

14
 Cluster Dendrogram
0 2 4 6 8
Height

Spa
Viladrea
Volvic
Aix
Evian
Cayranne
Ferita
Thonon
Ribes
Veri
Beckeris
Cristalline
Perrier
Ogeu
Laurier
St_Hyppo
Chambon
St_Cyr
Contrex
Hepar
Ondine
Vittel
Vichy_Celestins
Badoit
St_Yorre
data.dist
hclust (*, "ward.D2")

cluster <- cutree(data.ward, k = 4)

eaux$ClassesCAH <- cluster
fviz_pca_ind(acp.res.log,axes=c(1,2),select.ind = list(cos2 = 0.6),habillage=as.factor(cluster))

15
 Individuals − PCA
St_Yorre
3
Vichy_Celestins

Ogeu Groups
Dim2 (26.3%)

1
Viladrea
Spa a 1
Volvic Laurier Badoit a 2
Cristalline
Ribes
Veri a 3
0
a 4

−1 Ondine

Vittel
Aix
−2 Hepar
Contrex

−2 0 2 4
Dim1 (51.5%)

16