Je ne sais pas pour vous, mais moi, à chaque fois que j'assiste à une réunion de labo, il y a quasi systématiquement un graphique d'ACP pour montrer les données. Et à chaque fois, il s'agit d'un graphique de base, généré avec R, avec la fonction plot(), des couleurs qui piquent les yeux et des axes et légendes illisibles. La critique est facile me direz-vous, j'avoue avoir moi aussi présenté ce genre de graphique assez souvent. Mais au bout d'un moment, vu que les ACP ça devient très très routinier dans le travail d'un "data analyste", la nécessité de produire rapidement et facilement un graphique ACP se fait sentir.
À cette occasion, j'ai retroussé mes manches et j'ai écrit un petit bout de code maison qui génère tout seul les graphiques pour un nombre donné de composantes principales, et ce en une seule fonction. J'en ai profité pour utiliser ggplot2, dont je vous invite à lire la présentation, pour rendre ces graphiques un peu plus esthétiques. Je partage donc avec vous mon secret pour faire des graphiques sexy. Ce n'est pas très compliqué en soit, c'est juste des graphiques, mais je vous donne mon code si jamais ça peut aussi vous servir et vous épargner du temps.
Aperçu du type de graphique que produit le code qui va suivre
Première étape : chargement de la librairie
Comme je l'ai dit plus haut, nous allons utiliser la librairie ggplot2. Cette fonction permet d'installer automatiquement le paquet s'il n'est pas déjà installé.
|
1 |
if (!require("ggplot2")) install.packages("ggplot2") |
Deuxième étape : le screeplot
Lorsque l'on fait une ACP, on aime bien savoir quelles composantes principales (PC) accumulent le plus de variance (pour les détails sur ce qu'est une ACP, c'est par là). Pour aller vite, on peut simplement lancer la commande screeplot(), mais le rendu n'est pas hyper sexy. En plus il n'y a pas de légende sur l'axe des abscisses, ce qui est fâcheux quand on a plus de 20 PC. D'ailleurs, quand on lance une ACP sur un set de données avec plus de 100 échantillons (si si, ça arrive), on se retrouve avec un barplot tout serré où en réalité l'information qui nous intéresse est compressée sur la gauche du graphique.
Pour résoudre cela, la fonction suivante prend en entrée le résultat de l'ACP (calculée avec prcomp), et le nombre de PC à afficher. C'est pas le plus "hot" des graphiques mais ça peut être pratique.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 |
plot_percent_var <- function(pca, pc){ # Calcule du pourcentage de variance percent_var_explained <- (pca$sdev^2 / sum(pca$sdev^2))*100 # Préparation d'un tableau avec le numéro des composantes principales # et le pourcentage de variance qui lui est associé percent_var_explained <- data.frame( PC=1:length(percent_var_explained), percent_Var=percent_var_explained ) # Récupérer uniquement le nombre de PC indiqué en argument sub_percent_var_explained <- percent_var_explained[1:pc,] # Génère le graphique p <- ggplot(sub_percent_var_explained, aes(x=PC, y=percent_Var)) + # Génère un barplot geom_col()+ # Utilise le thème "black and white" theme_bw() + # Renomme l'axe des abscisses xlab("PCs") + # Renomme l'axe des ordonnées ylab("% Variance") + # Titre du graphique ggtitle("Screeplot")+ # Option de taille des éléments textuels theme( axis.text=element_text(size=16), axis.title=element_text(size=16), legend.text = element_text(size =16), legend.title = element_text(size =16 ,face="bold"), plot.title = element_text(size=18, face="bold", hjust = 0.5), # Astuce pour garder un graphique carré aspect.ratio=1 ) # Affiche le graphique print(p) } |
Troisième étape : l'ACP
Passons au nerf de la guerre : l'ACP en elle même. Pour apprécier la complexité d'un set de données, il est souvent nécessaire de regarder un peu plus loin que juste les 2 premières composantes. Ce qui me fatigue le plus (il ne me faut pas grand chose), c'est de relancer la commande de graphique pour chaque combinaison. Pour satisfaire ma fainéantise, il me faut donc une fonction avec le nombre de PC voulues et hop, ça me sort toutes les combinaisons possibles.
La fonction suivante prend en entrée le résultat d'une ACP (calculée avec prcomp), le nombre de PC à regarder, les conditions des échantillons (une liste qui fait la même taille que le nombre d'échantillons), et une palette de couleur en hexadécimale (par exemple: "#fb7072"), avec autant de couleurs que de conditions différentes (e.g. deux conditions, cas et contrôle, donc 2 couleurs).
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 |
plot_pca <- function(pca=pca, pc=pc, conditions=conditions, colours=colours){ # Transforme le nombre de PC en argument en nom de PC PCs <- paste("PC",1:pc, sep="") # Calcule le pourcentage de variance par PC percent_var_explained <- (pca$sdev^2 / sum(pca$sdev^2))*100 # Transforme le vecteur de conditions en un facteur cond <- factor(conditions) # Crée un autre facteur avec les conditions col <- factor(conditions) # Change les niveaux du facteur avec la palette de couleur pour attribuer # à chaque condition une couleur levels(col) <- colours # Re-transforme le facteur en vecteur col <- as.vector(col) # Récupère les scores pour le graphique scores <- as.data.frame(pca$x) # Génère toutes les combinaisons possibles de PC PCs.combinations <- combn(PCs,2) # Génère un graphique pour chaque combinaison # avec une boucle apply g <- apply( PCs.combinations, 2, function(combination) { p1 <- ggplot(scores, aes_string(x=combination[1], y=combination[2])) + # Dessine des points avec une bordure de 0.5 remplis avec une couleur geom_point(shape = 21, size = 2.5, stroke=0.5, aes(fill=cond)) + # Utilise le thème "black and white" theme_bw() + # Spécifie la palette de couleur et donne un titre vide à la légende scale_fill_manual( values=colours, name="" ) + # Renomme le titre des axes des abscisses et des ordonnées en "PCx (pourcentage de variance)" avec 3 chiffres après la virgule xlab(paste(combination[1], " (",round(percent_var_explained[as.numeric(gsub("PC", "", combination[1]))], digit=3),"%)", sep=""))+ ylab(paste(combination[2], " (",round(percent_var_explained[as.numeric(gsub("PC", "", combination[2]))], digit=3),"%)", sep=""))+ # Titre du graphique ggtitle("PCA")+ # Option de taille des éléments texte theme( axis.text=element_text(size=16), axis.title=element_text(size=16), legend.text = element_text(size =16), legend.title = element_text(size =16 ,face="bold"), plot.title = element_text(size=18, face="bold", hjust = 0.5), # Astuce pour garder un graphique carré aspect.ratio=1 ) # Affiche le graphique print(p1) } ) } |
Et enfin : la démo
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 |
# Génération d'un set de données aléatoires avec 3 groupes set.seed(12345) x <- c(rnorm(200, mean = -1), rnorm(200, mean = 1.5), rnorm(200, mean = 0.8)) y <- c(rnorm(200, mean = 1), rnorm(200, mean = 1.7), rnorm(200, mean = -0.8)) z <- c(rnorm(200, mean = 0.5), rnorm(200, mean = 7), rnorm(200, mean = 0)) data <- data.frame(x, y, z) # Définition des groupes group <- as.factor(rep(c(1,2,3), each=200)) # Définition de la palette de couleur (on peut aussi utiliser RColorBrewer ou tout autre palette déjà faite) palette <- c("#77b0f3", "#8dcf38", "#fb7072") # On lance le calcule de l'ACP pca <- prcomp(data, center=TRUE, scale=TRUE) # On affiche le graphique "Screeplot" (pourcentage de variance par composante principale) plot_percent_var(pca, 3) # On génère le graphique de l'ACP pour les 2 premières composantes principales plot_pca( pca=pca, pc=2, conditions=group, colours=palette ) |
Et voilà le résultat :
Résultat des commandes précédentes. À gauche le screeplot, à droite l'ACP avec les deux premières composantes et le pourcentage de variance entre parenthèse, le tout coloré en fonction des conditions (1 à 3).
Juste une dernière remarque : si vous spécifiez plus de 2 PC à afficher, ça va générer les graphiques chacun leur tour, donc je vous conseille d'appeler la fonction plot_pca() dans un pdf pour les sauvegarder dans un seul fichier. Peut-être un jour j'essayerai le paquet gridExtra pour afficher plusieurs graphiques sur une même page... Un jour...
Merci à Jnsll, Nico M., et Mat Blum pour leurs commentaires et relectures.
Guillaume Devailly
mai 12, 2017 à 11:35
Super article, merci 🙂
Pour les figures multi-panneaux, j'aime bien cowplot, que je trouve plus facile d'utilisation que juste gridExtra:
https://cran.r-project.org/web/packages/cowplot/vignettes/introduction.html
Sinon, le package ggfortify permet de ploter les objet issus de prcomp() directement:
https://cran.r-project.org/web/packages/ggfortify/vignettes/plot_pca.html
Il y a quelques features en plus, comme ploter la projection des axes d'origines sur le nouveau plan, et des features en moins, comme l'absence du pourcentage de variances sur les axes.
Bigey
mai 17, 2017 à 5:42
Bonjour et merci pour ton article.
Pour que le code fonctionne chez moi j'ai été obligé de modifier les deux lignes suivantes:
# pour le barplot:
geom_bar(stat="identity") au lieu de geom_col()
# pour la fonction apply():
scale_colour_manual(values = col, name = "") au lieu de scale_fill_manual(values=colours, name="")
ZaZo0o
mai 18, 2017 à 1:19
Merci pour ce commentaire.
J'ai en effet oublié de préciser que j'utilise ggplot2 version 2.2.1. Il peut donc y avoir de petites subtilités comme celle-ci si vous utilisez une version différente du paquet.
Ludo Cottret
août 30, 2018 à 2:58
Merci ! Pour compléter : l'excellent package explor (https://cran.r-project.org/web/packages/explor/vignettes/introduction_fr.html) permet d'"explorer" de façon interactive (interface shiny) les resultats d'acps (ou acms) provenant de differents packages R. Indispensable !
Rachel
mars 3, 2019 à 2:23
Bonjour ! J'ai utilisé votre script pour mes données qui consistent en un dataframe de 9 variables numériques qui peuvent prendre des valeurs nulles. Voici ce que je remplace au script pour composer mes groupes :
library(FactoMineR)
hcpc_db<-HCPC(db)
hcpc_db$data.clust
db_clust<-as.data.frame(hcpc_db$data.clust[,"clust"])
db.clust<-cbind(db,db_clust)
group<-as.factor(db.clust$clust)
et lors de l'application de la fonction du graphique ACP "plot_pca" voici ce que R me retourne :
Error: Aesthetics must be either length 1 or the same as the data (61): fill
J'ai épluché un peu les forums pour modifier le script en conséquence mais je n'ai pas réussi à me défaire de cette erreur. Je pense que je dois caler un length=length(db) quelque part...
Sauriez vous m'aider s'il vous plait ?
ZaZo0o
mars 4, 2019 à 10:11
Bonjour Rachel,
ce script est fait pour prendre le résultat d'une ACP calculée avec prcomp. HCPC a un format de sortie différent de prcomp, il va donc falloir adapter le code pour remplacer toutes les allusions au format de sorti de prcomp (tous les "pca$...") par leurs équivalents HCPC.
Si tu n'y arrives pas, je m'y collerai de mon côté pour t'aider.
Rachel
mars 4, 2019 à 11:36
Ah d'accord ! Merci beaucoup, je vais tenter et dans ce cas si j'ai encore un souci je reviendrais vers vous !