Analyse en Composantes Principales sous Python

factoriels. Nous utilisons la librairie Pandas. La vérification de la version de Pandas est

importante. Certaines options de read_excel() sont susceptibles de modifications.

#modification du dossier de travail

import os

os.chdir("… votre dossier de travail …")

#librairie pandas

import pandas

#version

print(pandas.  version  ) # 0.23.0

#chargement de la première feuille de données

X = pandas.read_excel("autos_acp_pour_python.xlsx",sheet_name=0,header=0,index_col=0)

Nous remarquons que :

 Le fichier est un classeur Excel nommé « autos_acp_pour_python.xlsx » ;

 Les données actives sont situées dans la première feuille (sheet_name = 0) ;

 La première ligne correspond aux noms des variables (header = 0)

 La première colonne aux identifiants des observations (index_col = 0).

Nous affichons la dimension de la matrice, nous récupérons le nombre d’observations (n =

18) et de variables (p = 6), enfin nous affichons les valeurs mêmes.

#dimension

print(X.shape) # (18, 6)

#nombre d'observations

n = X.shape[0]

#nombre de variables

p = X.shape[1]

#affichage des données

print(X)

CYL PUISS LONG LARG POIDS V_MAX

Modele

Alfasud TI

Audi 100

Simca 1300

1350

1588

1294

79 393

161

870

85

468

177 1110

68 424

168 1050

Citroen GS Club 1222

59 412

161

930

Fiat 132

Lancia Beta

Peugeot 504

Renault 16 TL

Renault 30

1585

1297

1796

1565

2664

98 439

164 1105

82 429

169 1080

79

449

169 1160

55 424

163 1010

128 452

173 1320

Toyota Corolla

1166

55 399

157

815

Alfetta-1.66

Princess-1800

Datsun-200L

Taunus-2000

Rancho

Mazda-9295

Opel-Rekord

Lada-1300

1570

1798

1998

1993

1442

1769

1979

1294

109 428

162 1060

82 445

172 1160

115 469

169 1370

98 438

170 1080

80 431

166 1129

83

440

165 1095

100 459

173 1120

68 404

161

955

165

160

152

151

165

160

154

140

180

140

175

158

160

167

144

165

173

140

2.2 Préparation des données

Nous devons explicitement centrer et réduire les variables pour réaliser une ACP normée

avec PCA. Nous utilisons la classe StandardScaler pour ce faire. Ici aussi, il est important de

vérifier la version de « scikit-learn » utilisée.

#scikit‐learn

import sklearn

#vérification de la version

print(sklearn.  version  ) # 0.19.1

Nous instancions l’objet et nous l’appliquons sur la matrice X. Nous obtenons une matrice Z

Où (cid:1876)̅(cid:1862) est la moyenne de la variable Xj, (cid:2026)(cid:1862)son écart‐type.

#classe pour standardisation

from sklearn.preprocessing import StandardScaler

#instanciation

sc = StandardScaler()

#transformation – centrage‐réduction

Z = sc.fit_transform(X)

print(Z)

[[-0.77509889 -0.28335818 -1.88508077 -1.09734528 -1.56900676 0.56976043]

[-0.12016326 0.01963869 1.60580955 2.0010414

0.23416142 0.14597168]

[-0.92920139 -0.83885242 -0.44217944 0.25819889 -0.21663062 -0.53209032]

[-1.12733318 -1.29334771 -1.00072189 -1.09734528 -1.11821472 -0.61684807]

[-0.12841875 0.67613189 0.25599862 -0.51639778 0.19659542 0.56976043]

[-0.9209459 -0.13185975 -0.20945342 0.45184806 0.0087654

0.14597168]

[ 0.45221746 -0.28335818 0.72145067 0.45184806 0.60982146 -0.36257482]

[-0.18345536 -1.49534562 -0.44217944 -0.71004695 -0.51715865 -1.54918332]

[ 2.84080623 2.19111619 0.86108628 1.22644473 1.81193359 1.84112668]

[-1.28143568 -1.49534562 -1.60580955 -1.87194195 -1.98223281 -1.54918332]

[-0.16969621 1.23162613 -0.25599862 -0.90369611 -0.14149861 1.41733793]

[ 0.45772112 -0.13185975 0.53526985 1.03279556 0.60982146 -0.02354382]

[ 1.0080872

1.53462299 1.65235475 0.45184806 2.18759363 0.14597168]

[ 0.99432805 0.67613189 0.20945342 0.64549722 0.0087654

0.73927593]

[-0.5219305 -0.2328587 -0.11636301 -0.12909944 0.37691224 -1.21015232]

[ 0.37791804 -0.08136027 0.30254383 -0.32274861 0.12146341 0.56976043]

[ 0.95580242 0.77713084 1.18690271 1.22644473 0.30929343 1.24782243]

[-0.92920139 -0.83885242 -1.37308353 -1.09734528 -0.9303847 -1.54918332]]

