Cours Analyse des données : Analyse en Composantes Principales (ACP)
Dhafer Malouche
Ecole Supérieure de la Statistique et de l’Analyse de l’Information
Ecole Supérieure de la Statistique et de l’Analyse de l’Information
Janvier, 2016
- Introduction et problématique
- Calcul des Composantes Principales (ACP).
- Mise en oeuvre sous
R.
Website: http://essai.academia.edu/DhaferMalouche
email : dhafer.malouche@yale.edu
Introduction et problématique
Exemples de données,
Exemple 1 decathlon
- JO athènes 23/24 Août 2004
Decastar 25/26 Septembre 2004
Les variables
10 scores des athlètes (épreuve du décathlon)
Total des points, classement des athlètes
Variable compétition.
12 Variables quantitatives, 1 Variable qualitative.
> library(FactoMineR)
> data(decathlon)
> head(decathlon)
100m Long.jump Shot.put High.jump 400m 110m.hurdle Discus
SEBRLE 11.04 7.58 14.83 2.07 49.81 14.69 43.75
CLAY 10.76 7.40 14.26 1.86 49.37 14.05 50.72
KARPOV 11.02 7.30 14.77 2.04 48.37 14.09 48.95
BERNARD 11.02 7.23 14.25 1.92 48.93 14.99 40.87
YURKOV 11.34 7.09 15.19 2.10 50.42 15.31 46.26
WARNERS 11.11 7.60 14.31 1.98 48.68 14.23 41.10
Pole.vault Javeline 1500m Rank Points Competition
SEBRLE 5.02 63.19 291.7 1 8217 Decastar
CLAY 4.92 60.15 301.5 2 8122 Decastar
KARPOV 4.92 50.31 300.2 3 8099 Decastar
BERNARD 5.32 62.77 280.1 4 8067 Decastar
YURKOV 4.72 63.44 276.4 5 8036 Decastar
WARNERS 4.92 51.77 278.1 6 8030 Decastar
> summary(decathlon)
100m Long.jump Shot.put High.jump
Min. :10.44 Min. :6.61 Min. :12.68 Min. :1.850
1st Qu.:10.85 1st Qu.:7.03 1st Qu.:13.88 1st Qu.:1.920
Median :10.98 Median :7.30 Median :14.57 Median :1.950
Mean :11.00 Mean :7.26 Mean :14.48 Mean :1.977
3rd Qu.:11.14 3rd Qu.:7.48 3rd Qu.:14.97 3rd Qu.:2.040
Max. :11.64 Max. :7.96 Max. :16.36 Max. :2.150
400m 110m.hurdle Discus Pole.vault
Min. :46.81 Min. :13.97 Min. :37.92 Min. :4.200
1st Qu.:48.93 1st Qu.:14.21 1st Qu.:41.90 1st Qu.:4.500
Median :49.40 Median :14.48 Median :44.41 Median :4.800
Mean :49.62 Mean :14.61 Mean :44.33 Mean :4.762
3rd Qu.:50.30 3rd Qu.:14.98 3rd Qu.:46.07 3rd Qu.:4.920
Max. :53.20 Max. :15.67 Max. :51.65 Max. :5.400
Javeline 1500m Rank Points
Min. :50.31 Min. :262.1 Min. : 1.00 Min. :7313
1st Qu.:55.27 1st Qu.:271.0 1st Qu.: 6.00 1st Qu.:7802
Median :58.36 Median :278.1 Median :11.00 Median :8021
Mean :58.32 Mean :279.0 Mean :12.12 Mean :8005
3rd Qu.:60.89 3rd Qu.:285.1 3rd Qu.:18.00 3rd Qu.:8122
Max. :70.52 Max. :317.0 Max. :28.00 Max. :8893
Competition
Decastar:13
OlympicG:28
Exemple 2 : temperature
Ces données peuvent être téléchargées sur le site
[http://www.factominer.free.fr/book/temperature.csv]
> temperature <- read.csv("~/Documents/Teaching/Analyse_des_donnees_2016/temperature.csv", sep=";")
> head(temperature)
X January February March April May June July August September
1 Amsterdam 2.9 2.5 5.7 8.2 12.5 14.8 17.1 17.1 14.5
2 Athens 9.1 9.7 11.7 15.4 20.1 24.5 27.4 27.2 23.8
3 Berlin -0.2 0.1 4.4 8.2 13.8 16.0 18.3 18.0 14.4
4 Brussels 3.3 3.3 6.7 8.9 12.8 15.6 17.8 17.8 15.0
5 Budapest -1.1 0.8 5.5 11.6 17.0 20.2 22.0 21.3 16.9
6 Copenhagen -0.4 -0.4 1.3 5.8 11.1 15.4 17.1 16.6 13.3
October November December Annual Amplitude Latitude Longitude Area
1 11.4 7.0 4.4 9.9 14.6 52.2 4.5 West
2 19.2 14.6 11.0 17.8 18.3 37.6 23.5 South
3 10.0 4.2 1.2 9.1 18.5 52.3 13.2 West
4 11.1 6.7 4.4 10.3 14.4 50.5 4.2 West
5 11.3 5.1 0.7 10.9 23.1 47.3 19.0 East
6 8.8 4.1 1.3 7.8 17.5 55.4 12.3 NorthCette base contient les enregistrements des températures des capitales européennes de Janvier ? Décembre avec les variables + Les coordonnées GPS de chaque ville. + Amplitude thermale : Différence entre les températures maximales et minimales + Moyenne Annuelle. + Une variable qualitative : la direction (S, N, O, E).
Exemple 3 : Données génomiques chicken
> chicken <- read.csv("~/Documents/Teaching/Analyse_des_donnees_2016/chicken.csv", sep=";")
> dim(chicken)
[1] 7406 44Description : 43 poulets ayant subis 6 régimes : régime normal (N), Je?ne pendant 16h (F16), Je?ne pendant 16h et puis se réalimenter pendant 5h (F16R5), (F16R16), (F48), (F48R24)
Variables : Après le régime, on a effectué une analyse des gènes utilisant une puce ADN : 7407 expressions de gènes.
Objectif :
- Voir si les gènes s’expriment différemment selon le niveau de stress.
- Combien de temps lui faut-il le poulet pour revenir à la situation normale ?
Exemple 4 : Données banques
> banques <- read.csv("~/Documents/Teaching/Analyse_des_donnees_2016/banques.csv")
> head(banques)
Classe_PV Age_PV Cadre Non_Cadre Diplome Age_Moy Anc_Moy Surface_
1 3 29 14 10 10 41,71 5,54 12,75
2 3 21 6 14 6 34,2 4,7 11,8
3 5 13 3 5 3 31,13 6,63 12,25
4 3 18 5 7 4 36,58 5,33 16,42
5 4 17 2 5 2 44,86 4,14 27,86
6 4 16 3 3 2 40,33 5 40
Type_concep Nbr_reclam Age_RPV Anc_RPV Anc_RPV_PV Qualt_client
1 A 16 38 21 20 0,37
2 A 7 41 28 9 0,53
3 A 7 46 66 24 0,49
4 A 12 38 28 9 0,69
5 A 2 47 7 7 0,44
6 A 8 37 66 7 0,4
> dim(banques)
[1] 102 14C’est une base de données qui concerne 102 agences bancaires en retenant plusieurs informations sur chaque variable dont une variable qui mesure la performance de chaque agence (point de vente).
L’objectif étant d’expliquer la peformance de chaque point de vente en fonction des ses caractéristiques.
Les données : Individus et variables
\(E=\{w_1,\ldots,w_n\}\) un échantillon de taille \(n\).
\(x^1,\ldots,x^d\), \(d\) variables quantitatives.
\(y\) une variable quantitative supplémentaire, \(y=f(x^1,\ldots,x^d)\), .
\(A\) une variable qualitative supplémentaire à \(r\) modalités : \(a_1,\ldots,a_r\).
les observations :
- quantitatives : \(x_i^j=x^j(w_i)\), \(\forall\, i=1\ldots n\) et \(j=1\ldots d\).
- supplémentaire quantitative : \(y_i=y(w_i)\), \(\forall\, i=1\ldots n\).
- supplémentaire qualitative : \(a_l=A(w_i),\) \(\forall\, i=1\ldots n\) et \(l=1\ldots r\).
\(x^1,\ldots,x^d\) variables actives, \(y\) et \(A\) variables supplémentaires et \(X=(x_i^j)\) est la matrice \(n\times d\) des variables actives.
Les éléments de l’analyse
- \(p_i\) est le poids de chaque individu \(w_i\) dans l’échantillon \(E\) : \[p_i\in [0,1],\; \sum_{i=1}^n p_i=1\].
Matrice des poids \[D=\left(\begin{array}{ccc} p_1 & & 0 \\ & \ddots & \\ 0 & & p_n \\ \end{array}\right)\]
Nuage des individus : chaque individu \(w_i\) est représenté par \(x_i=\,^t(x_i^1,\ldots,x_i^d)\in \mathbb{R}^d\). \(\mathcal{N}_I=\{ \, x_i\in \mathbb{R}^d ,\,i=1\ldots n\}.\)
Nuage des variables : chaque variable \(x^j\) est représentée par \(x^j=\,^t(x_1^j,\ldots,x_n^j)\in \mathbb{R}^n\). \(\mathcal{N}_V=\{ \, x^j\in \mathbb{R}^n ,\,j=1\ldots d\}.\)
Géométries dans \(\mathbb{R}^n\) et \(\mathbb{R}^d\)
Centre de gravité de \(\mathcal{N}_I\), \(g=\,^t(\overline{x}^1,\ldots,\overline{x}^d)\in \mathbb{R}^d\) où \[\overline{x}^j=\sum_{i=1}^n p_ix_i^j \mbox{ est la moyenne de } x^j\]
Pour la suite \(X\) est une matrice \(\Rightarrow\) \(g=0\).
- \(\mathbb{R}^d\) est un espace euclidien muni du produit scalaire de métrique \(M\) : \[<x_i,x_{i'}>_M=\,(x_i)'Mx_{i'}=\sum_{j=1}^d x_i^jx_{i'}^j.\] \(M\) est une matrice diagonale \(d\times d\).
\(\mathbb{R}^n\) est un espace euclidien muni du produit scalaire : \[<x^j,x^{j'}>_M=\,\,(x^j)'\,D\,x^{j'}=\sum_{i=1}^np_i\, x_i^jx_{i}^{j'}.\] C’est un produit scalaire de métrique \(D\).
Variances, Covariances et Corrélations
Variance de la variable \({x}^j\) : \[\mbox{Var}(x^j)=\sum_{i=1}^n p_i (x^j_i)^2\] \(\Rightarrow\) \(\mbox{Var}(x^j)=\parallel x^j \parallel^2\) dans \(\mathbb{R}^n\).
Covariance entre les variables \(x^j\) et \(x^{j'}\) \[\mbox{Cov}(x^j,x^{j'})= \sum_{i=1}^n p_i x_i^jx_{i}^{j'}=<x^j,x^{j'}>.\]
Corrélation entre les variables \({x}^j\) et \({x}^{j'}\) \[\hspace{-.5cm}\mbox{Cor}(x^j,x^{j'})=\displaystyle\frac{\mbox{Cov}(x^j,x^{j'})}{\sqrt{\mbox{Var}(x^j)\,\mbox{Var}(x^{j'})}} =\displaystyle\frac{<x^j,x^{j'}>}{\parallel x^j \parallel\, \parallel x^{j'} \parallel}=\cos \alpha\] où \(\alpha\) est l’angle entre les vecteurs \(x^j\) et \(x^{j'}\) dans \(\mathbb{R}^n\).
Problématique
Trouver un sous-espace vectoriel de dimension faible (2 ou 3) pour faire une représentation des graphiques des variables et des individus.
Cette représentation a pour objectif de bien montrer
les relations existantes entre les variables : les corrélations.
les similarités et les éventuels sous-groupes dans l’échantillon.
Chercher une suite de variables artificielles \(C^1,\ldots,C^k\), \(k<<d\), qui
sont combinaisons linéaires des variables \(x^1,\ldots,x^d\),
non-corrélées entre elles
et captent le maximum de variance du nuage des individus \(\mathcal{N}_I\).
Principe d’une ACP
Calcul des Composantes Principales (ACP).
Objectif
Rappelons qu’un jeu de donnés peut présenté comme une matrice \(X\) de dimension \(n\times p\) composées d’observations de \(p\) variables quantitatives sur un échantillon de taille \(n\).
- Nous allons montrer dans cette section comment
- trouver un espace de dimension faible (2 ou 3) qui permet d’observer la variation de ses données ?
- détecter soit des groupes d’individus homogènes ou quelques individus qui présentant des observations abérantes par rapport au jeu de données ?
- détecter les variables qui sont les plus corrélées entre elles ?
- Pour répondre à toutes ces questions on va se placer dans
- \(\mathbb{R}^p\) muni de la métrique \(M\) si on s’intéresse à une représentation des individus (lignes de \(X\)). Les éléments de cette ensemble seront représentés par des “points”.
- \(\mathbb{R}^n\) muni de la métrique \(D\) si on s’intéresse à une représentation des variables (colonnes de \(X\)). Les éléments de cette ensemble seront représntés par des “vecteurs”.
- Donc on cherchera
- (Prob 1.) un vecteur \(u\in \mathbb{R}^p\) unitaire (de norme égale à 1) tel que la projection du nuage des points \(\left\{x_i,\,i=1,\dots,n\right\}\) soit d’inertie maximale.
- (Prob 2.) un vecteur \(C\in \mathbb{R}^n\) qui soit combinaison linéaire des variables originales \(x^1,\ldots,x^p\) tel que la variance de \(C\) (considérée comme variable artificielle) soit maximale.
Calcul de la première composante
Réponse à (Prob 1.) :
calculons l’inertie du nuage de points projection du nuage des points \(\mathcal{N}_I=\left\{x_i,\,i=1,\dots,n\right\}\) sur une droite vectorielle engendrée par un vecteur unitaire \(u\).
Notons ce nuage par \(\mathcal{N}^\perp(u)\).
- Pour tout \(i=1,\ldots,n\), \(x_i^\perp=<x_i,u>_Mu=\left[(x_i)'Mu\right]u\). Donc \[||x_i^\perp||^2_M=\left((x_i)'Mu\right)^2.\]
- On a \[\mathcal{I}(\mathcal{N}^\perp(u))= \displaystyle\sum_{i=1}^n d_i ||(x_i)^\perp||_M^2=\displaystyle\sum_{i=1}^n d_i \left[(x_i)'Mu\right]^2=u'(XM)'D(XM)u\]
- Donc on aura à chercher le maximum de la fonction suivante \[u\longmapsto u'(XM)'D(XM)u,\; \mbox{ quand }||u||_M=1,\] où \(V=X'DX\). Ce qui est équivalent à chercher le maximum de \[\phi_\lambda\,:\,u\longmapsto u'MVMu-\lambda (u'Mu-1),\; \mbox{ quand }||u||_M=1,\]
- Pour tout \(\lambda\) le gradiant de \(\phi_\lambda\) est \[\nabla \phi_\lambda (u)=2(MVM)u-2Mu\]
- Donc \(\nabla \phi_\lambda (u)=0\) si et seulement si \((MV)(Mu)=\lambda Mu\). Donc si \(v\) est un vecteur propre associé à \(\lambda\) de la matrice \(VM\) alors \(u=M^{-1}v\) est extremum. Comme \(||u||_M=1\). Donc il faut que \[ 1=u'Mu=v'(M^{-1})MM^{-1} v=v'M^{-1}v.\] Donc \(v\) est \(M^{-1}-\)unitaire.
- Comme \(\phi_\lambda(u)=\lambda\), si \(Mu\) est un vecteur propre \(M^{-1}-\)unitaire de \(MV\). Donc \(\lambda\) est la plus grande valeur propre de \(MV\)
- Si \(\lambda_1\geq \lambda_2\geq \ldots \geq \lambda_p\) sont les valeurs propres de \(MV\) (La matrice \(MV\) est définie positive). Si \(v^1\) est un vecteur propre \(M^{-1}-\)unitaire associé à \(\lambda_1\), alors \(u^1=M^{-1}v^1\) un maximum de \(\phi_\lambda\) et \(C^1=XMu^1\) est un vecteur de \(\mathbb{R}^n\) contentant les coordonnées des projections de \(\mathcal{N}\) sur \(\Delta_{u^1}\).
Terminologie
- \(C^1\) est appelée première composante principale
- La droite engendrée par le vecteur \(u^1\) est appelée première axe principal
- \(v^1=Mu^1\) est le premier facteur principal (
loading factor).
Réponse à (Prob 2.) :
Toute variable artificielle \(C\) combinaison linéaire des vecteurs \(x^1,\ldots,x^d\) s’écrit : \[ C=\displaystyle\sum_{i=1}^d v^j x^j=Xv\] ouù \(v=(v^1,\ldots,v^d)'\in \mathbb{R}^d\) et la variance de \(C\) s’exprime \[\mbox{Var}(C)=C'DC=v'X'DXv=v'Vv\]
Exercice Montrer que chercher la variable \(C\) de variance maximale est équivalent au probleème (Prob 1).
Calcul des autres composantes principales
Maintenant nous allons un deuxième axe orthogonal au premier axe principal \(\Delta_{u^1}=<u^1>\) telle que la projection de \(\mathcal{N}\) l’inertie du nuage de points soit maximale. Evidement pour n’importe quel vecteur unitaire \(u\) l’inertie du nuage de points projection du nuage des points \(\mathcal{N}\) sur la droite vectorielle \(\Delta_u\) s’exprime par \[\mathcal{I}(\mathcal{N}^\perp(u))=u'(XM)'D(XM)u\] Ainsi pour estimer le deuxième axe principal il faut chercher le maximum de la fonction de \(\phi_\lambda\) sous les contraines \[u'Mu=1\;\mbox{ et }\; u'Mu^1=0\] Donc la solution est le vecteur \(u^2\), tel que \(v^2=Mu^2\) est le vecteur propre \(M^{-1}-\)orthogonal ? \(v^1\) et associé à la deuxième valeur propre \(\lambda_1\). On \(u^2\) et \(u^1\) sont aussi \(M-\)orthogonaux, en effet, \[(u^2)'Mu^2=(v^2)'M^{-1}MM^{-1}v^2=(v^2)'M^{-1}v^2=0.\] Ainsi la deuxième composante principale \(C^2=XMu^2=Xv^2\)…
On constinue ainsi à calculer toutes les composantes principales, les axes et les facteurs principaux :
- Diagonaliser \(MV\) : \(\lambda_1\geq \ldots \geq \lambda_p\)
- \(v^1,\ldots,v^p\) une base \(M^{-1}-\)orthornomale de \(\mathbb{R}^p\) : les facteurs principaux de l’analyse en composantes principales \((X,D,M)\)
- \(u^1=M^{-1}v^1,\ldots,u^p=M^{-1}v^p\) une base \(M-\)orthornomale de \(\mathbb{R}^p\) : les axes principaux de l’analyse en composantes principales \((X,D,M)\)
- \(C^1=Xv^1=XMu^1,\ldots,C^p=Xv^p=XMu^p\) une famille de vecteurs de \(\mathbb{R}^n\) : les composantes principales.
Représentation graphiques des individus et/ou variables
Les individus
On a déjà vu les coordonnées des indivdus sur le \(k-\)i?me axe principal est donnée dans la \(k-\)ième composante principal. Donc la coordonn?e du \(i-\)i?me est \(C_i^k\).
Les variables
Reconstitution des donn'ees : \(\mathbf{C}=X \mathbf{u}\;\iff\; X= \mathbf{C}\, ^t\mathbf{u}\)
\(\Rightarrow\) \(\forall,\, j=1,\ldots,d\), $x^j=_{k=1}^d u_j^k C^k $.
Comme \(\{C^1,\ldots,C^d\}\) est une famille de vecteurs orthogonaux 2 `a 2 dans \(\mathbb{R}^n\) muni de la m'etrique \(D\), \(\forall\,k\not=k'\), \[\mbox{Cov}(C^k,C^{k'})=<C^k, C^{k'}>=\, ^t(C^k) D (C^{k'})\mbox{ et } \parallel C^k\parallel =\sqrt{\lambda_k}.\] On consid`ere la famille \(\{Z^1,\ldots,Z^d\}\) o? \(\forall\,k=1,\ldots,d\) \[Z^k= \displaystyle\frac{C^k}{\sqrt{\lambda_k}}\] \(<Z^k>\) est appel'e le \(k-\)i`eme axe unitaire.
On a déjà vu les vecteurs \(C^k\) forment une famille \(D-\)orthoganaux. Posons alors pour tout \(k=1,\ldots,p\) les vecteurs de \(\mathbb{R}^n\) \[Z^k=\displaystyle\frac{C^k}{\sqrt{\lambda_k}}\] sont des vecteurs \(D-\)orthonnormés. Ainsi pour une variable \(x^j\) sa coordonnée sur l’axe engendré par \(Z^k\) est donnée par \[<x^j,Z^k>=(x^j)'DZ^k=\mbox{Cor}(x^j,C^k)\times \sqrt{\mbox{Var}(x^j)}.\]
Exercice
Montrez que
les composantes principales sont les vecteurs propres de la matrice \(W=XMX'D\). Quelles sont les valeurs propres ?
les facteurs principaux \(u^1,\ldots,u^p\) sont les vecteurs propres d’une matrice qu’on d?terminera. Quelles sont les valeurs propres ?
Représentation des variables/indiviuds supplémentaires
Individu supplémentaire
Supposons qu’on dispose d’une observation supplémentaire et qu’on voudrait la représenter sur un \(k-\)ième axe principal. Pour cela il nous suffit de projeter ce point sur l’axe engendré par \(u^k\). Si \(x_s=\left(x^1_s,\ldots,x^p_s\right)'\), la coordonnée de \(x_s\) sur \(<u^k>\) est \[<x_s,u^k>_M=(x_s)'Mu^k\]
Variable qualitative supplémentaire
Supposons qu’on dispose d’une variable qualitative supplémentaire et qu’on voudrait la représenter. Sachant que cette variable divise notre echantillon en sous-groupes, nous calculons alors les centres de gravités de chaque sous-groupe et nous traitons chaque centre de gravité comme une observation supplémentaire.
Variable quantitative supplémentaire
Supposons qu’on dispose d’une variable quantitative supplémentaire et qu’on voudrait la représenter sur un \(k-\)ième axe unitaire engendré par \(<Z^k>\). Pour cela il nous suffit de projeter ce vecteur sur l’axe engendr? par \(Z^k\). Si \(x^s=\left(x_1^s,\ldots,x_n^s\right)'\), la coordonnée de \(x^s\) sur \(<u^k>\) est \[<x^s,Z^k>_D=(x^s)'Dz^k=\mbox{Cor}(x^s,C^k)\times \sqrt{\mbox{Var}(x^s)}\]
Exercice
On considère \(x^{1}=\left(x_{1}^{1},\ldots,x_{n}^{1}\right)'\) et \(x^{2}=\left(x_{1}^{2},\ldots,x_{n}^{2}\right)\) deux séries d’observations de deux variables centr?es \(x^{1}\) et \(x^{2}\). On a calculé la matrice variance-covariance empirique entre ces deux variables : \[ S=\left(\begin{array}{cc} 1 & a\\ a & 1 \end{array}\right) \]
où \(a>0.\) On note par \(\lambda_{1}\) et \(\lambda_{2}\) les deux valeurs propres de \(S\) telles que \(\lambda_{1}\geq\lambda_{2}.\) On considère dans \(\mathbb{R}^{2}\) le produit scalaire usuel.
- Calculer \(\lambda_{1}\) et \(\lambda_{2}\).
- Calculer deux vecteurs propres orthogonaux \(u^{1}\) et \(u^{2}\) de normes égales 1 de la matrice \(S.\)
- Calculer les deux composantes principales \(C^{1}\) et \(C^{2}\) en fonction de \(x^{1}\), \(x^{2}\) et \(a\).
- En déduire la moyenne et la variance de \(C^{1}\) et de \(C^{2}\) en fonction de \(a\).
- On note par \(C\) la matrice de colonnes \(C^{1}\) et \(C^{2}\), par \(X\) la matrice de colonnes \(x^{1}\) et \(x^{2}\) et par \(u\) la matrice orthogonale des vecteurs propres \(u^{1}\) et \(u^{2}\) de \(S\). Montrer que \(C^{1}\) et \(C^{2}\) sont deux variables non correl?es.
Indiquer si les propositions suivantes sont exactes ou non. Justifier :\%
- \(C'C=\left(\begin{array}{cc} \lambda_{1} & 0\\ 0 & \lambda_{2} \end{array}\right)\)
- \(C'C=\left(\begin{array}{cc} n\lambda_{1} & 0\\ 0 & n\lambda_{2} \end{array}\right)\)
- \(S=X'X\)
- \(S=(1/n)XX'\)
- \(S=(1/n)X'X\)
- \(S=u'u\)
- \(C=Xu\)
- \(u\) est une matrice \(2\times2\)
Qualité de la représentation.
Des individis
Contribution d’un individu \(w_i\) à l’inertie totale \[\mbox{Ctr}(x_i)=\displaystyle\frac{p_i \, ^t(x_i) x_i}{ I}=\displaystyle\frac{\sum_{k=1}^d (C_i^k)^2 }{\sum_{k=1}^d \lambda_k}\] \(\mbox{Ctr}(x_i)\) permet d’indiquer la présence d’une observation aberrante.
Contribution d’un individu \(w_i\) à l’inertie du \(k-\)ième axe principal \(<u^k>\) \[\mbox{Ctr}_k(x_i)=\displaystyle\frac{(C_i^k)^2}{\sum_{l=1}^d \lambda_l}\]
\(\cos^2\) de l’angle entre \(x_i\) et \((x_i)^\bot_k\) mesure la qualité de la représentation de \(x_i\) sur \(<u^k>\) \[\cos^2_k(x_i)=\displaystyle\frac{(C_i^k)^2}{\,^t(x_i)x_i}\]
Des variables
- Rappelons que
\[\forall,\, j=1,\ldots,d, \;\;x^j=\sum_{k=1}^d u_j^k \sqrt{\lambda_k} \,Z^k \] \[\forall,\, k=1,\ldots,d, \;\;C^k=\sum_{j=1}^d u_j^k x^j \]
\(\Rightarrow\) \(\mbox{Var}(x^j)=\parallel x^j\parallel^2=\sum_{k=1}^d (u_j^k)^2\lambda_k\)
[\(\Rightarrow\)] \(\mbox{Ctr}_k(x^j\sim Z^k)=\displaystyle\frac{(u_j^k)^2\lambda_k}{\mbox{Var}(x^j)}\) est la contribution de \(Z^k\) dans l’explication de \(x^j\).
- Calculons la corr'elation entre \(x^j\) et \(C^k\) ?
\[\mbox{Cor}(x^j,C^k)=\mbox{Cor}(x^j,Z^k)=\displaystyle\frac{<x^j,Z^k>}{\sqrt{\mbox{Var}(x^j)\,\mbox{Var}(Z^k)}} =\displaystyle\frac{u_j^k\sqrt{\lambda_k}}{\sqrt{\mbox{Var}(x^j)}}.\] Si \(X\) est une matrice centr'ee-r'eduite (i.e., \(\mbox{Var}(x^j)=1,\) \(\forall\, j\). \[\mbox{Cor}(x^j,C^k)=u_j^k\,\sqrt{\lambda_k}.\]
On calcule de \(\cos^2\) de l’angle entre \(x^j\) et \((x^j)_k^\bot\) : \[\cos^2\left(x^j,(x^j)_k^\bot\right)=\displaystyle\frac{\parallel (x^j)_k^\bot \parallel^2}{\parallel x^j \parallel^2}=\displaystyle\frac{(u_j^k)^2\lambda_k}{\mbox{Var}(x^j)}=\left(\mbox{Cor}(x^j,C^k)\right)^2.\]
Dans le cas d’une ACP normée, \(\mbox{Var}(x^j)=\parallel x^j\parallel^2=1,\) \(\forall\, j\).
\(\parallel (x^j)_k^\bot \parallel^2\leq \parallel x^j\parallel^2=1\) \(\Rightarrow\) la projection de \(x^j\) se trouve à l’intérieur d’un cercle de rayon :
\(\parallel (x^j)_k^\bot \parallel^2=\cos^2\left(x^j,(x^j)_k^\bot\right)=\left(\mbox{Cor}(x^j,C^k)\right)^2=(u_j^k)^2\lambda_k\)
Si \((x^j)^\bot\) est la projection de \(x^j\) sur le plan \(<Z^1,Z^2>\), alors \((x^j)^\bot\) se trouve à l’intérieure du cercle de corrélation.
Mise en oeuvre sous R.
Etapes d’une mise en oeuvre d’une ACP
Identifier les variables actives dans une ACP.
Possibilités de diviser l’echantillon en deux sous-echantillons
- Echantillon d’apprentissage
- Echantillon Validation/Test
Construire la matrice active \(X\) :
- Centrer \(X\).
- Faut-il réduire \(X\) ?
- Centrer \(X\).
Choisir l’outil
R:PCAdu packageFactoMineRdudi.pcadu packageade4prcomp()ouprincomp()
Faire les calculs sans aucun package
Utilisation de FactoMineR
1. Installer FactoMineR
> install.packages("FactoMineR")
> install.packages("devtools")
> library("devtools")
> install_github("kassambara/factoextra")2. Charger les packages
> require(FactoMineR)
> require(devtools)
Loading required package: devtools
> require(factoextra)
Loading required package: factoextra
Loading required package: ggplot2
Loading required package: gridLes variables actives : 10 scores du decathlon :
> colnames(decathlon)[1:10]
[1] "100m" "Long.jump" "Shot.put" "High.jump" "400m"
[6] "110m.hurdle" "Discus" "Pole.vault" "Javeline" "1500m" Les variables supplémentaires (Quantitatives) :
> colnames(decathlon)[11:12]
[1] "Rank" "Points"La variable supplémentaire (Qualitative) : Competition
> res.pca=PCA(X = decathlon,scale.unit = T,ncp = 3,quanti.sup = 11:12,
+ quali.sup = 13,graph=F)3. Les objets dans res.pca
> names(res.pca)
[1] "eig" "var" "ind" "svd" "quanti.sup"
[6] "quali.sup" "call" Les valeurs propres sont des variances
> res.pca$eig
eigenvalue percentage of variance
comp 1 3.2719055 32.719055
comp 2 1.7371310 17.371310
comp 3 1.4049167 14.049167
comp 4 1.0568504 10.568504
comp 5 0.6847735 6.847735
comp 6 0.5992687 5.992687
comp 7 0.4512353 4.512353
comp 8 0.3968766 3.968766
comp 9 0.2148149 2.148149
comp 10 0.1822275 1.822275
cumulative percentage of variance
comp 1 32.71906
comp 2 50.09037
comp 3 64.13953
comp 4 74.70804
comp 5 81.55577
comp 6 87.54846
comp 7 92.06081
comp 8 96.02958
comp 9 98.17773
comp 10 100.000004. Eboulis des valeurs propres/Scree plot
Avec facto_extra
> fviz_screeplot(res.pca, ncp = 10) + theme_classic()
Avec rbokeh
Utilisation du package rbokeh pour la repreésentation de l’eébouli des valeurs propres.
> eigs = cbind.data.frame(paste("Comp", 1:nrow(res.pca$eig), sep = " "), res.pca$eig)
> colnames(eigs)[1] = "Comps"
> colnames(eigs)[3] = "Percentage"
> colnames(eigs)[4] = "Cumulative"
> xlim = eigs$Comps
> library(rbokeh)
> figure(title = "Eboulis des valeurs propres", xlim = xlim) %>% ly_bar(Comps,
+ Percentage, data = eigs) %>% ly_lines(Comps, Percentage, data = eigs) %>%
+ ly_points(Comps, Percentage, data = eigs, hover = list(Comps, Percentage,
+ Cumulative)) %>% x_axis(label = "") %>% y_axis(label = "Pourcentage de variance") %>%
+ theme_axis("x", major_label_orientation = 45)5. Représentation des individus
Avec FactoMineR
> plot.PCA(res.pca, choix = "ind")
Avec facto_extra
Avec les noms des indidivus.
> fviz_pca_ind(res.pca, geom = "text") + theme_classic()
En colorant selon le cos2
- On calcul d’abord le
cos2sur le plan 1-2
> cos2 = rowSums(res.pca$ind$cos2[, 1:2])
> median(cos2)
[1] 0.3374411- On peut maintenant changer l’échelle des couleurs selon la médiane du
cos2sur le plan 1-2
> fviz_pca_ind(res.pca, geom = "text", col.ind = "cos2") + scale_color_gradient2(low = "white",
+ mid = "blue", high = "red", midpoint = median(cos2), space = "Lab") + theme_classic()
- Colorer selon le groupe.
> fviz_pca_ind(res.pca, geom = "text", habillage = decathlon$Competition) + theme_classic()
- Ajouter des ellipses de confiance
> fviz_pca_ind(res.pca, geom = "text", habillage = decathlon$Competition, addEllipses = TRUE,
+ ellipse.level = 0.95) + theme_classic()
Avec ggplot2 : Creér ces propres graphiques
- Tracer les individus comme des cercles proportionnels à leurs
cos2
> library(ggplot2)
> dt = cbind.data.frame(res.pca$ind$coord[, 1:2], cos2, decathlon$Competition,
+ rownames(decathlon))
> colnames(dt) = c("PC1", "PC2", "Cos2", "Competition", "Athlets")
> p <- ggplot(data = dt, aes(x = PC1, y = PC2, col = Competition, size = Cos2)) +
+ geom_hline(yintercept = 0, alpha = 0.4) + geom_vline(xintercept = 0, alpha = 0.4) +
+ geom_point(alpha = 0.5) + xlab(paste("Axis 1 (", round(res.pca$eig[1, 2],
+ 1), "%)", sep = "")) + ylab(paste("Axis 2 (", round(res.pca$eig[2, 2], 1),
+ "%)", sep = ""))
> p + theme_classic()
- Mettre les noms des athlètes
> p <- ggplot(data = dt, aes(x = PC1, y = PC2, col = Competition, label = Athlets)) +
+ geom_hline(yintercept = 0, alpha = 0.4) + geom_vline(xintercept = 0, alpha = 0.4) +
+ geom_text() + xlab(paste("Axis 1 (", round(res.pca$eig[1, 2], 1), "%)",
+ sep = "")) + ylab(paste("Axis 2 (", round(res.pca$eig[2, 2], 1), "%)", sep = ""))
> p + theme_classic()
- Eviter les chevauchements dans les noms et utilisation du package
ggrepel
> library(ggrepel)
> p <- ggplot(data = dt, aes(x = PC1, y = PC2, col = Competition, label = Athlets)) +
+ geom_hline(yintercept = 0, alpha = 0.4) + geom_vline(xintercept = 0, alpha = 0.4) +
+ geom_point() + geom_text_repel() + xlab(paste("Axis 1 (", round(res.pca$eig[1,
+ 2], 1), "%)", sep = "")) + ylab(paste("Axis 2 (", round(res.pca$eig[2, 2],
+ 1), "%)", sep = ""))
> p + theme_classic()
- Mettre en évidence les groupes avec des convexes de Hull
> library(ggplot2)
> library(scales)
> library(grid)
> library(plyr)
> library(gridExtra)
> find_hull <- function(X) X[chull(X$PC1, X$PC2), ]
> hulls <- ddply(dt, "Competition", find_hull)
> p <- ggplot(data = dt, aes(x = PC1, y = PC2, col = Competition, fill = Competition)) +
+ geom_hline(yintercept = 0, alpha = 0.4) + geom_vline(xintercept = 0, alpha = 0.4) +
+ geom_polygon(data = hulls, alpha = 0.2) + geom_point() + xlab(paste("Axis 1 (",
+ round(res.pca$eig[1, 2], 1), "%)", sep = "")) + ylab(paste("Axis 2 (", round(res.pca$eig[2,
+ 2], 1), "%)", sep = ""))
> p <- p + theme_classic()
> p
Avec plotly
> library(plotly)
Attaching package: 'plotly'
The following object is masked from 'package:ggplot2':
last_plot
The following object is masked from 'package:graphics':
layout
> ggplotly(p)avec rbokeh
> library(rbokeh)
> v1 = paste("Axis 1 (", round(res.pca$eig[1, 2], 1), "%)", sep = "")
> v2 = paste("Axis 2 (", round(res.pca$eig[2, 2], 1), "%)", sep = "")
> p <- figure(width = 600, height = 600) %>% ly_points(PC1, PC2, data = dt, color = Competition,
+ hover = list(PC1, PC2, Cos2, Athlets)) %>% ly_polygons(PC1, PC2, group = Competition,
+ data = hulls, color = Competition) %>% x_axis(label = v1) %>% y_axis(label = v2)
> p6. Représentation des variables
Avec FactoMineR
> plot.PCA(res.pca, choix = "var")
Avec facto_extra
> fviz_pca_var(res.pca, geom = "text") + theme_classic()
avec des flèches
> fviz_pca_var(res.pca, geom = c("text", "arrow")) + theme_classic()
en colorant selon le cos2
> fviz_pca_var(res.pca, geom = c("text", "arrow"), col.var = "cos2") + theme_classic()
On peut choisir de représenter les variables ayant un cos2 supérieure à une certaine valeur.
> fviz_pca_var(res.pca, geom = c("text", "arrow"), select.var = list(cos2 = 0.5)) +
+ theme_classic()
On peu aussi garder les 3 variables les mieux repreésentées
> fviz_pca_var(res.pca, geom = c("text", "arrow"), select.var = list(cos2 = 3)) +
+ theme_classic()
Avec ggplot2
On construit le cercle
> circle <- function(center = c(0, 0), npoints = 100) {
+ r = 1
+ tt = seq(0, 2 * pi, length = npoints)
+ xx = center[1] + r * cos(tt)
+ yy = center[1] + r * sin(tt)
+ return(data.frame(x = xx, y = yy))
+ }
> corcir = circle(c(0, 0), npoints = 100)On récupère les coordonneées des variables sur les axes principaux
> correlations = as.data.frame(res.pca$var$coord)
> colnames(correlations) = c("PC1", "PC2", "PC3")Les données pour construire les flèches
> d = nrow(correlations)
> arrows = data.frame(x1 = rep(0, d), y1 = rep(0, d), x2 = correlations[, 1],
+ y2 = correlations[, 2])Tracer le cercle des corrélations
> cercor <- ggplot() + geom_path(data = corcir, aes(x = x, y = y), colour = "gray65") +
+ geom_segment(data = arrows, aes(x = x1, y = y1, xend = x2, yend = y2), colour = "gray65") +
+ geom_text_repel(data = correlations, aes(x = PC1, y = PC2, label = rownames(correlations))) +
+ geom_hline(yintercept = 0, colour = "gray65") + geom_vline(xintercept = 0,
+ colour = "gray65") + xlim(-1.1, 1.1) + ylim(-1.1, 1.1) + labs(x = v1, y = v2) +
+ labs(title = "Cercle de correlations")
> cercor + theme_classic()
6. Le biplot des variables et des individus
On repreésente les variables et les individus sur le meême graphe en gardant les labels des variables
> fviz_pca_biplot(res.pca, label = "var") + theme_classic()
Avec un habillage
> fviz_pca_biplot(res.pca, label = "var", habillage = decathlon$Competition, addEllipses = TRUE,
+ ellipse.level = 0.95) + theme_classic()