ACP dans R: prcomp vs princomp

kassambara | 15/10/2017 | 37029 | Commentaires (3) | Méthodes des Composantes Principales dans R: Guide Pratique

Ce tutoriel R décrit comment faire l’Analyse en Composantes Principales (ACP) en utilisant les fonctions de base de R: prcomp() et princomp(). Vous apprendrez comment prédire les coordonnées de nouveaux individus et variables en utilisant l’ACP. Nous expliquerons également la théorie derrière les résultats de l’ACP.

Apprendre plus sur l’ACP: ACP - Analyse en Composantes Principales avec R: L’Essentiel.

ACP dans R: prcomp vs princomp

Contents:

Méthodes générales concernant l’ACP
Fonctions prcomp() et princomp()
Package pour la visualisation de l’ACP
Données
ACP dand R avec prcomp()
Accéder aux résultats de l’ACP
Prediction avec l’ACP
- Individus supplémentaires
- Variables supplémentaires
Théory des résultats de l’ACP
- Résultats de l’ACP pour les variables
- Résultats de l’ACP pour les individus

Livre (en anglais):

Practical Guide to Principal Component Methods in R

Méthodes générales concernant l’ACP

Il existe deux méthodes générales pour effectuer l’ACP dans R:

La “décomposition spectrale” qui examine les covariances / corrélations entre les variables
La “décomposition de la valeur singulière” qui examine les covariances / corrélations entre les individus

La fonction princomp() utilise l’approche de la décomposition spectrale. Les fonctions prcomp() et PCA() [FactoMineR] utilisent la décomposition en valeur singulière.

Selon la documentation de R, la décomposition en valeur singulière a une meilleure précision numérique. Par conséquent, la fonction prcomp() est préférée par rapport à princomp().

Fonctions prcomp() et princomp()

Format simplifié des deux fonctions:

prcomp(x, scale = FALSE)
princomp(x, cor = FALSE, scores = TRUE)

Arguments pour prcomp():

x: matrice ou data frame
scale: une valeur logique indiquant si les variables doivent être standardisées

Arguments pour princomp():

x: matrice ou data frame
cor: une valeur logique. Si TRUE, les données seront standardisées avant l’analyse
scores: une valeur logique. Si TRUE, les coordonnées de chaque composante principale sont calculées

Résultats renvoyés par les deux fonctions:

prcomp() nom	princomp() nom	Description
sdev	sdev	l’écart type des composantes principales
rotation	loadings	Matrice. loading des variables (les colonnes sont les vecteurs propres)
center	center	Moyenne des variables
scale	scale	Ecart type des variables
x	scores	Coordonnées des individus (observations).

Dans les sections suivantes nous allons nous focaliser sur la fonction prcomp()

Package pour la visualisation de l’ACP

Package R: factoextra

Installation:

install.packages("factoextra")

Chargez le package:

library(factoextra)

Données

Jeu de données: decathlon2 [factoextra].

Contenu des données:

Individus actifs (lignes 1 à 23) et variables actives (colonnes 1 à 10). Utilisés pour effectuer l’analyse en composantes principales
Individus supplémentaires (lignes 24 à 27) et variables supplémentaires (colonnes 11 à 13). Coordonnées prédites à l’aide de l’information de l’ACP et des paramètres obtenus avec les individus / variables actifs.

Analyse en composantes principales. Format des données

Extraction des individus et variables actifs:

library("factoextra")
data(decathlon2)
decathlon2.active <- decathlon2[1:23, 1:10]
head(decathlon2.active[, 1:6])

##           X100m Long.jump Shot.put High.jump X400m X110m.hurdle
## SEBRLE     11.0      7.58     14.8      2.07  49.8         14.7
## CLAY       10.8      7.40     14.3      1.86  49.4         14.1
## BERNARD    11.0      7.23     14.2      1.92  48.9         15.0
## YURKOV     11.3      7.09     15.2      2.10  50.4         15.3
## ZSIVOCZKY  11.1      7.30     13.5      2.01  48.6         14.2
## McMULLEN   10.8      7.31     13.8      2.13  49.9         14.4

ACP dand R avec prcomp()

Chargement du package factoextra pour la visualisation

library(factoextra)

Calculer l’ACP

res.pca <- prcomp(decathlon2.active, scale = TRUE)