Vérifions, par acquit de conscience, les propriétés du nouvel ensemble de données. Les

moyennes sont maintenant nulles (aux erreurs de troncature près) :

#vérification ‐ librairie numpy

import numpy

#moyenne

print(numpy.mean(Z,axis=0))

[-2.22044605e-16 -1.41861831e-16 0.00000000e+00 1.86270752e-15

5.73615229e-16 5.55111512e-16]

Et les écarts-type unitaires.

#écart‐type

print(numpy.std(Z,axis=0,ddof=0))

[1. 1. 1. 1. 1. 1.]

Nous sommes maintenant parés pour lancer l’ACP.

2.3 Analyse en composantes principales avec PCA de ‘’scikit-learn’’

Instanciation et lancement des calculs

Il faut instancier l’objet PCA dans un premier temps, nous affichons ses propriétés.

#classe pour l'ACP

from sklearn.decomposition import PCA

#instanciation

acp = PCA(svd_solver='full')

#affichage des paramètres

print(acp)

PCA(copy=True, iterated_power='auto', n_components=None, random_state=None,

svd_solver='full', tol=0.0, whiten=False)

Le paramètre (svd_solver = ‘full’) indique l’algorithme utilisé pour la décomposition en valeurs

singulières. Nous choisissons la méthode « exacte », sélectionnée de toute manière par défaut

pour l’appréhension des bases de taille réduite. D’autres approches sont disponibles pour le

traitement des grands ensembles de données. Le nombre de composantes (K) n’étant pas

spécifié (n_components = None), il est par défaut égal au nombre de variables (K = p).

Nous pouvons lancer les traitements dans un second temps. La fonction fit_transform()

renvoie en sortie les coordonnées factorielles Fik que nous collectons dans la variable coord.

Nous affichons le nombre de composantes générées (K), il est bien égal à p = 6.

#calculs

coord = acp.fit_transform(Z)

#nombre de composantes calculées

print(acp.n_components_) # 6

Valeurs propres et scree plot

(valeurs propres, λk) associées aux axes factoriels.

#variance expliquée

print(acp.explained_variance_)

[4.68090853 0.90641889 0.39501114 0.22650574 0.09826011 0.04583676]

#valeur corrigée

eigval = (n‐1)/n*acp.explained_variance_

print(eigval)

[4.42085806 0.85606229 0.37306608 0.21392209 0.09280121 0.04329027]

Là, tout rentre dans l’ordre

Nous aurions pu obtenir les bonnes valeurs propres en passant par les valeurs singulières

.singular_values_ issues de la factorisation de la matrice des données centrées et réduites

#ou bien en passant par les valeurs singulières

print(acp.singular_values_**2/n)

[4.42085806 0.85606229 0.37306608 0.21392209 0.09280121 0.04329027]

PCA fournit également les proportions de variance associées aux axes. Il n’est pas nécessaire

d’effectuer une correction dans ce cas.

#proportion de variance expliquée

print(acp.explained_variance_ratio_)

[0.73680968 0.14267705 0.06217768 0.03565368 0.01546687 0.00721505]

La première composante accapare 73.68% de l’information disponible. Il y a un fort « effet

taille » dans nos données . Nous disposons de 87.94% avec les deux premiers facteurs. Les

suivants semblent anecdotiques.

Nous disposons des éléments permettant de construire le graphique ‘’Scree plot’’ (éboulis des

valeurs propres) (Figure 2).

#scree plot

plt.plot(numpy.arange(1,p+1),eigval)

plt.title("Scree plot")

plt.ylabel("Eigen values")

plt.xlabel("Factor number")

plt.show()

Figure 2 - Scree plot

Le graphique du cumul de variance restituée selon le nombre de facteurs peut être

intéressant également (Figure 3).

#cumul de variance expliquée

plt.plot(numpy.arange(1,p+1),numpy.cumsum(acp.explained_variance_ratio_))

