Analyse en Composantes Principales (ACP)

Données iris

Nous utilisons pour exemple les données iris, déjà présentes dans R. Elle regroupe les mesures (longueur et largeur des sépales et des pétales) de 150 iris, de trois espèces différentes (setosa, versicolor et virginica). Les voici présentées ci-dessous

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(DT)
datatable(iris)

Pour réaliser l’ACP, nous utilisons le package FactoMineR, très complet et bien documenté. Il existe d’autres versions, mais moins intéressantes à utiliser.

library(FactoMineR)

Application de la méthode sur les données

Pour mettre en oeuvre une ACP, il suffit d’utiliser la fonction PCA() (pour Principal Components Analysis), qui renvoie un objet complxe, comme nous le verrons plus tard. Nous devons bien évidemment supprimer la variable Species. L’opion graph = FALSE évite la création de deux graphiques (cercle des corrélations et premier plan factoriel, que nous créérons par la suite).

acp = PCA(iris[,-5], graph = FALSE)

L’objet complexe (liste nommée) ainsi créé contient toutes les informations nécessaires. Son affichage permet de lister l’ensemble des sous-objets contenus, ainsi qu’une brève description.

acp
## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 150 individuals, described by 4 variables
## *The results are available in the following objects:
## 
##    name               description                          
## 1  "$eig"             "eigenvalues"                        
## 2  "$var"             "results for the variables"          
## 3  "$var$coord"       "coord. for the variables"           
## 4  "$var$cor"         "correlations variables - dimensions"
## 5  "$var$cos2"        "cos2 for the variables"             
## 6  "$var$contrib"     "contributions of the variables"     
## 7  "$ind"             "results for the individuals"        
## 8  "$ind$coord"       "coord. for the individuals"         
## 9  "$ind$cos2"        "cos2 for the individuals"           
## 10 "$ind$contrib"     "contributions of the individuals"   
## 11 "$call"            "summary statistics"                 
## 12 "$call$centre"     "mean of the variables"              
## 13 "$call$ecart.type" "standard error of the variables"    
## 14 "$call$row.w"      "weights for the individuals"        
## 15 "$call$col.w"      "weights for the variables"

Pour avoir un certain nombre d’informations en une fois, on peut afficher le résumé de cet objet, comme ci-dessous.

summary(acp)
## 
## Call:
## PCA(X = iris[, -5], graph = FALSE) 
## 
## 
## Eigenvalues
##                        Dim.1   Dim.2   Dim.3   Dim.4
## Variance               2.918   0.914   0.147   0.021
## % of var.             72.962  22.851   3.669   0.518
## Cumulative % of var.  72.962  95.813  99.482 100.000
## 
## Individuals (the 10 first)
##                  Dist    Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3
## 1            |  2.319 | -2.265  1.172  0.954 |  0.480  0.168  0.043 | -0.128
## 2            |  2.202 | -2.081  0.989  0.893 | -0.674  0.331  0.094 | -0.235
## 3            |  2.389 | -2.364  1.277  0.979 | -0.342  0.085  0.020 |  0.044
## 4            |  2.378 | -2.299  1.208  0.935 | -0.597  0.260  0.063 |  0.091
## 5            |  2.476 | -2.390  1.305  0.932 |  0.647  0.305  0.068 |  0.016
## 6            |  2.555 | -2.076  0.984  0.660 |  1.489  1.617  0.340 |  0.027
## 7            |  2.468 | -2.444  1.364  0.981 |  0.048  0.002  0.000 |  0.335
## 8            |  2.246 | -2.233  1.139  0.988 |  0.223  0.036  0.010 | -0.089
## 9            |  2.592 | -2.335  1.245  0.812 | -1.115  0.907  0.185 |  0.145
## 10           |  2.249 | -2.184  1.090  0.943 | -0.469  0.160  0.043 | -0.254
##                 ctr   cos2  
## 1             0.074  0.003 |
## 2             0.250  0.011 |
## 3             0.009  0.000 |
## 4             0.038  0.001 |
## 5             0.001  0.000 |
## 6             0.003  0.000 |
## 7             0.511  0.018 |
## 8             0.036  0.002 |
## 9             0.096  0.003 |
## 10            0.293  0.013 |
## 
## Variables
##                 Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3    ctr
## Sepal.Length |  0.890 27.151  0.792 |  0.361 14.244  0.130 | -0.276 51.778
## Sepal.Width  | -0.460  7.255  0.212 |  0.883 85.247  0.779 |  0.094  5.972
## Petal.Length |  0.992 33.688  0.983 |  0.023  0.060  0.001 |  0.054  2.020
## Petal.Width  |  0.965 31.906  0.931 |  0.064  0.448  0.004 |  0.243 40.230
##                cos2  
## Sepal.Length  0.076 |
## Sepal.Width   0.009 |
## Petal.Length  0.003 |
## Petal.Width   0.059 |