Visualiser les valeurs propres. Montre le pourcentage de variances expliquées par chaque axe principal.

fviz_eig(res.pca)

Graphique des individus. Coloration en fonction du cos2 (qualité de représentation). Les individus similaires sont groupés ensemble.

fviz_pca_ind(res.pca,
             col.ind = "cos2", # Colorer par le cos2
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE     
             )

Graphique des variables. Coloration en fonction de la contribution des variables. Les variables corrélées positivement sont du même côté du graphique. Les variables corrélées négativement sont sur des côtés opposés du graphique.

fviz_pca_var(res.pca,
             col.var = "contrib", 
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE     
             )

Biplot des individus et des variables

fviz_pca_biplot(res.pca, repel = TRUE,
                col.var = "#2E9FDF", 
                col.ind = "#696969"  
                )

Accéder aux résultats de l’ACP

library(factoextra)
# Valeurs propres
eig.val <- get_eigenvalue(res.pca)
eig.val
  
# Résultats des variables
res.var <- get_pca_var(res.pca)
res.var$coord          # Coordonnées
res.var$contrib        # Contributions aux axes
res.var$cos2           # Qualité de représentation 
# Résultats des individus
res.ind <- get_pca_ind(res.pca)
res.ind$coord          # Coordonnées
res.ind$contrib        # Contributions aux axes
res.ind$cos2           # Qualité de représentation

Prediction avec l’ACP

Dans cette section, nous allons montrer comment prédire les coordonnées des individus et des variables supplémentaires en utilisant uniquement les informations fournies par l’ACP précédemment effectuée.

Individus supplémentaires

Données: lignes 24 à 27 et colonnes 1 à 10 [jeu de données decathlon2]. Les nouvelles données doivent contenir les colonnes (variables) portant les mêmes noms et dans le même ordre que les données actives utilisées pour calculer l’ACP.

# Individus supplémentaires
ind.sup <- decathlon2[24:27, 1:10]
ind.sup[, 1:6]

##         X100m Long.jump Shot.put High.jump X400m X110m.hurdle
## KARPOV   11.0      7.30     14.8      2.04  48.4         14.1
## WARNERS  11.1      7.60     14.3      1.98  48.7         14.2
## Nool     10.8      7.53     14.3      1.88  48.8         14.8
## Drews    10.9      7.38     13.1      1.88  48.5         14.0

Redire les coordonnées des individus. Fonction R: predict():

ind.sup.coord <- predict(res.pca, newdata = ind.sup)
ind.sup.coord[, 1:4]

##            PC1    PC2   PC3    PC4
## KARPOV   0.777 -0.762 1.597  1.686
## WARNERS -0.378  0.119 1.701 -0.691
## Nool    -0.547 -1.934 0.472 -2.228
## Drews   -1.085 -0.017 2.982 -1.501

Graphique des individus incluant les individus supplémentaires:

# Individus actifs
p <- fviz_pca_ind(res.pca, repel = TRUE)
# Ajouter les individus supplementaires
fviz_add(p, ind.sup.coord, color ="blue")

Les coordonnées prédites des individus peuvent être calculées manuellement comme suit:

Centrer-réduire les nouveaux individus en utilisant les informations de l’ACP
Calculer les coordonnées prédites en multipliant les valeurs centrées-réduites par les vecteurs propres des axes principaux.

Code R:

# Centrer-réduire les individus supplémentaires
ind.scaled <- scale(ind.sup, 
                    center = res.pca$center,
                    scale = res.pca$scale)
# Coordonnées des individus
coord_func <- function(ind, loadings){
  r <- loadings*ind
  apply(r, 2, sum)
}
pca.loadings <- res.pca$rotation
ind.sup.coord <- t(apply(ind.scaled, 1, coord_func, pca.loadings ))
ind.sup.coord[, 1:4]

##            PC1    PC2   PC3    PC4
## KARPOV   0.777 -0.762 1.597  1.686
## WARNERS -0.378  0.119 1.701 -0.691
## Nool    -0.547 -1.934 0.472 -2.228
## Drews   -1.085 -0.017 2.982 -1.501

Variables supplémentaires

Variables catégorielles

Le jeu de données decathlon2 contient une variable qualitative supplémentaire, la colonne 13 correspondant aux types de compétition.