plt.title("Explained variance vs. # of factors")

plt.ylabel("Cumsum explained variance ratio")

plt.xlabel("Factor number")

plt.show()

Figure 3 - Variance expliquée vs. Nombre de facteurs

Détermination du nombre de facteur à retenir

Les ‘’cassures’’ dans les graphiques ci-dessus (Figure 2, Figure 3) sont souvent évoquées

(règle du coude) pour identifier le nombre de facteurs K à retenir. La solution (K = 2)

semble s’imposer ici.

Représentation des individus – Outils pour l’interprétation

Coordonnées factorielles. Les coordonnées factorielles (Fik) des individus ont été

collectées dans la variable coord . Nous les positionnons dans le premier plan

factoriel avec leurs labels pour situer et comprendre les proximités entre les véhicules.

Je ferai deux commentaires au préalable :

1. L’ajout d’une étiquette dans un graphique nuage de points n’est pas très pratique sous

Python (librairie Matplotlib),

2. Les outils graphiques calculent souvent automatiquement les échelles en fonction des

plages de valeurs. Ce n’est pas une bonne idée en ce qui concerne l’ACP. En effet, les axes

n’ont pas la même importance (% de variance restituée). Pour ne pas fausser la

perception des proximités, il est très important de veiller à ce que les échelles soient

identiques en abscisse et en ordonnée

#positionnement des individus dans le premier plan

fig, axes = plt.subplots(figsize=(12,12))

axes.set_xlim(‐6,6) #même limites en abscisse

axes.set_ylim(‐6,6) #et en ordonnée

#placement des étiquettes des observations

for i in range(n):

plt.annotate(X.index[i],(coord[i,0],coord[i,1]))

#ajouter les axes

plt.plot([‐6,6],[0,0],color='silver',linestyle='‐',linewidth=1)

plt.plot([0,0],[‐6,6],color='silver',linestyle='‐',linewidth=1)

#affichage

plt.show()

Représentation des individus dans le premier plan factoriel

Qualité de représentation – Les COS² (cosinus carré).

Pour calculer la qualité de représentation des individus sur les axes, nous devons d’abord

calculer les carrés des distances à l’origine des individus, qui correspondent également à

leur contribution dans l’inertie totale

#contribution des individus dans l'inertie totale

di = numpy.sum(Z**2,axis=1)

print(pandas.DataFrame({'ID':X.index,'d_i':di}))

Concrètement, la Renault 30 et la Toyota Corolla sont les deux véhicules qui se démarquent

le plus des autres, et on les retrouve aux deux extrémités du premier axe factoriel qui porte

73.68% de l’information disponible (Figure 4).

0

1

2

ID

d_i

Alfasud TI

8.225176

Audi 100

6.673755

Simca 1300

2.159327

3 Citroen GS Club

6.780145

4

5

6

7

8

9

10

11

12

13

14

15

16

17

Fiat 132

1.169124

Lancia Beta

1.134950

Peugeot 504

1.512793

Renault 16 TL

5.636826

Renault 30 21.789657

Toyota Corolla 16.290143

Alfetta-1.66

4.456770

Princess-1800

1.952513

Datsun-200L 11.112624

Taunus-2000

2.452986

Rancho

1.963373

Mazda-9295

0.684521

Opel-Rekord

6.083119

Lada-1300

7.922198

Nous pouvons alors déduire la qualité de représentation des individus sur l’axe n°k avec :

#qualité de représentation des individus ‐ COS2

cos2 = coord**2

for j in range(p):

(cid:0)

cos2[:,j] = cos2[:,j]/di

print(pandas.DataFrame({'id':X.index,'COS2_1':cos2[:,0],'COS2_2':cos2[:,1]}))

Les COS² pour les deux premiers facteurs sont affichés [ACP, page 30 ; TUTO_R, page 9]

id

COS2_1

COS2_2

Alfasud TI 0.556218 0.387670

Audi 100 0.365334 0.349406

Simca 1300 0.580284 0.210694

Citroen GS Club 0.976992 0.001879

Fiat 132 0.156579 0.413826

Lancia Beta 0.081555 0.033900

Peugeot 504 0.309202 0.575488

Renault 16 TL 0.673539 0.170535

Renault 30 0.892431 0.051920

Toyota Corolla 0.975219 0.003426

Alfetta-1.66 0.042978 0.820652

Princess-1800 0.530947 0.362855