Focus sur la part de variance expliquée

Le sous-objet eig contient donc les valeurs propres (première colonne) mais aussi (et surtout) les pourcentages de variance expliquée (simple et cumulative) pour chaque composante.

Dans notre cas, on remarque que le premier axe explique à lui seul 73% de l’information. Les deux premières composantes expliquent à elles deux presque 96% de l’information.

acp$eig
##        eigenvalue percentage of variance cumulative percentage of variance
## comp 1 2.91849782             72.9624454                          72.96245
## comp 2 0.91403047             22.8507618                          95.81321
## comp 3 0.14675688              3.6689219                          99.48213
## comp 4 0.02071484              0.5178709                         100.00000

Comme un graphique est parfois plus agréable à lire qu’un tableau, on peut représenter cette information via un diagramme en barres. Pour cela, nous devons modifier légèrement l’objet pour le rendre compatible avec ggplot(), ce que fait le code ci-après.

eig = tibble(data.frame(acp$eig) %>% rownames_to_column("comp"))
eig
## # A tibble: 4 × 4
##   comp   eigenvalue percentage.of.variance cumulative.percentage.of.variance
##   <chr>       <dbl>                  <dbl>                             <dbl>
## 1 comp 1     2.92                   73.0                                73.0
## 2 comp 2     0.914                  22.9                                95.8
## 3 comp 3     0.147                   3.67                               99.5
## 4 comp 4     0.0207                  0.518                             100

On voit nettement la part importante de la première composante.

Le seuil de 25% est utilisé, car si chaque composante fournissait la même quantité d’information, chacune aporterait 25% de celle-ci (100 divisé par 4 composantes, car 4 variables).