Les variables qualitatives / catégorielles peuvent être utilisées pour colorer les individus par groupes. La variable catégorielle doit être de même longueur que le nombre d’individus actifs (ici 23).

groups <- as.factor(decathlon2$Competition[1:23])
fviz_pca_ind(res.pca,
             col.ind = groups, # colorer par groupes
             palette = c("#00AFBB",  "#FC4E07"),
             addEllipses = TRUE, # Ellipse de concentration
             ellipse.type = "confidence",
             legend.title = "Groups",
             repel = TRUE
             )

Calculer les coordonnées pour les différentes catégories de la variable catégorielle. Les coordonnées d’un groupe donné sont calculées comme les coordonnées moyennes des individus du groupe.

library(magrittr) # pour le pipe %>%
library(dplyr)   # tout le reste
# 1. Coordonnées des individus
res.ind <- get_pca_ind(res.pca)
# 2. Coordonnées des groupes
coord.groups <- res.ind$coord %>%
  as_data_frame() %>%
  select(Dim.1, Dim.2) %>%
  mutate(competition = groups) %>%
  group_by(competition) %>%
  summarise(
    Dim.1 = mean(Dim.1),
    Dim.2 = mean(Dim.2)
    )
coord.groups

## # A tibble: 2 x 3
##   competition Dim.1  Dim.2
##           
## 1    Decastar -1.31 -0.119
## 2    OlympicG  1.20  0.109

Variables quantitatives

Données: colonnes 11:12. Doit être de la même longueur que le nombre d’individus actifs (ici 23)

quanti.sup <- decathlon2[1:23, 11:12, drop = FALSE]
head(quanti.sup)

##           Rank Points
## SEBRLE       1   8217
## CLAY         2   8122
## BERNARD      4   8067
## YURKOV       5   8036
## ZSIVOCZKY    7   8004
## McMULLEN     8   7995

Les coordonnées d’une variable quantitative donnée sont calculées comme la corrélation entre les variables quantitatives et les composantes principales.

# Predire les coordonées et calculer le cos2
quanti.coord <- cor(quanti.sup, res.pca$x)
quanti.cos2 <- quanti.coord^2
# Graphique des variables incluants les variables supplementaires
p <- fviz_pca_var(res.pca)
fviz_add(p, quanti.coord, color ="blue", geom="arrow")

Théory des résultats de l’ACP

Résultats de l’ACP pour les variables

Calcul des résultats (coordonnées, cos2 et contributions) des variables:

var.coord = loadings * l’écart type des axes
var.cos2 = var.coord^2
var.contrib. En pourcentage : (var.cos2 * 100) / (cos2 total de l’axe)

# Fonction d'aide
#::::::::::::::::::::::::::::::::::::::::
var_coord_func <- function(loadings, comp.sdev){
  loadings*comp.sdev
}
# Coordonnées
#::::::::::::::::::::::::::::::::::::::::
loadings <- res.pca$rotation
sdev <- res.pca$sdev
var.coord <- t(apply(loadings, 1, var_coord_func, sdev)) 
head(var.coord[, 1:4])

##                 PC1     PC2    PC3     PC4
## X100m        -0.851  0.1794 -0.302  0.0336
## Long.jump     0.794 -0.2809  0.191 -0.1154
## Shot.put      0.734 -0.0854 -0.518  0.1285
## High.jump     0.610  0.4652 -0.330  0.1446
## X400m        -0.702 -0.2902 -0.284  0.4308
## X110m.hurdle -0.764  0.0247 -0.449 -0.0169

# Cos2
#::::::::::::::::::::::::::::::::::::::::
var.cos2 <- var.coord^2
head(var.cos2[, 1:4])

##                PC1      PC2    PC3      PC4
## X100m        0.724 0.032184 0.0909 0.001127
## Long.jump    0.631 0.078881 0.0363 0.013315
## Shot.put     0.539 0.007294 0.2679 0.016504
## High.jump    0.372 0.216424 0.1090 0.020895
## X400m        0.492 0.084203 0.0804 0.185611
## X110m.hurdle 0.584 0.000612 0.2015 0.000285

# Contributions
#::::::::::::::::::::::::::::::::::::::::
comp.cos2 <- apply(var.cos2, 2, sum)
contrib <- function(var.cos2, comp.cos2){var.cos2*100/comp.cos2}
var.contrib <- t(apply(var.cos2,1, contrib, comp.cos2))
head(var.contrib[, 1:4])