Datsun-200L 0.778390 0.028137

Taunus-2000 0.704819 0.096496

Rancho 0.243273 0.410469

Mazda-9295 0.217336 0.185337

Opel-Rekord 0.861900 0.001790

Lada-1300 0.926052 0.002607

0

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

Conformément à la théorie, pour chaque individu, la somme des COS² sur l’ensemble des

facteurs est égale à 1.

#vérifions la théorie ‐ somme en ligne des cos2 = 1

print(numpy.sum(cos2,axis=1))

[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]

Contribution des individus aux axes (CTR). Elles permettent de déterminer les

individus qui pèsent le plus dans la définition de chaque facteur.

#contributions aux axes

ctr = coord**2

0

for j in range(p):

1

id

Alfasud TI 0.057493 0.206933

Audi 100 0.030640 0.151329

CTR_1

CTR_2

ctr[:,j] = ctr[:,j]/(n*eigval[j])

Simca 1300 0.015746 0.029525

2

print(pandas.DataFrame({'id':X.index,'CTR_1':ctr[:,0],'CTR_2':ctr[:,1]}))

3 Citroen GS Club 0.083244 0.000827

4

5

6

7

8

9

10

11

12

13

14

15

16

17

Fiat 132 0.002300 0.031398

Lancia Beta 0.001163 0.002497

Peugeot 504 0.005878 0.056499

Renault 16 TL 0.047711 0.062384

Renault 30 0.244369 0.073419

Toyota Corolla 0.199640 0.003622

Alfetta-1.66 0.002407 0.237357

Princess-1800 0.013028 0.045978

Datsun-200L 0.108701 0.020292

Taunus-2000 0.021727 0.015361

Rancho 0.006002 0.052300

Mazda-9295 0.001870 0.008233

Opel-Rekord 0.065888 0.000707

Lada-1300 0.092194 0.001340

Sans surprises, ce sont la Renault 30 et la Toyota Corolla qui sont déterminants pour le

premier axe ; pour le second, nous avons l’Alfetta-1.66, l’Alfasud TI et l’Audi 100.

Les sommes en ligne sont égales à l’unité ici :

#vérifions la théorie

print(numpy.sum(ctr,axis=0))

[1. 1. 1. 1. 1. 1.]

Représentation des variables – Outils pour l’aide à l’interprétation

Nous avons besoin des vecteurs propres pour l’analyse des variables. Ils sont fournis par le

champ .components_

#le champ components_ de l'objet ACP

print(acp.components_)

[[ 0.42493602 0.42179441 0.42145993 0.38692224 0.43051198 0.35894427]

[ 0.12419108 0.41577389 -0.41181773 -0.446087

-0.24267581 0.6198626 ]

[-0.35361252 -0.18492049 0.06763394 0.60486812 -0.48439601 0.48547226]

[ 0.80778648 -0.35779199 -0.27975231 0.21156941 -0.30171136 -0.0735743 ]

[ 0.15158003 -0.29373465 0.73056903 -0.47819008 -0.30455842 0.18865511]

[-0.05889517 -0.63303302 -0.19029153 -0.10956624 0.5808122

0.45852167]]

#racine carrée des valeurs propres

sqrt_eigval = numpy.sqrt(eigval)

#corrélation des variables avec les axes

corvar = numpy.zeros((p,p))

for k in range(p):

corvar[:,k] = acp.components_[k,:] * sqrt_eigval[k]

#afficher la matrice des corrélations variables x facteurs

print(corvar)

Les variables sont maintenant en ligne, les facteurs en colonne :

[[ 0.89346354 0.1149061 -0.21598347 0.37361508 0.04617627 -0.01225391]

[ 0.88685803 0.38468911 -0.11294784 -0.16548492 -0.08948124 -0.13171084]

[ 0.88615477 -0.38102873 0.04131023 -0.12939024 0.22255537 -0.03959265]

[ 0.81353638 -0.4127359

0.36944822 0.09785447 -0.14567244 -0.0227967 ]

[ 0.90518746 -0.22453248 -0.29586489 -0.13954667 -0.09277852 0.12084561]

[ 0.75471037 0.57351941 0.29652226 -0.03402937 0.05747056 0.09540146]]

Si l’on s’en tient spécifiquement aux deux premiers facteurs :

#on affiche pour les deux premiers axes

