L'ACP, ou Analyse en Composantes Principales, est une méthode d'exploration de données qui consiste à réduire la dimensionnalité du problème pour en extraire l'essentiel. Par une projection dans un espace plus petit, on réduit le nombre de variables, et si on réduit suffisamment on peut en faire un outil de diagnostic graphique. Comme c'est une projection, il est important de comprendre qu'on perd de l'information dans le processus, mais cela permet d'interpréter plus facilement les données.
Table des matières :
La méthode
En pratique, on part d'un tableau avec N lignes représentant les individus mesurés — par exemple 1000 souris — et P colonnes représentant les variables mesurées — par exemple la taille, le poids, l'âge etc. des souris. L'idée, c'est qu'il est facile de représenter sur un graphique une situation à deux variables : la taille sur l'axe des x, le poids sur l'axe des y, un point par souris. On pourrait imaginer ajouter un 3ème axe pour l'âge et visualiser les points dans un "cube" en 3D. Mais à partir de 4 variables, la visualisation devient impossible. Pour réduire la dimension, on va essayer de ne garder que l'essentiel. L'ACP va transformer le tableau pour qu'il ne compte plus que deux (nouvelles) colonnes tout en conservant un maximum d'information.
Sans entrer dans les détails, la transformation implique d'abord une rotation des axes de coordonnées : le premier axe suivra la direction la plus "allongée" (imaginez un poisson ; l'axe va de la tête à la queue), le second une direction perpendiculaire (du ventre vers le dos), le troisième perpendiculaire aux deux premiers (l'épaisseur du poisson), etc. Cela sert à maximiser la variabilité des points le long du premier axe, puis la variabilité du reste le long du second, etc. de sorte qu'il n'y ait presque plus de variation le long des axes suivants. Ces nouveaux axes s'appellent composantes principales ; ce sont P vecteurs de taille P. Mathématiquement, ce sont les vecteurs propres de la matrice de covariance de notre table, et on fait un changement de base.
Ensuite on fait une sélection des composantes : plus on en garde, plus complètement on peut recomposer l'objet (si on n'a gardé que les deux premières composantes du poisson, on ne peut reconstruire qu'un poisson tout plat).
L'important dans tout ca, c'est que si les deux premières composantes cumulent une grande portion de la variabilité totale, on peut ignorer les autres et les opposer dans un graphique.
Un exemple avec R
Ici nous partirons du jeu de données "iris" déjà inclus dans la plupart des distributions de R. La table ressemble à ça :
|
1 2 3 4 5 6 7 8 |
> head(iris) Sepal.Length Sepal.Width Petal.Length Petal.Width Species 1 5.1 3.5 1.4 0.2 setosa 2 4.9 3.0 1.4 0.2 setosa 3 4.7 3.2 1.3 0.2 setosa 4 4.6 3.1 1.5 0.2 setosa 5 5.0 3.6 1.4 0.2 setosa 6 5.4 3.9 1.7 0.4 setosa |
Il y a 150 individus et 4 dimensions/variables. Les 150 sont groupés en 3 espèces setosa, versicolor et virginica. Ce qu'on voudrait savoir, c'est s'il y a une différence importante entre les espèces, et s'il y en a une, quelles sont les variables qui contribuent le plus à expliquer cette différence.
Il y a deux commandes à choix pour faire l'ACP en R (j'expliquerai plus tard pourquoi) : prcomp et princomp. Ici je vais utiliser prcomp ; avec princomp le résultat est presque le même, et l'utilisation est identique aux noms des attributs près.
|
1 |
pca = prcomp(iris[,1 :4]) |
Les composantes principales (PCx, pour "Principal Component") sont les colonnes de pca$rotation :
|
1 2 3 4 5 6 |
> pca$rotation PC1 PC2 PC3 PC4 Sepal.Length 0.36138659 -0.65658877 0.58202985 0.3154872 Sepal.Width -0.08452251 -0.73016143 -0.59791083 -0.3197231 Petal.Length 0.85667061 0.17337266 -0.07623608 -0.4798390 Petal.Width 0.35828920 0.07548102 -0.54583143 0.7536574 |
On voit par exemple que l'axe "PC1" est une combinaison de 0.86 x Petal.Length, 0.36 x de Sepal.Length et Petal.Width, et ‑0.08 x Sepal.Width - c'est-à-dire que la direction de la PC1 est à peu près alignée avec celle de Petal.Length et quasi perpendiculaire à celle de Sepal.Width. On appelle aussi ces nombres "loadings".
pca$sdev donne la fraction d'information (sdev="standard deviation") contenue dans chacune des composantes :
|
1 2 3 4 5 6 7 8 |
> pca$sdev # Les écarts-types "bruts" [1] 2.0562689 0.4926162 0.2796596 0.1543862 > 100 * pca$sdev^2 / sum(pca$sdev^2) # Proportion de variance [1] 92.4618723 5.3066483 1.7102610 0.5212184 # de chaque composante > sum(100 * (pca$sdev^2)[1 :2] / sum(pca$sdev^2)) # Variance totale de [1] 97.76852 # PC1 + PC2, en % |
Ainsi les deux premières composantes à elles seules contiennent 97.77% de l'information, il est donc raisonnable de garder seulement deux variables en ne perdant que 2.33% d'info, ou même une seule (ce que je ferais en pratique, mais ici on en gardera deux pour plus de généralité). On peut aussi rendre compte de ces proportions avec un bar plot des variances :
|
1 |
plot(pca) |
Je ne vous cache pas qu'on peut en produire un plus joli soi-même.
Les données transformées par rotation sont dans pca$x ; chaque colonne est la projection des points sur l'une des composantes principales. Ainsi si on garde seulement les deux premières, on a la projection des points dans le plan PC1-PC2 et on peut en faire un graphique qu'on appelle "biplot" :
|
1 |
plot(pca$x[,1 :2], col=iris[,5]) |
L'axe horizontal est notre PC1, l'axe vertical la PC2. En coloriant les points selon l'espèce de plante, on voit que la PC1 les distingue presque parfaitement :
En fait on a vu que la composante 2 n'est pas très utile, et on pourrait très bien choisir de réduire encore la dimension et séparer nos points selon un seul axe :
|
1 |
plot(pca$x[,1], rep(0,nrow(iris)), col=iris[,5]) |
L'interprétation ensuite est la suivante : d'abord l'espèce setosa (à gauche) est nettement plus reconnaissable que les deux autres entre elles. Ensuite, comme c'est la composante 1 qui sépare les espèces, on regarde de quoi elle est faite. On avait vu qu'elle était composée principalement de Petal.Length ; on en déduit que c'est principalement à la longueur des pétales qu'on reconnaît les espèces. A vrai dire l'interprétation est souvent difficile, et on se contente alors de remarquer que les groupes sont bien différenciés/groupés comme on l'attend.
Il existe aussi une commande R qui fait le biplot toute seule mais… chacun ses goûts. Les flèches qui apparaissent sur le graphique sont les anciens axes de coordonnées (il y a aura donc P flèches — et on ne peut pas en réduire le nombre, ce qui peut devenir très lourd). Il est aussi centré en zéro et réduit (ce n'est pas important pour l'interprétation).
|
1 |
biplot(pca) |
Autres applications
On peut se servir de l'ACP par exemple en RNA-seq pour regarder si l'expression globale des gènes est très différente d'un groupe de patients à un autre, ou si les réplicats de l'expérience sont suffisamment semblables — dans le cas contraire ça permet d'éliminer les mauvais échantillons et rendre l'analyse plus robuste. On considère alors ces fameuses tables avec un gène par ligne et une colonne par échantillon, indiquant dans chaque case l'expression d'un gène dans un échantillon particulier.
De façon similaire, on peut mesurer la présence ou l'absence de chaque variation du génome (SNP) de 3000 Européens, et comparer le graphe de l'ACP à la carte de la région, c'est plutôt spectaculaire. C'est aussi un bon exemple de cas où l'interprétation découle directement de la disposition des points et pas du contenu des composantes principales :
Dans des contextes différents, on peut l'appliquer pour la compression d'images, la reconnaissance de mouvements, etc.
Finalement, c'est aussi un moyen de trouver la "meilleure" droite passant par un nuage de points, un peu comme on le fait en régression linéaire, avec l'avantage que le résultat est le même si on échange les axes (ce n'est pas le cas de la régression !). Cette droite passe par le centre du nuage dans la direction du vecteur PC1.
En bonus, le code pour ce dernier graphique :
|
1 2 3 4 5 6 7 8 9 10 |
library(MASS) x <- mvrnorm(1000, c(0,0), matrix(c(5,1,3,2), 2,2)) plot(x, xlim=c(-8,8), ylim=c(-8,8), xlab="X", ylab="Y") pca = prcomp(x) r = pca$rotation abline(0, r[2,1]/r[1,1], col="red") lin.model = lm(x[,1] ~ x[,2]) abline(0, lin.model$coefficients[2], col="blue") lin.model = lm(x[,2] ~ x[,1]) abline(0, lin.model$coefficients[2], col="blue") |
Prcomp et princomp
La différence entre ces deux commandes R se trouve dans la méthode mathématique utilisée. La décomposition en composantes principales peut se faire soit par SVD (Décomposition en Valeurs Singulières), soit par recherche des vecteurs propres de la matrice de covariance. Ainsi, on remarque que prcomp (SVD) est plus rapide, et même que si on essaye avec la matrice transposée (voir plus bas), princomp va se plaindre car il n'existe pas suffisamment de vecteurs propres :
|
1 2 3 |
> princomp(t(iris[,1 :4])) Error in princomp.default(t(iris[,1 :4])) : 'princomp' can only be used with more units than variables |
alors que pour prcomp tout va très bien :
|
1 |
prcomp(t(iris[,1 :4])) |
De plus, le calcul de la SVD est plus stable numériquement. Les noms des attributs sont aussi différents. Voici la correspondance si on y tient :
prcomp -> princomp
pca$x -> pca$scores
pca$rotation -> pca$loadings[,1:4]
pca$sdev -> pca$sdev
Pour aller plus loin
- On peut aussi toujours s'intéresser au problème inverse en transposant la table des données. Par exemple, si on a une table avec une ligne par gène et une colonne par condition expérimentale, l'ACP telle que nous l'avons présentée ressortira un graphique avec un point par gène, et on espère voir des "clusters" de gènes en fonction des échantillons. Si on transpose la table pour avoir une ligne par condition et une colonne par gène, on obtiendra un graphique avec un point par condition, et on espère voir les réplicats ensemble et les malades bien séparés des contrôles.
- Lorsque les variables sont, comme dans l'exemple "iris", du même type (des longueurs, en cm), tout va bien. Mais si on commence à mélanger des tonnes et des microns, alors les échelles seront différentes et l'ACP sera biaisée en faveur des nombres les plus grands. Pour y remédier, on peut choisir de normaliser les variables avec l'option "scale=TRUE", ce qui divisera chaque colonne par sa variance. Si cette normalisation s'avère trop forte (c'est le cas avec les expressions de gènes), on peut simplement appliquer un logarithme.
- La théorie part du principe que les variables initiales sont indépendantes (axes orthogonaux). En pratique on ignore joyeusement cette contrainte, comme par exemple dans le cas de l'expression des gènes. Un inconvénient, c'est que les composantes principales contiennent un mélange d'un peu toutes les variables à la fois, rendant l'interprétation difficile. Il existe une adaptation de l'ACP pour variables corrélées qu'on appelle "ACP éparse" ("sparse PCA"), exploitant la technique du lasso en interprétant l'ACP comme un problème de régression. Alors un grand nombre de coefficients des composantes principales deviennent nuls et on peut plus facilement extraire les variables significatives. Pour plus de détails, voir H. Zou, T. Hastie, R. Tibshirani, "Sparse principal component analysis", 2004. Les auteurs ont implémenté la méthode dans le package R elasticnet.
- Il existe des packages R comprenant des fonctions pour une étude plus poussée de l'ACP, d'ailleurs développés par des Français : FactomineR, Ade4 and amap. Je ne les ai pas testés, mais je vous invite à regarder par exemple ce lien.
Laisser un commentaire