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 soi, ce sont juste des graphiques, mais je vous donne mon code si jamais ça peut aussi vous servir et vous épargner du temps.
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. Ce n'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 :
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.
Laisser un commentaire