print(pandas.DataFrame({'id':X.columns,'COR_1':corvar[:,0],'COR_2':corvar[:,1]}))

id

COR_1

COR_2

0

CYL 0.893464 0.114906

1 PUISS 0.886858 0.384689

2

3

LONG 0.886155 -0.381029

LARG 0.813536 -0.412736

4 POIDS 0.905187 -0.224532

5 V_MAX 0.754710 0.573519

Nous pouvons dessiner maintenant le cercle des corrélations (Figure 5).

#cercle des corrélations

fig, axes = plt.subplots(figsize=(8,8))

axes.set_xlim(‐1,1)

axes.set_ylim(‐1,1)

#affichage des étiquettes (noms des variables)

for j in range(p):

plt.annotate(X.columns[j],(corvar[j,0],corvar[j,1]))

#ajouter les axes

plt.plot([‐1,1],[0,0],color='silver',linestyle='‐',linewidth=1)

plt.plot([0,0],[‐1,1],color='silver',linestyle='‐',linewidth=1)

#ajouter un cercle

cercle = plt.Circle((0,0),1,color='blue',fill=False)

axes.add_artist(cercle)

#affichage

plt.show()

Cercle des corrélations

On perçoit clairement l’effet taille sur le premier axe : les voitures puissantes et rapides sont

aussi les plus lourdes et imposantes, la relation globale entre les variables est en réalité

déterminée par la cylindrée (CYL).

3 Traitement des individus et variables illustratifs

3.1 Individus supplémentaires

Nous souhaitons positionner deux véhicules supplémentaires, des Peugeot, par rapport aux

existantes [TUTO_R, page 14].

Modele

Peugeot 604

Peugeot 304 S

CYL

PUISS

LONG

LARG

POIDS

V_MAX

2664

1288

136

74

472

414

177

157

1410

915

180

160

Figure 6 - Individus supplémentaires

Nous les chargeons avec read_excel() de Pandas, elles sont situées dans le seconde feuille

du classeur Excel (sheet_name = 1).

#chargement des individus supplémentaires

indSupp = pandas.read_excel("autos_acp_pour_python.xlsx",sheet_name=1,header=0,index_col=0)

print(indSupp)

CYL PUISS LONG LARG POIDS V_MAX

Modele

Peugeot 604

2664

Peugeot 304 S 1288

472

177

1410

414

157

915

180

160

Nous devons centrer et réduire les variables des individus supplémentaires à l’aide des

paramètres (moyennes et écarts-type) des données actives ayant servi à construire le repère

factoriel.

#centrage‐réduction avec les paramètres des individus actifs

ZIndSupp = sc.transform(indSupp)

print(ZIndSupp)

[[ 2.84080623 2.59511201 1.79199036 2.0010414

2.48812166 1.84112668]

[-0.94571238 -0.53585556 -0.90763148 -1.87194195 -1.23091273 0.14597168]]

Il ne reste plus qu’à faire calculer par la fonction .transform() leurs coordonnées.

#projection dans l'espace factoriel

coordSupp = acp.transform(ZIndSupp)

print(coordSupp)

[[ 5.56329226 0.33860928 -0.46428878 0.40214608 -0.38981076 -0.08102064]

[-2.21224139 1.25777905 -0.09304388 -0.35370189 0.648528 0.12473042]]

Et à les représenter dans le premier plan factoriel parmi les observations actives.

#positionnement des individus supplémentaires dans le premier plan

fig, axes = plt.subplots(figsize=(12,12))

axes.set_xlim(‐6,6)

axes.set_ylim(‐6,6)

#étiquette des points actifs

for i in range(n):

plt.annotate(X.index[i],(coord[i,0],coord[i,1]))

#étiquette des points supplémentaires (illustratifs) en bleu ‘b’

for i in range(coordSupp.shape[0]):

plt.annotate(indSupp.index[i],(coordSupp[i,0],coordSupp[i,1]),color='b')

#ajouter les axes

plt.plot([‐6,6],[0,0],color='silver',linestyle='‐',linewidth=1)

plt.plot([0,0],[‐6,6],color='silver',linestyle='‐',linewidth=1)

#affichage

plt.show()

Figure 7 - Positionnement des individus supplémentaires dans le premier plan factoriel

La Peugeot 604 se rapproche plutôt de la Renault 30, la Peugeot 304 de l’Alfasud TI. Pour qui

connaît un peu les voitures de ces années-là, tout cela est parfaitement cohérent.