ggplot(eig, aes(comp, percentage.of.variance)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  geom_hline(yintercept = 25, linetype = "dashed", color = "gray") +
  annotate("text", x = 4.5, y = 27, label = "25 %") +
  labs(x = "Composantes", y = "% de variance expliquée") +
  theme_minimal()

Représentation des individus

La représentation directe des individus via la fonction plot() n’est pas optimal, mais permet de rapidement voir ce qu’il se passe.

On remarque ici deux groupes d’iris.

Notez qu’il est d’usage de représenter sur le plan factoriel le pourcentage de variance expliquée.

plot(acp, choix = "ind")

Les coordonnées obtenues sont contenues dans le sous-dataframe coord du sous-objet ind. On va s’en servir pour réaliser un graphique plus agréable à lire.

datatable(round(acp$ind$coord, 2))

En ajoutant l’information sur les espèces (contenue donc dans la variable Species du dataframe original), on peut utiliser un nuage de points, avec une couleur en fonction de l’espèce justement.

On remarque que la méthode sépare très clairement les setosa des deux autres espèces (versicolor et virginica), qui semble être légèrement séparées, mais un peu recouvrantes tout de même.

ind = tibble(data.frame(acp$ind$coord)) %>% mutate(Species = iris$Species)
ggplot(ind, aes(Dim.1, Dim.2, color = Species)) +
  geom_point() +
  theme_minimal()

Représentation des variables

Le cercle des corrélations est un outil précieux pour caractériser les axes, afin de mieux comprendre ce qui différencient les éventuels groupes d’individus.

On voit ici que trois variables sont très proches (longueur de sépale ainsi que largeur et longueur de pétale).

plot(acp, choix = "var")

La fonction dimdesc() permet de tester les liens entre chaque variable et chaque composante, avec un test de correlation (seule la p-value est affichée).

Les 4 variables sont liées au premier axe, seul deux sont liées au deuxième et au troisième (pas les mêmes).

dimdesc(acp)
## $Dim.1
## 
## Link between the variable and the continuous variables (R-square)
## =================================================================================
##              correlation       p.value
## Petal.Length   0.9915552 3.369916e-133
## Petal.Width    0.9649790  6.609632e-88
## Sepal.Length   0.8901688  2.190813e-52
## Sepal.Width   -0.4601427  3.139724e-09
## 
## $Dim.2
## 
## Link between the variable and the continuous variables (R-square)
## =================================================================================
##              correlation      p.value
## Sepal.Width    0.8827163 2.123801e-50
## Sepal.Length   0.3608299 5.731933e-06
## 
## $Dim.3
## 
## Link between the variable and the continuous variables (R-square)
## =================================================================================
##              correlation      p.value
## Petal.Width    0.2429827 0.0027349555
## Sepal.Length  -0.2756577 0.0006395628

Caractérisation des axes

A l’aide des informations ci-dessus, il est intéressant de décrire ce que chaque axe explique.

Ici, nous ferions les conclusions suivantes :

  • Le premier axe oppose, à droite, les grands iris (grands pétales et longs sépales, mais avec une faible largeur de sépale) aux petits iris à gauche (petits pétales et cours sépales, mais sépales très larges) ;
  • Le deuxième axe oppose, en haut, les iris avec des grands sépales aux iris avec petits sépales (en bas) ;
  • Le troisième axe expliquant très peu de variance n’est pas à prendre en considération.

Analyse Factorielle des Correspondances (AFC)

Données USAccDeaths

Nous allons travailler sur des données accidents aux Etats-Unis sur la période 1973-1978, déjà présentes elles-aussi dans R. C’est le nombre de mort par mois, sur la période. Elles sont dans un format spécifique, qu’il faut modifier avant de pouvoir les utiliser avec FactoMineR.

Acc = data.frame(matrix(USAccDeaths, ncol = 12, byrow = TRUE, dimnames = list(1973:1978, month.name)))
datatable(Acc, options = list(scrollX = TRUE))

Application de la méthode sur les données

Pour appliquer l’AFC sur le tableau de données, nous utilisons ici la fonction CA() (pour Correspondance Analysis). Ici aussi, nous choississons de ne pas afficher les graphiques au début.

afc = CA(Acc, graph = FALSE)

Comme précédemment, l’objet ainsi créé contient beaucoup d’informations, que nous allons utiliser.

afc
## **Results of the Correspondence Analysis (CA)**
## The row variable has  6  categories; the column variable has 12 categories
## The chi square of independence between the two variables is equal to 442.9585 (p-value =  5.010401e-62 ).
## *The results are available in the following objects:
## 
##    name              description                   
## 1  "$eig"            "eigenvalues"                 
## 2  "$col"            "results for the columns"     
## 3  "$col$coord"      "coord. for the columns"      
## 4  "$col$cos2"       "cos2 for the columns"        
## 5  "$col$contrib"    "contributions of the columns"
## 6  "$row"            "results for the rows"        
## 7  "$row$coord"      "coord. for the rows"         
## 8  "$row$cos2"       "cos2 for the rows"           
## 9  "$row$contrib"    "contributions of the rows"   
## 10 "$call"           "summary called parameters"   
## 11 "$call$marge.col" "weights of the columns"      
## 12 "$call$marge.row" "weights of the rows"

Comme pour PCA(), le résumé permet de voir certains éléments.

summary(afc)
## 
## Call:
## CA(X = Acc, graph = FALSE) 
## 
## The chi square of independence between the two variables is equal to 442.9585 (p-value =  5.010401e-62 ).
## 
## Eigenvalues
##                        Dim.1   Dim.2   Dim.3   Dim.4   Dim.5
## Variance               0.000   0.000   0.000   0.000   0.000
## % of var.             53.754  26.524  10.826   5.384   3.511
## Cumulative % of var.  53.754  80.278  91.105  96.489 100.000
## 
## Rows
##             Iner*1000    Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3
## 1973      |     0.103 |  0.017 14.738  0.538 | -0.013 17.228  0.310 |  0.004
## 1974      |     0.094 | -0.014  8.375  0.335 | -0.014 18.584  0.366 |  0.011
## 1975      |     0.170 |  0.030 39.023  0.866 | -0.001  0.053  0.001 | -0.010
## 1976      |     0.117 |  0.003  0.507  0.016 |  0.025 52.096  0.826 |  0.009
## 1977      |     0.068 | -0.012  5.894  0.326 |  0.011 10.236  0.279 | -0.005
## 1978      |     0.148 | -0.027 31.464  0.800 | -0.004  1.802  0.023 | -0.010
##              ctr   cos2  
## 1973       4.384  0.032 |
## 1974      27.497  0.221 |
## 1975      20.961  0.094 |
## 1976      18.552  0.120 |
## 1977       5.256  0.058 |
## 1978      23.350  0.120 |
## 
## Columns (the 10 first)
##             Iner*1000    Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3
## January   |     0.044 |  0.023 11.075  0.937 |  0.003  0.370  0.015 | -0.005
## February  |     0.121 |  0.030 16.116  0.501 |  0.024 22.211  0.341 |  0.014
## March     |     0.034 |  0.019  7.094  0.784 | -0.001  0.015  0.001 |  0.006
## April     |     0.020 | -0.007  0.899  0.166 | -0.003  0.487  0.044 |  0.013
## May       |     0.061 |  0.018  7.507  0.465 |  0.003  0.411  0.013 | -0.019
## June      |     0.032 |  0.014  4.462  0.523 | -0.012  6.752  0.391 | -0.001
## July      |     0.045 | -0.007  1.162  0.098 |  0.015 11.713  0.487 | -0.007
## August    |     0.017 |  0.001  0.054  0.012 | -0.011  5.832  0.655 |  0.000
## September |     0.064 | -0.015  4.658  0.275 | -0.020 17.147  0.499 | -0.004
## October   |     0.031 | -0.014  4.273  0.526 | -0.008  2.642  0.161 |  0.008
##              ctr   cos2  
## January    2.166  0.037 |
## February  18.900  0.118 |
## March      3.794  0.085 |
## April     18.438  0.687 |
## May       40.804  0.509 |
## June       0.148  0.003 |
## July       5.902  0.100 |
## August     0.004  0.000 |
## September  1.526  0.018 |
## October    7.368  0.183 |

Focus sur la part de variance expliquée

On peut reprendre le code précédent, pour l’appliquer directement et voir la part d’inertie expliquée par chaque axe.

Le premier axe explique un peu plus de la moitié de l’information, et les deux premiers permettent d’expliquer 80% en tout.

eig = tibble(data.frame(afc$eig) %>% rownames_to_column("comp"))
eig
## # A tibble: 5 × 4
##   comp  eigenvalue percentage.of.variance cumulative.percentage.of.variance
##   <chr>      <dbl>                  <dbl>                             <dbl>
## 1 dim 1  0.000376                   53.8                               53.8
## 2 dim 2  0.000186                   26.5                               80.3
## 3 dim 3  0.0000758                  10.8                               91.1
## 4 dim 4  0.0000377                   5.38                              96.5
## 5 dim 5  0.0000246                   3.51                             100

La notion de seuil est moins pertinent en AFC, et on se contentera très souvent (voire tout le temps) de ne garder que deux axes.

ggplot(eig, aes(comp, percentage.of.variance)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  labs(x = "Composantes", y = "% de variance expliquée") +
  theme_minimal()

Représentation des modalités

Dans une AFC, le but est de représenter graphiquement les modalités de la variable en ligne (ici les années), ainsi que les modalités de la variable en colonnes (ici les mois). On peut créer les graphiques séparément (en rendant l’autre variable invisible).

On a une opposition 1975 vs 1978 sur l’axe 1, et une opposition 1976 vs 1973/1974 sur l’axe 2.

plot(afc, invisible = "col")

Le mois de décembre semble avoir un comportement différent des autres mois, ainsi que Février et Juillet, dans une moindre mesure.

plot(afc, invisible = "row")

Mais l’idéal est de représenter les deux directement sur le même graphique afin de relier les modalités des variables entre elles.

Les mois de Février et Décembre 1976 ont l’air particuliers, ainsi que le mois de juillet 1977 (et 1976).

plot(afc)

### Si on veut aller plus loin sur ce jeu de données

A FAIRE

Données Wine – ACP

Nous allons travailler sur des données concernant 3 types de vin. Elles sont disponibles sur cette page de l’UCI MLR. Il s’agit de 178 vins, réparties en 3 classes donc, et décrit par 13 variables quantitatives (lire la description dans le fichier wine.names pour plus d’informations).

Le code suivant permet de charger les données, et de nommer correctement les variables.

wine = read_csv("https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data", col_names = FALSE)
## Rows: 178 Columns: 14
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (14): X1, X2, X3, X4, X5, X6, X7, X8, X9, X10, X11, X12, X13, X14
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
names(wine) = c("class", "Alcohol", "Malic acid", "Ash", "Alcalinity of ash", "Magnesium", 
                "Total phenols", "Flavanoids", "Nonflavanoid phenols", "Proanthocyanins", 
                "Color intensity", "Hue", "OD280/OD315 of diluted wines", "Proline")
datatable(wine, options = list(scrollX = TRUE))

Vous devez donc réaliser les étapes suivantes :

  • Décrire les données
  • Réaliser une ACP centrée ou normée (choix à justifier)
  • Produire les graphiques nécessaires à l’interprétation
  • Identifier les classes sur le plan factoriel
  • Que peut-on dire globalement ?

Données CSP 1968 et 2014 – AFC

Nous disposons sur cette page de la répartition de la population active de 25 à 54 ans selon la catégorie socioprofessionnelle et la position vis à vis de l’emploi, par commune et département, de 1968 et de 2014.

Les codes ci-dessous vous montrent comment les importer (en supprimant les données sur les chômeurs).

regions1968 = read_csv("https://fxjollois.github.io/donnees/regions-csp/regions-csp-1968.csv") %>%
  select(!contains("Chômeurs")) %>%
  column_to_rownames("Région")
## Rows: 13 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): Région
## dbl (12): Agriculteurs Actifs ayant un emploi RP1968, Agriculteurs Chômeurs ...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
names(regions1968) = sub(" RP1968", "", names(regions1968))
datatable(regions1968, options = list(scrollX = TRUE))
regions2014 = read_csv("https://fxjollois.github.io/donnees/regions-csp/regions-csp-2014.csv") %>%
  select(!contains("Chômeurs")) %>%
  column_to_rownames("Région")
## Rows: 13 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): Région
## dbl (12): Agriculteurs Actifs ayant un emploi RP2014, Agriculteurs Chômeurs ...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
names(regions2014) = sub(" RP2014", "", names(regions2014))
datatable(round(regions2014), options = list(scrollX = TRUE))
  • Réaliser une AFC sur le croisement entre la catégorie socioprofessionnelle et la région
    • pour 1968
    • pour 2014
  • A l’aide des outils à votre disposition (tableaux et graphiques),
    • Décrire la situation pour chaque année,
    • Comparer les deux situations