Analyse en composantes principales
et autres réductions dimensionnelles
en langage R
Parlons d'abord statistiques multivariées.
Lorsque l'on a plusieurs variables, que faire ? Comment simplifier leur approche ? Faisons un bref récapitulatif :
L'analyse en composante principale (ACP) permet de réduire un jeu de données si toutes les données sont numériques. cela permet de les visualiser. Cette réduction n'est pas une obligation, on peut après tout faire de la régression multiple.
L'analyse factorielle des correspondances (AFC) lorsque l'on souhaite comprendre la nature du lien entre deux variables catégorielles.
L'Analyse des Correspondances Multiples (ACM) est une méthode d’analyse factorielle adaptée aux données qualitatives (aussi appelées catégorielles)
L’Analyse Factorielle des Données Mixtes (AFMD), qui est une extension de l’ACP (Analyse en Composantes Principales) pour les données mixtes.
Le tSNE conserve la distance entre les individus sans reproductibilité et sans permettre l'ajout ultérieur d'individus.
L'essentiel de cette page !
L'analyse en composante principale (ACP ou PCA en anglais) permet de réduire le nombre de dimensions d'un jeu de données multidimensionnel.
Cette page explique comment réaliser une ACP avec R mais aussi comment visualiser les résultats.
Elle dirige vers d'autres pages pour la mise en forme.
Voici une démarche simplifiée pour (bien) mener une ACP :
1) Il n'est pas obligatoire de faire l'ACP sur toutes les variables disponibles. Il est bon de tester les corrélations/khi2/effet sur la moyenne AVANT afin de ne retenir QUE les variables qui semblent avoir un effet (effet putatif) sur le paramètre qui nous intéresse.
2) Penser à centrer-réduire les données lors de l'ACP si les variables mesurées ne sont pas de même nature
3) Interprétation du diagramme de pareto : les 2 composantes (si on en souhaite 2) doivent intégrer 80% de la variabilité de l'ensemble du jeu de données.
4) Biplot = graphique classique "la composante 2 en fonction de la composante 1" (y = f(x)) + "CERCLE" . Le cercle avec les flèches permettent grosso modo de voir la participation de chaque variable à chaque composante. 2 flèches de directions opposées semblent avoir des effets opposés. Les composantes sont des sortes de super-variables qui regroupent la variabilité de toutes les variables.
5) Interprétation de la distribution des points : on peut faire 1) du regroupement ou 2) des régressions : l'interprétation, le "comme comprendre ces graph" dépend de l'objectif souhaité.
Si en revanche, l'ACP ne donne pas des résultats satisfaisants, il existe d'autres solutions ! Les plus basiques sont expliquées dans cette page.
Si vous voulez être formé (en Aquitaine) aux autres solutions, plus performantes, n'hésitez-pas à demander un devis !
L'analyse en composantes principales permet, d'après Wikipedia (septembre 2017) de transformer des variables liées entre elles (dites « corrélées » en statistique) en nouvelles variables décorrélées les unes des autres. Ces nouvelles variables sont nommées « composantes principales », ou axes principaux. Elle permet à l'analyste de réduire le nombre de variables et de rendre l'information moins redondante.
Cette page s'intéresse ici non pas à l'ACP en elle-même (présentation succincte rubrique 1) mais à la représentation graphiques des composantes dans le cadre d'une étude ou d'une publication.
1. Réaliser une analyse en composantes principales
L'analyse des données dans le cadre des statistiques multivariées se fait par plusieurs méthodes dont l'ACP (Analyse en Composantes Principales).
Voici comment réaliser une telle étude.
Réaliser une ACP (exemple d'ACP pas à pas avec R par défaut ou les librairies ade4 et FactoMineR)
Étape 1 - Ouvrir des données
Ouverture des données - données cancers.txt
data = read.table("cancers.txt",sep="\t",header=T)
colnames(data)
rownames(data) = data[,1] # On donne à chaque ligne le nom de son cancer (Cancer 1, Cancer 2)
x = data[1:7,-1] ;x # On exclue la colonne 1 qui contenait les noms des cancers.
Dans cet exemple, on ouvre des données simulées
présentant des résultats moyens de concentration de 10 marqueurs sanguins (M1 à M10) pour 7 catégories de cancers (Cancer 1 à 7).
Chaque catégorie de cancer présente ainsi un taux moyen de marqueurs variable d'une catégorie à l'autre.
On cherche à voir s'il existe des cancers qui présenterait une signature sanguine caractéristique.
Pour cela on réalise une ACP.
Étape 2 - Réaliser l'ACP
Plusieurs librairies permettent de réaliser des ACP avec R.
1) La commande par défaut prcomp() (installé sur R par défaut), renvoie une acp dont la description est délicate à obtenir
2) ade4 (dudi.pca()) permet d'obtenir facilement une acp et sa description. On peut facilement voir un diagramme de Pareto pour sélectionner les composantes qui ont du poids
3) factomineR (PCA()) est un package très performant qui permet d'aller bien au-delà de l'ACP.
Méthode par défaut de R (prcomp())
acp1 <- prcomp(x, scale = F)
acp1
Avec la librairie ade4
install.packages("ade4")
library(ade4)
acp2 <- dudi.pca(x , scannf= F,scale=FALSE,,nf=3)
acp2
Avec la librairie FactoMineR
install.packages("FactoMineR")
library(FactoMineR)
acp3 = PCA(x, scale.unit=F, ncp=5, graph=T)
acp3
Paramétrage de la commande prcomp() (R par défaut)
Paramétrage de la commande dudi.pca() (librairie ade4)
Paramétrage de la commande PCA() (librairie FactoMineR)
Étape 3 - Contrôler la pertinence de l'ACP par un diagramme de Pareto
Le principe de Pareto suppose que la variabilité globale se concentre sur les premières composantes.
Afin de pouvoir travailler sur seulement quelques composantes, il faut s'assurer que ce principe de Pareto est respecté.
Il est possible de contrôler la pertinence de l'ACP, non déformation du nuage en s'assurant que les composantes retenues couvrent bien 95%, 97.5% ou plus selon son choix de la variabilité du nuage.
La démarche est le suivante :
1) Il faut identifier où se positionnent les résultats des poids de chaque composante dans la sortie d'ACP.
2) Il faut récupérer une fonction de réalisation d'un diagramme de Pareto
3) Il faut tracer le graphique
Dans cet exemple, comme les 2 premières composantes couvrent 95% de la variabilité, on peut réduire ainsi le nombre de dimensions des données à seulement 2 au lieu de 10.
Pareto de contrôle de l'ACP
Les 2 premières composantes couvrent ici plus de 95% de la variabilité.Méthode par défaut de R
# Installer mon package {KefiR} pour avoir la fonction pareto()
pareto(acp1$sdev)
# Attention : il ne s'agit pas des variances mais des déviations standards : prcomp() ne permet pas une analyse pertinente
Avec la librairie ade4
# Installer mon package {KefiR} pour avoir la fonction pareto()
pareto(acp2$eig, h=95)
# h = 95 : afficher un seuil de cumul de la variance à 95%
Avec la librairie FactoMineR
# Récupérer et copier-coller avant dans la console la fonction pareto()
pareto(acp3$eig[,2], h=95)
# h = 95 : afficher un seuil de cumul de la variance à 95%
Étape 4 - Afficher les données en 2 dimensions - Ex : Axe1 = Composante 1 & Axe2 = Composante2
Une fois que l'on sait sur quelles composantes s'attarder, la première envie est de voir comment vont se répartir les points sur un nombre réduit de dimensions.
Chaque librairie possède ses propres outils de visualisations. Tous sont très limités dans la mise en forme.
Le mieux est donc d'utiliser les commandes de base R comme plot().
Vous trouverez sur ce site une page dédiée aux représentations graphiques avec ade4 (commande s.label()) ou FactoMineR.
On peut aussi se référer aux aides de ces différents packages.
Courbe sur les composantes 1 et 2 (on peut aussi bien représenter d'autres composantes en remplaçant les 1 et 2 par 3 et 4 dans axeX et axeY
# Avec des données ACP issues de prcomp() de R par défaut
axeX <- acp1$x[,1] ; axeY <- acp1$x[,2]
# Avec des données ACP issues de dudi.pca() de ade4
axeX <- acp2$li[,1] ; axeY <- acp2$li[,2]
# Avec des données ACP issues de PCA() de FactoMineR
axeX <- acp3$ind$coord[,1] ; axeY <- acp3$ind$coord[,2]
Tracer le graphique
plot(axeX,axeY,pch=16);grid()
text(axeX,axeY,rownames(x))
Étape 5 - (Optionnelle) - Visualiser la corrélation et l'influence des différentes variables avec les composantes (diagramme circulaire, radial, en radiant)
Un tel diagramme permet de calculer comment chaque variable intervient dans une composante.
Il est utile pour comprendre le lien entre variables et va donner des résultats équivalents à un clustering sur variables.
Dans cet exemple qui met en relation composante 1 (horizontale) et composante 2 (verticale), on voit que les composantes M4 et M5, longues, jouent fortement.
M7 se superpose à la composante 1 sans en être la seule part.
A) Récupérer les données nécessaire à tracer le diagramme
Méthode par défaut de R
# Récupérer la description donnant la corrélation de chaque variable par composante
# 1 et 2 pour composantes 1 et 2
poids_var_x <- acp1$rotation[,1]
poids_var_y <- acp1$rotation[,2]
vecnames <- rownames(acp1$rotation)
Avec la librairie ade4
# Récupérer la description donnant la corrélation de chaque variable par composante
poids_var_x <- acp2$c1[,1]
poids_var_y <- acp2$c1[,2]
vecnames <- rownames(acp2$c1)
Avec la librairie FactoMineR
# Récupérer la description donnant la corrélation de chaque variable par composante
# Pas de contenu prêt à l'emploi mais des fonctions font directement le diagramme
B) Optionnel : filtrer les variables pour n'afficher que celles qui ont un poids suffisant
#OPTIONNEL
poids_var_global <- sqrt(poids_var_x^2+poids_var_y^2)
barplot(poids_var_global,names.arg=vecnames)
# FILTRER : Ex : les variables dont le poids est > 0.25
seuil = 0.25 ; abline(h=seuil,col="red",lwd=2)
poids_var_x <- poids_var_x[poids_var_global>seuil]
poids_var_y <-poids_var_y[poids_var_global>seuil]
vecnames <- vecnames[poids_var_global>seuil]
On voit ici que le seuil de 0.25 va éliminer les variables M2, M6 et M9 de l'affichage.
C) Tracer le diagramme
# Tracer un fond blanc orthonormé
plot(c(-1,1),c(-1,1),type="n",xlab="",ylab="",asp=1,axes=F) ; grid()
# Tracer le cercle
install.packages('plotrix') ; library("plotrix") ;draw.circle(0,0,1,lwd=3)
# Tracer les flèches
arrows(0,0,poids_var_x,poids_var_y,length=0.1,angle=17,lwd=2)
# Calculer la position du nom de chaque vecteur // aux pointes des flèches
poids_x_bin <- poids_var_x ; poids_x_bin[poids_x_bin<0]<- -1 ; poids_x_bin[poids_x_bin>0]<- 1
poids_y_bin <- poids_var_y ; poids_y_bin[poids_y_bin<0]<- -1 ; poids_y_bin[poids_y_bin>0]<- 1
poids_x <- poids_var_x^2 ; poids_y <- poids_var_y^2; poids_z <- poids_x + poids_y
poids_x <- poids_x / poids_z * poids_x_bin * 0.08 ; poids_y <- poids_y / poids_z * poids_y_bin * 0.08
# Tracer les noms des vecteurs
text(poids_var_x+poids_x,poids_var_y+poids_y,vecnames,cex=0.75,col="red")
Étape 6 - Combiner les étapes 4 et 5 pour visualiser la distributions des résultats selon 2 composantes et superposer la projection de chaque variable. (biplot)
On peut obtenir directement ce graphique avec la fonction biplt() de {KefiR}.
Il suffit tout simplement de réaliser le graphique de l'étape 4 puis de faire le graphique de l'étape 5 en séparant les 2 de la ligne suivante :
par(new = T)
# Graphique issu de l'étape 4 - données de prcomp()
axeX <- acp1$x[,1] ; axeY <- acp1$x[,2]
plot(axeX,axeY,pch=16,xlim=c(-100,70),ylim=c(-40,70),cex=3,col=c(1:100),
xlab="Composante 1",ylab="Composante 2");grid()
text(axeX,axeY-5,rownames(x))
# Pour superposition des 2 graphiques
par(new = T)
# Graphique issu de l'étape 5
plot(c(-1,1),c(-1,1),type="n",xlab="",ylab="",asp=1,axes=F)
arrows(0,0,poids_var_x,poids_var_y,length=0.1,angle=17,lwd=2)
draw.circle(0,0,1,lwd=3)
text(poids_var_x+poids_x,poids_var_y+poids_y,vecnames,cex=0.75,col="red")
Le même diagramme, nettement mois esthétique, peut s'obtenir directement avec la fonction biplot()
biplot(prcomp(dataframe))
Étape 7 - Tracer des ellipses, colorer par catégories
Si les catégories ne sont pas encore définies, on peut réaliser un clustering pour identifier ces catégories (knn, k-means, clustering hiérarchique ascendant (CAH)...).
Sinon, on peut les colorer ou les entourer directement avec une ellipse.
Reprenons les étapes 1 à 6 sur une nouvelle étude, complète.
Voici des résultats obtenus par différents étudiants n'ayant pas tous le même bac. On essaye de voir s'il y a un résultats aux examens et type de bac.
# Étape 1 et 2 Code pour faire l'ACP - data_bac.txt
data_bac = read.table(file.choose(),sep="\t",header=T)
colnames(data_bac)
library(ade4)
acp <- dudi.pca(data_bac[,-c(1,2)] , scannf= F,scale=T,center=T,nf=3) ; acp
Biplot maison, mise en forme fine. Tout est paramétrable manuellement.
Exemple d'ACP cf. aide ci-dessus. - données data_bac.txt
# Visualiser les étudiants par bac
plot(acp$li[,1],acp$li[,2],pch=16,col=data_bac$bac,xlab="Composante 1",ylab="Composante 2")
Regrouper des points pour définir des catégories : appliquons un clustering hiérarchique. cf. aide - clustering hiérarchique.
mydi = dist(data_bac[,-c(1,2)])
mytree = cutree((hclust(mydi)), k=5) # 3 catégories ; print(mytree)
Colorer les groupes de points par catégorie - Si les catégories sont déjà définies dans une autre liste : voir aussi ce lien.
# Il suffit de déclarer les catégories dans le paramètres col
plot(acp$li[,1],acp$li[,2], col = mytree,pch=16,cex=1.2,
main="Regroupement des différents types d'étudiant",xlab="Composante 1",ylab="Composante 2")
legend(x="topright", legend=unique(mytree), col=unique(mytree), pch=16,bg="white")
Entourer les points dans une ellipse de confiance - cf. aide - tracer des ellipses
x_ell <- acp$li[,1][mytree==3] ; y_ell <- acp$li[,2][mytree==3] ; xy_ell <- data.frame(x_ell,y_ell)
install.packages("ellipse") ; library("ellipse")
# Ajouter une ellipse (intervalle de confiance à 95%)
lines(ellipse(cov(xy_ell),centre=colMeans(xy_ell),level=0.95),type="l", lty=1, col="green")
Superposer le diagramme de corrélation pour indiquer la description des composantes
# Une fois le graphique tracé
par(new = T) # pour ajouter par-dessus les vecteurs
plot(c(-1,1),c(-1,1),type="n",xlab="",ylab="",asp=1,axes=F) # pour tracer un graphique blanc orthonormé
# optionnel : tracer un cercle
install.packages('plotrix') ; library("plotrix") ;
draw.circle(0,0,1)
# Tracer les vecteurs ==> Remarque acp$c1 contient ces descriptions pour les acp faites avec ade4.
arrows(0,0,acp$c1[,1],acp$c1[,2],length=0.1,angle=17)
# Récupérer les noms des vecteurs
vecnames <- rownames(acp$c1)
# Tracer les noms des vecteurs - Remarque, cf. ETAPE 4 ci-dessus pour bien localiser la légende.
text(acp$c1[,1]+0.07,acp$c1[,2]+0.07,vecnames,cex=0.75)
2. Interagir avec des résultats d'ACP en 2D ou 3D
2.0. Prenons un nouvel exemple plus complexe qui justifierait une meilleure interaction
Il s'agit de reprendre l'exemple d'ACP réalisé Étape 7 ci-dessus.
Voici des résultats obtenus par différents étudiants n'ayant pas tous le même bac. On essaye de voir s'il y a un résultats aux examens et type de bac.
# Code pour faire l'ACP
data_bac = read.table(file.choose(),sep="\t",header=T)
colnames(data_bac)
library(ade4)
acp <- dudi.pca(data_bac[,-c(1,2)] , scannf= F,scale=FALSE,,nf=3) ; acp
Une catégorie particulière d'étudiants est ainsi entourée en vert, celle des bons élèves qui n'auront pas besoin de soutien.
2.1. Visualiser les résultats en 3D ou établir le lien des composantes avec une autre variable
Étape A - Graphique 3D à 3 axes
Illustrer 3 composantes en 3D, entourer d'un ellipsoïde les points correspondant à une catégorie. C'est possible !
cf. Réalisation de graphiques 3D avec R et de gif animés.
# Début du graphique
# Installer la librairie 3D nécessaire à la visualisation
install.packages("rgl") ; library(rgl) # installer puis lancer la librairie
# Tracer le graphique 3D (les bac S)
plot3d(acp$li[,1][data_bac$bac=="S"],acp$li[,2][data_bac$bac=="S"],acp$li[,3][data_bac$bac=="S"],
col="green",type="p",xlab="Dim1",ylab="Dim2",zlab="Dim3")
# Fin du graphique
# Colorer et ajouter d'autres catégories de bacs
points3d(acp$li[,1][data_bac$bac=="STAV"],acp$li[,2][data_bac$bac=="STAV"],acp$li[,3][data_bac$bac=="STAV"],col="red",radius=0.02)
points3d(acp$li[,1][data_bac$bac=="ST2S"],acp$li[,2][data_bac$bac=="ST2S"],acp$li[,3][data_bac$bac=="ST2S"],col="blue",radius=0.02)
spheres3d(acp$li[,1][data_bac$bac=="STL"],acp$li[,2][data_bac$bac=="STL"],acp$li[,3][data_bac$bac=="STL"],col="black",radius=0.02)
# Tracer un ellipsoïde noir pour le bac STL
x = acp$li[,1][data_bac$bac=="STL"] ; y = acp$li[,2][data_bac$bac=="STL"] ; z = acp$li[,3][data_bac$bac=="STL"]
ellipse <- ellipse3d(cov(cbind(x,y,z)), centre=c(mean(x), mean(y), mean(z)), level = 0.95)
plot3d(ellipse, col = "black", alpha = 0.1, add = TRUE, type="shade")
# Tracer un ellipsoïde rouge pour le bac STAV
x = acp$li[,1][data_bac$bac=="STAV"] ; y = acp$li[,2][data_bac$bac=="STAV"] ;z = acp$li[,3][data_bac$bac=="STAV"]
ellipse <- ellipse3d(cov(cbind(x,y,z)), centre=c(mean(x), mean(y), mean(z)), level = 0.95)
plot3d(ellipse, col = "red", alpha = 0.1, add = TRUE, type="shade")
# Réaliser une série de 90 images pour faire un gif animé avec un logiciel comme GIMP, photofiltre, Easy GIF Animator ou gif animator ou Photoscape.
for (i in 1:90) { rgl.viewpoint(i, 20) ;
filename <- paste("pic", formatC(i, digits = 1, flag = "0"), ".png", sep = "")
rgl.snapshot(filename) } ;getwd()
Étape B - Régression 3D couplée au graphique pour mettre en relation une variable y avec les composantes x1 et x2
cf. L'utilisation de la librairie plotly. Régression 3D, nappes et graphiques 3D interactifs : faire appel à plotly !!!
2.2. Construire sa propre interface pour cliquer sur les points pour les identifier, etc.
Développer une interface pour pouvoir afficher la légende d'un point en cliquant dessus.
Il est relativement simple de reprendre le script disponible aux liens ci-dessous pour développer une interface pour lire ses données d'ACP par interaction avec les graphiques.
# à partir d'une acp faite avec ade4
plot(acp$li[,1],acp$li[,2],pch=16, col="#AA1256")
loc_punct(acp$li[,1],acp$li[,2])
La commande loc_punt n'est pas installée sous R, il faut aller récupérer son code :
cf. La commande de localisation d'un point par simple click.
# à partir d'une acp faite avec ade4
plot_id_point(acp$li[,1],acp$li[,2],pch=16, col="#12AA56",col_point="red",labels=c(1:100))
La commande plot_id_point n'est pas installée sous R, il faut aller récupérer son code :
2.3. Utiliser la librairie plotly pour créer automatiquement des graphiques html interactifs
cf. Le lien suivant. plotly permet aussi de rendre les ggplot interactif.
3. Aller plus loin dans la mise en forme des résultats d'ACP...
Simple pour réaliser des ACP, la librairie ade4 ne permet pas de faire de belles mise en forme.
Si factoextra ne permet pas de faire des ACP, elle permet de reprendre les résultats d'autres commandes pour réaliser de magnifiques graphiques facilement mais dont le paramétrage sera forcément limité.
Mises en forme fines avec la librairie FactoMineR
La librairie FactoMineR semble offrir des fonctionnalités équivalentes à factoextra (cf. présentation) avec des intérêts notoires en particulier pour gérer les données manquantes avec missMDA. Une aide hybride entre factoextra et FactoMineR est disponible ici.
Mise en forme des données de façon performante mais complexe avec ggplot2 et plotly
Cette rubrique reste à rédiger. ggplot et plotly ont un peu le propre langage au sein du langage R mais sont très performants pour réaliser tous les types de graphiques automatiquement (croisement de paramètres, régressions, 3D, nappes, ellipses, interactifs...)
La rotation varimax est une sorte de bidouillage qui va permettre de faire tourner les résultats de l'ACP pour que les "composantes" se superposent approximativement à des variables afin de simplifier l'interprétation des résultats.
le but recherché de cette rotation est une facilitation de l'interprétation. on voit ici qu'il sera possible de définir les personnes selon 2 axes (motivation/enthousiasme versus rigueur/connaissances/savoir-faire). (Exemple artificiel)
4. Lorsque l'ACP trouve ses limites : utiliser d'autres méthodes de réduction des dimensions supervisées ou non
Lorsque l'ACP ne donne pas des résultats satisfaisants, ce n'est pas toujours de la faute des résultats. L'ACP consiste à réduire les dimensions en écrasant un nuage de points.
Comme si on écrasait une mouche, si on s'y prend bien, on conserve 2 composantes principales sur la mouche : la longueur et la largeur.
Mais si on l'écrase de travers ? Le résultat d'ACP peut devenir incompréhensible, comme si la tête se retrouvait au niveau des pattes.
Et si on écrase un nuage de points en forme de serpent enroulé ? Idem.
Ou encore, si la variabilité semble plus forte au sein d'une catégorie identifiée qu'entre des catégories ?
Il existe alors d'autres méthodes dont le tSNE semble le plus prometteur, le cmdscale qui repose sur la distance, la projection de Fisher (fonction LDA()) pour créer une distance entre les catégories...
En dernier lieu, si vous n'obtenez toujours pas des résultats satisfaisants, n'hésitez pas à faire appel à nous pour être formés à d'autres outils non disponibles sur ce site ! Demande de devis gratuit.
L'ACM
Sorte d'ACP pour des données non-numériques (qualitatives), l'ACM (Analyse en Composante Multiple) se fait de la façon suivante :
library(ade4)
acm <- dudi.acm(data)
On retrouve un exemple d'ACM réalisée avec {FactoMineR} dans la page dédiée à l'AFC.
cmdscale
Faire une ACP à partir de la matrice de distance entre les points.
Permet de travailler avec des données numériques et qualitatives.
data(iris)
mydi <- dist(iris[,-5])
my_cmd <- cmdscale(mydi)
plot(my_cmd,col=iris[,5])
Faire une ACP qui maximise la différence entre les catégories. Attention : méthode supervisée !
Une ACP exige parfois de la classification/Catégorisation...
Comment faire une réduction dimensionnelle sous python (ACP, tSNE).
