Représenter rapidement une ACP avec R et ggplot2

Je ne sais pas pour vous, mais moi, à chaque fois que j'assiste à une réunion de labo, il y a qua­si sys­té­ma­ti­que­ment un gra­phique d'ACP pour mon­trer les don­nées. Et à chaque fois, il s'agit d'un gra­phique de base, géné­ré avec R, avec la fonc­tion plot(), des cou­leurs qui piquent les yeux et des axes et légendes illi­sibles. La cri­tique est facile me direz-vous, j'avoue avoir moi aus­si pré­sen­té ce genre de gra­phique assez sou­vent. Mais au bout d'un moment, vu que les ACP ça devient très très rou­ti­nier dans le tra­vail d'un "data ana­lyste", la néces­si­té de pro­duire rapi­de­ment et faci­le­ment un gra­phique ACP se fait sen­tir.

À cette occa­sion, j'ai retrous­sé mes manches et j'ai écrit un petit bout de code mai­son qui génère tout seul les gra­phiques pour un nombre don­né de com­po­santes prin­ci­pales, et ce en une seule fonc­tion. J'en ai pro­fi­té pour uti­li­ser ggplot2, dont je vous invite à lire la pré­sen­ta­tion, pour rendre ces gra­phiques un peu plus esthé­tiques. Je par­tage donc avec vous mon secret pour faire des gra­phiques sexy. Ce n'est pas très com­pli­qué en soit, c'est juste des gra­phiques, mais je vous donne mon code si jamais ça peut aus­si vous ser­vir et vous épar­gner du temps.

Aper­çu du type de gra­phique que pro­duit le code qui va suivre

Première étape : chargement de la librairie

Comme je l'ai dit plus haut, nous allons uti­li­ser la librai­rie ggplot2. Cette fonc­tion per­met d'installer auto­ma­ti­que­ment le paquet s'il n'est pas déjà ins­tal­lé.

Deuxième étape : le screeplot

Lorsque l'on fait une ACP, on aime bien savoir quelles com­po­santes prin­ci­pales (PC) accu­mulent le plus de variance (pour les détails sur ce qu'est une ACP, c'est par là). Pour aller vite, on peut sim­ple­ment lan­cer la com­mande scree­plot(), mais le ren­du n'est pas hyper sexy. En plus il n'y a pas de légende sur l'axe des abs­cisses, ce qui est fâcheux quand on a plus de 20 PC. D'ailleurs, quand on lance une ACP sur un set de don­nées avec plus de 100 échan­tillons (si si, ça arrive), on se retrouve avec un bar­plot tout ser­ré où en réa­li­té l'information qui nous inté­resse est com­pres­sée sur la gauche du gra­phique.

Pour résoudre cela, la fonc­tion sui­vante prend en entrée le résul­tat de l'ACP (cal­cu­lée avec prcomp), et le nombre de PC à affi­cher. C'est pas le plus "hot" des gra­phiques mais ça peut être pra­tique.

Troisième étape : l'ACP

Pas­sons au nerf de la guerre : l'ACP en elle même. Pour appré­cier la com­plexi­té d'un set de don­nées, il est sou­vent néces­saire de regar­der un peu plus loin que juste les 2 pre­mières com­po­santes. Ce qui me fatigue le plus (il ne me faut pas grand chose), c'est de relan­cer la com­mande de gra­phique pour chaque com­bi­nai­son. Pour satis­faire ma fai­néan­tise, il me faut donc une fonc­tion avec le nombre de PC vou­lues et hop, ça me sort toutes les com­bi­nai­sons pos­sibles.

La fonc­tion sui­vante prend en entrée le résul­tat d'une ACP (cal­cu­lée avec prcomp), le nombre de PC à regar­der, les condi­tions des échan­tillons (une liste qui fait la même taille que le nombre d'échantillons), et une palette de cou­leur en hexa­dé­ci­male (par exemple : "#fb7072"), avec autant de cou­leurs que de condi­tions dif­fé­rentes (e.g. deux condi­tions, cas et contrôle, donc 2 cou­leurs).

Et enfin : la démo

Et voilà le résultat :