##                PC1     PC2   PC3     PC4
## X100m        17.54  1.7505  7.34  0.1376
## Long.jump    15.29  4.2904  2.93  1.6249
## Shot.put     13.06  0.3967 21.62  2.0141
## High.jump     9.02 11.7716  8.79  2.5499
## X400m        11.94  4.5799  6.49 22.6509
## X110m.hurdle 14.16  0.0333 16.26  0.0348

Résultats de l’ACP pour les individus

ind.coord = res.pca$x
Cos2 des individus. Deux étapes:
- Calculer le carré de la distance entre chaque individu et le centre de gravité de l’ACP: d2 = [(var1_ind_i - mean_var1)/sd_var1]^2 + …+ [(var10_ind_i - mean_var10)/sd_var10]^2 + …+..
- Cosinus carré des individus: ind.coord^2/d2
Contributions des individus: 100 * (1 / number_of_individuals)*(ind.coord^2 / comp_sdev^2).

# Coordonnées
#::::::::::::::::::::::::::::::::::
ind.coord <- res.pca$x
head(ind.coord[, 1:4])

##              PC1    PC2    PC3     PC4
## SEBRLE     0.191 -1.554 -0.628  0.0821
## CLAY       0.790 -2.420  1.357  1.2698
## BERNARD   -1.329 -1.612 -0.196 -1.9209
## YURKOV    -0.869  0.433 -2.474  0.6972
## ZSIVOCZKY -0.106  2.023  1.305 -0.0993
## McMULLEN   0.119  0.992  0.844  1.3122

# Cos2 des individus
#:::::::::::::::::::::::::::::::::
# 1. Carré de la distance entre chaque individu 
# et le centre de gravité de l'ACP
center <- res.pca$center
scale<- res.pca$scale
getdistance <- function(ind_row, center, scale){
  return(sum(((ind_row-center)/scale)^2))
  }
d2 <- apply(decathlon2.active,1,getdistance, center, scale)
# 2. Cos2. La somme de chaque ligne est 1
cos2 <- function(ind.coord, d2){return(ind.coord^2/d2)}
ind.cos2 <- apply(ind.coord, 2, cos2, d2)
head(ind.cos2[, 1:4])

##               PC1    PC2     PC3     PC4
## SEBRLE    0.00753 0.4975 0.08133 0.00139
## CLAY      0.04870 0.4570 0.14363 0.12579
## BERNARD   0.19720 0.2900 0.00429 0.41182
## YURKOV    0.09611 0.0238 0.77823 0.06181
## ZSIVOCZKY 0.00157 0.5764 0.23975 0.00139
## McMULLEN  0.00218 0.1522 0.11014 0.26649

# Contributions des individus
#:::::::::::::::::::::::::::::::
contrib <- function(ind.coord, comp.sdev, n.ind){
  100*(1/n.ind)*ind.coord^2/comp.sdev^2
}
ind.contrib <- t(apply(ind.coord, 1, contrib, 
                       res.pca$sdev, nrow(ind.coord)))
head(ind.contrib[, 1:4])

##              PC1    PC2    PC3     PC4
## SEBRLE    0.0385  5.712  1.385  0.0357
## CLAY      0.6581 13.854  6.460  8.5557
## BERNARD   1.8627  6.144  0.135 19.5783
## YURKOV    0.7969  0.443 21.476  2.5794
## ZSIVOCZKY 0.0118  9.682  5.975  0.0523
## McMULLEN  0.0148  2.325  2.497  9.1353

2 Notes

Enjoyed this article? Give us 5 stars (just above this text block)! Reader needs to be STHDA member for voting. I’d be very grateful if you’d help it spread by emailing it to a friend, or sharing it on Twitter, Facebook or Linked In.

Show me some love with the like buttons below... Thank you and please don't forget to share and comment below!!

Avez vous aimé cet article? Donnez nous 5 étoiles (juste au dessus de ce block)! Vous devez être membre pour voter. Je vous serais très reconnaissant si vous aidiez à sa diffusion en l'envoyant par courriel à un ami ou en le partageant sur Twitter, Facebook ou Linked In.

Montrez-moi un peu d'amour avec les like ci-dessous ... Merci et n'oubliez pas, s'il vous plaît, de partager et de commenter ci-dessous!