Screeplot et ACP
Résul­tat des com­mandes pré­cé­dentes. À gauche le scree­plot, à droite l'ACP avec les deux pre­mières com­po­santes et le pour­cen­tage de variance entre paren­thèse, le tout colo­ré en fonc­tion des condi­tions (1 à 3).

 

Juste une der­nière remarque : si vous spé­ci­fiez plus de 2 PC à affi­cher, ça va géné­rer les gra­phiques cha­cun leur tour, donc je vous conseille d'appeler la fonc­tion plot_​pca() dans un pdf pour les sau­ve­gar­der dans un seul fichier. Peut-être un jour j'essayerai le paquet gri­dEx­tra pour affi­cher plu­sieurs gra­phiques sur une même page… Un jour…

 

Mer­ci à Jnsll, Nico M., et Mat Blum pour leurs com­men­taires et relec­tures.



Pour continuer la lecture :


Commentaires

7 réponses à “Représenter rapidement une ACP avec R et ggplot2”

  1. Super article, mer­ci 🙂

    Pour les figures mul­ti-pan­neaux, j'aime bien cow­plot, que je trouve plus facile d'utilisation que juste gri­dEx­tra :
    https://cran.r‑project.org/web/packages/cowplot/vignettes/introduction.html

    Sinon, le package ggfor­ti­fy per­met de plo­ter les objet issus de prcomp() direc­te­ment :
    https://cran.r‑project.org/web/packages/ggfortify/vignettes/plot_pca.html
    Il y a quelques fea­tures en plus, comme plo­ter la pro­jec­tion des axes d'origines sur le nou­veau plan, et des fea­tures en moins, comme l'absence du pour­cen­tage de variances sur les axes.

  2. Bon­jour et mer­ci pour ton article.

    Pour que le code fonc­tionne chez moi j'ai été obli­gé de modi­fier les deux lignes sui­vantes :

    # pour le bar­plot :
    geom_bar(stat="identity") au lieu de geom_​col()

    # pour la fonc­tion apply():
    scale_colour_manual(values = col, name = "") au lieu de scale_fill_manual(values=colours, name="")

    1. Mer­ci pour ce com­men­taire.
      J'ai en effet oublié de pré­ci­ser que j'utilise ggplot2 ver­sion 2.2.1. Il peut donc y avoir de petites sub­ti­li­tés comme celle-ci si vous uti­li­sez une ver­sion dif­fé­rente du paquet.

  3. Mer­ci ! Pour com­plé­ter : l'excellent package explor (https://cran.r‑project.org/web/packages/explor/vignettes/introduction_fr.html) per­met d'"explorer" de façon inter­ac­tive (inter­face shi­ny) les resul­tats d'acps (ou acms) pro­ve­nant de dif­fe­rents packages R. Indis­pen­sable !

  4. Avatar de Rachel

    Bon­jour ! J'ai uti­li­sé votre script pour mes don­nées qui consistent en un data­frame de 9 variables numé­riques qui peuvent prendre des valeurs nulles. Voi­ci ce que je rem­place au script pour com­po­ser 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 fonc­tion du gra­phique ACP "plot_​pca" voi­ci ce que R me retourne :
    Error : Aes­the­tics must be either length 1 or the same as the data (61): fill

    J'ai éplu­ché un peu les forums pour modi­fier le script en consé­quence mais je n'ai pas réus­si à me défaire de cette erreur. Je pense que je dois caler un length=length(db) quelque part…

    Sau­riez vous m'aider s'il vous plait ?

    1. Bon­jour Rachel,

      ce script est fait pour prendre le résul­tat d'une ACP cal­cu­lée avec prcomp. HCPC a un for­mat de sor­tie dif­fé­rent de prcomp, il va donc fal­loir adap­ter le code pour rem­pla­cer toutes les allu­sions au for­mat de sor­ti de prcomp (tous les "pca$…") par leurs équi­va­lents HCPC.

      Si tu n'y arrives pas, je m'y col­le­rai de mon côté pour t'aider.

      1. Avatar de Rachel

        Ah d'accord ! Mer­ci beau­coup, je vais ten­ter et dans ce cas si j'ai encore un sou­ci je revien­drais vers vous !

Laisser un commentaire