M1 BBS - TP Maths - Analyse en Composantes Principales

1 ACP : Premier exemple illustrant le changement de repère et la variance portée par les axes

1.1 Jeu de données et préliminaires

Jeu de données utilisé en TDB TP2 pour la régression linéaire :

library(tidyverse)

## ── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0      ✔ purrr   0.3.5 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.4.1 
## ✔ readr   2.1.3      ✔ forcats 0.5.2 
## ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

genomes=read.table("http://silico.biotoul.fr/site/images/d/de/Bacterial_genomes.txt", sep="\t", header=TRUE)

Exploration du contenu

summary(genomes)

##   Genome_size      ORF_number  
##  Min.   :  920   Min.   : 744  
##  1st Qu.: 2250   1st Qu.:2310  
##  Median : 3608   Median :3475  
##  Mean   : 4023   Mean   :3831  
##  3rd Qu.: 4842   3rd Qu.:4634  
##  Max.   :10467   Max.   :9501

La taille des génomes en kilobases et le nombre de gènes/ORF sont du même ordre. De manière très grossière, on observe généralement 1 gène par kilobase. Variances :

apply(genomes, 2, var)

## Genome_size  ORF_number 
##     4842536     3621794

Visualisation sous forme de boites à moustaches:

par(mfrow=c(1,2))
boxplot(genomes$Genome_size, ylim=c(0,10**4), main='Genome size')
boxplot(genomes$ORF_number, ylim=c(0,10**4), main='ORF number')

Les variances sont-elles homogènes ?

var.test(genomes$Genome_size, genomes$ORF_number)

## 
##  F test to compare two variances
## 
## data:  genomes$Genome_size and genomes$ORF_number
## F = 1.3371, num df = 98, denom df = 98, p-value = 0.1523
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.8977926 1.9912333
## sample estimates:
## ratio of variances 
##           1.337054

Dans ce cas, on pourrait réaliser une ACP non normée. En pratique, pour ce jeu de données, ça ne fait presque aucune différence. On va faire avec la normalisation pour décomposer l’ensemble des opérations.

1.2 ACP normée avec `princomp`

Réalisation de l’ACP normée avec la fonction princomp :

acp = princomp(genomes, cor=T)

Visualisation et comparaison du repère original vs. projection dans celui des 2 premières composantes de l’ACP (les coordonnées des individus dans ce repère se trouvent dans scores) :

par(mfrow=c(1,2))

plot(genomes, pch="", main="Scatterplot + régression linéaire\n + centre du nuage")
abline(lm(genomes$ORF_number ~ genomes$Genome_size), lty='dotted', col='blue')
points(mean(genomes$Genome_size), mean(genomes$ORF_number), pch='×', col='red', cex=3)
text(genomes$Genome_size, genomes$ORF_number, labels= 1:nrow(genomes), cex=.5)

#plot(acp$scores[,1], -acp$scores[,2], xlab='PC1', ylab='PC2', pch="", ylim=2*c(-2000,2000), main='ACP')
plot(acp$scores[,1], -acp$scores[,2], xlab='PC1', ylab='PC2', pch="", ylim=c(-2,3), main='ACP')
text(acp$scores[,1], -acp$scores[,2], labels=1:nrow(genomes), cex=.5)
abline(h=0, lty="dashed")
abline(v=0, lty="dashed")

Variances portées par les axes via summary

summary(acp)

## Importance of components:
##                           Comp.1     Comp.2
## Standard deviation     1.3987612 0.20848734
## Proportion of Variance 0.9782665 0.02173348
## Cumulative Proportion  0.9782665 1.00000000

Discuter de la variance portée sur les deux premières composantes principales.

De manière visuelle avec la fonction screeplot :

screeplot(acp)

1.3 Analyse du résultat de `princomp`

Il est toujours conseillé de consulter l’aide (même si parfois, elle n’aide pas vraiment) pour savoir ce qui est renvoyé par une fonction :

?princomp

Pour ce qui est renvoyé, on s’intéressera à la section Value :

Value
princomp returns a list with class "princomp" containing the following components:

sdev    the standard deviations of the principal components.

loadings    the matrix of variable loadings (i.e., a matrix whose columns contain the eigenvectors). This is of class "loadings": see loadings for its print method.

center  the means that were subtracted.

scale   the scalings applied to each variable.

n.obs   the number of observations.

scores  if scores = TRUE, the scores of the supplied data on the principal components. These are non-null only if x was supplied, and if covmat was also supplied if it was a covariance list. For the formula method, napredict() is applied to handle the treatment of values omitted by the na.action.

call    the matched call.

na.action   If relevant.

Pour les initié·e·s, on pourra noter :

sdev : les écart-types des composantes principales (correspondent également aux valeurs propres de la matrice de covariance)
loadings : les vecteurs propres de la matrice de covariance
scores : les coordonnées individus (déjà utilisé plus tôt)
center et scale : valeurs utilisées pour centrer et réduire les variables

1.3.1 Rappel sur les étapes d’une ACP

Centrage des individus (changement d’origine pour la placer au centre du nuage de points)
Optionnel : Normalisation des variables
Calcul de la matrice de variance-covariance
Calcul des vecteurs et valeurs propres associés de la matrice de covariance

Les vecteurs propres (composantes principales) correspondent à la combinaison linéaire des variables de départ. Les valeurs propres associées correspondent aux variances de ces nouveaux axes.

1.3.2 Centrage

apply(genomes, 2, mean)

## Genome_size  ORF_number 
##    4023.051    3831.273

acp$center

## Genome_size  ORF_number 
##    4023.051    3831.273

Ainsi, la matrice centrée s’obtient :

Uc= apply(genomes, 2, function(x) x-mean(x)) %>% 
  as.matrix
head(Uc)

##         Genome_size   ORF_number
## Aaci     -1005.0505  -799.272727
## AbauD01   -119.0505     9.727273
## Acap       103.9495  -339.272727
## Acel     -1580.0505 -1565.272727
## Acfe     -1865.0505 -1739.272727
## Amir      4224.9495  3345.727273

Les moyennes sont à 0 :

apply(Uc, 2, mean)

##   Genome_size    ORF_number 
##  1.609487e-13 -2.251392e-13

1.3.3 Réduction

Pour la normalisation (des variables), on calcule les écart-types (au niveau de la population) soit la formule \(\sigma = \sqrt{ \frac{\sum_{i=1}^n (x_i - \overline{x})^2}{n}}\) avec ici, les \(\overline{x}=0\), on obtient \(\sigma = \sqrt{ \frac{\sum_{i=1}^n (x_i \cancel{- \overline{x}})^2}{n}} = \sqrt{ \frac{\sum_{i=1}^n x_i^2}{n}}\) que l’on peut écrire sous forme matricielle : \(\sqrt{\frac{1}{n}X^tX}\) :

sd.n = function(x) sqrt( 1/length(x) * t(x) %*% x)

Appliqué à chacune des variables

apply(Uc, 2, sd.n )

## Genome_size  ORF_number 
##    2189.434    1893.465

acp$scale

## Genome_size  ORF_number 
##    2189.434    1893.465

Ainsi, la matrice (déjà centrée) réduite s’obtient :

Ucn = t( t(Uc) / acp$scale )
head(Ucn)

##         Genome_size   ORF_number
## Aaci    -0.45904579 -0.422121681
## AbauD01 -0.05437501  0.005137286
## Acap     0.04747779 -0.179180859
## Acel    -0.72167074 -0.826670964
## Acfe    -0.85184136 -0.918565970
## Amir     1.92969934  1.766986379

Les écarts-types sont de 1 :

apply(Ucn, 2, sd.n )

## Genome_size  ORF_number 
##           1           1

1.3.4 Variance-covariance

mcov = 1/nrow(Ucn) * t(Ucn) %*% Ucn
mcov

##             Genome_size ORF_number
## Genome_size    1.000000   0.956533
## ORF_number     0.956533   1.000000

1.3.5 Vecteurs propres et valeurs propres

e = eigen(mcov)

Les vecteurs propres sont dans loadings :

e$vectors

##           [,1]       [,2]
## [1,] 0.7071068 -0.7071068
## [2,] 0.7071068  0.7071068

acp$loadings

## 
## Loadings:
##             Comp.1 Comp.2
## Genome_size  0.707  0.707
## ORF_number   0.707 -0.707
## 
##                Comp.1 Comp.2
## SS loadings       1.0    1.0
## Proportion Var    0.5    0.5
## Cumulative Var    0.5    1.0

Et les valeurs propres dans sdev (à mettre au carré pour avoir les variances) :

e$values

## [1] 1.95653303 0.04346697

acp$sdev**2

##     Comp.1     Comp.2 
## 1.95653303 0.04346697

Les proportions de variances :

e$values / sum(e$values)

## [1] 0.97826652 0.02173348

acp$sdev**2 / sum(acp$sdev**2)

##     Comp.1     Comp.2 
## 0.97826652 0.02173348

summary(acp)

## Importance of components:
##                           Comp.1     Comp.2
## Standard deviation     1.3987612 0.20848734
## Proportion of Variance 0.9782665 0.02173348
## Cumulative Proportion  0.9782665 1.00000000

1.3.6 Projection dans le nouveau repère

Il reste à calculer les coordonnées des individus dans le nouveau repère en utilisant les combinaisons linéaires (des variables de départ) fournie par les vecteurs propres.

Pour le 1er individu surla 1ère composante :

Ucn[1,]

## Genome_size  ORF_number 
##  -0.4590458  -0.4221217

Avec le 1er vecteur :

e$vectors[,1]

## [1] 0.7071068 0.7071068

Sa coordonnée sur le 1er axe est donc : -0.4590458 * 0.7071068 + -0.4221217 * 0.7071068 = -0.6230795

Et donc pour tous les individus :

Ucnp = Ucn %*% e$vectors
head(Ucnp)

##                [,1]        [,2]
## Aaci    -0.62307950  0.02610929
## AbauD01 -0.03481633  0.04208155
## Acap    -0.09312813 -0.16027187
## Acel    -1.09484292 -0.07424637
## Acfe    -1.25186703 -0.04718142
## Amir     2.61395154 -0.11505544

head(acp$scores)

##              Comp.1      Comp.2
## Aaci    -0.62307950 -0.02610929
## AbauD01 -0.03481633 -0.04208155
## Acap    -0.09312813  0.16027187
## Acel    -1.09484292  0.07424637
## Acfe    -1.25186703  0.04718142
## Amir     2.61395154  0.11505544

Visualisation

par(mfrow=c(1,2))
plot(acp$scores[,1], -acp$scores[,2], xlab='PC1', ylab='PC2', pch="", main='ACP princomp')
text(acp$scores[,1], -acp$scores[,2], labels=1:nrow(genomes), cex=.4)
abline(h=0, lty="dashed")
abline(v=0, lty="dashed")

plot(Ucnp, , xlab='PC1', ylab='PC2', pch="", main='ACP pas à pas')
text(Ucnp, labels=1:nrow(genomes), cex=.4)
abline(h=0, lty="dashed")
abline(v=0, lty="dashed")

2 Application sur le jeu de données iris

Dans ce TP, nous utilisons le jeu de données sur les iris. Description sur wikipedia http://en.wikipedia.org/wiki/Iris_flower_data_set

Et dans R, avec la commande

?iris

Chargement des données

iris %>%
  as_tibble

Bien sûr, on se posera certainement les questions habituelles : combien de variables ? de quelles natures ? combien d’individus ? pour chaque variétés, … Pour cela, ggpairs permet un apperçu rapide (s’il n’y a pas trop de variables) :

library(GGally)
ggpairs(iris, progress=F, aes(color=Species, alpha=.2))

Réaliser une ACP sur le jeu de données iris pas à pas (sans utiliser princomp) et discuter les résultats obtenus.

On factorise le code précédent sous forme de fonction :

pca = function(df, normalize=FALSE) {
  n = nrow(df)
  # centered
  U = as.matrix( apply(df, 2, function(x) x-mean(x)) )
  # normalized if asked
  scale = rep(1, ncol(df))
  if (normalize) {
    scale = apply(U, 2, function(x) sqrt( 1/length(x) * t(x) %*% x))
    U = t( t(U)/scale )
  }
  # variance-covariance & eigen
  mcov = 1/n * t(U) %*% U
  e = eigen(mcov)
  names(e$values) = paste0('Comp.', 1:length(e$values))
  colnames(e$vectors) = paste0('Comp.', 1:length(e$values))
  rownames(e$vectors) = colnames(df)
  # projection
  Up = U %*% e$vectors
  colnames(Up) = paste0('Comp.', 1:length(e$values))
  # format results as a list for returning a princomp list
  res = list(
    sdev = sqrt(e$values),
    loadings = e$vectors,
    center = apply(df, 2, mean),
    scale = scale,
    n.obs = n,
    scores = Up,
    call = sys.call()
  )
  class(res) = c('princomp')
  return(res)
}

acp.iris = pca(iris[,-5], normalize = T)

On devait obtenir ceci pour une ACP normée:

summary(acp.iris)

## Importance of components:
##                           Comp.1    Comp.2     Comp.3      Comp.4
## Standard deviation     1.7083611 0.9560494 0.38308860 0.143926497
## Proportion of Variance 0.7296245 0.2285076 0.03668922 0.005178709
## Cumulative Proportion  0.7296245 0.9581321 0.99482129 1.000000000

Discuter de la variance portée par les axes.

Illustrations

library(factoextra)

## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

fviz_screeplot(acp.iris, addlabels=T)

fviz_pca_biplot(acp.iris, habillage = iris$Species, addEllipses = T, label='var')

Même chose mais sans normalisation

acp.iris = pca(iris[,-5], normalize=F)

On devait obtenir ceci :

summary(acp.iris)

## Importance of components:
##                           Comp.1     Comp.2     Comp.3      Comp.4
## Standard deviation     2.0494032 0.49097143 0.27872586 0.153870700
## Proportion of Variance 0.9246187 0.05306648 0.01710261 0.005212184
## Cumulative Proportion  0.9246187 0.97768521 0.99478782 1.000000000

Discuter de la variance portée par les axes.

Illustrations

fviz_screeplot(acp.iris, addlabels=T)

fviz_pca_biplot(acp.iris, habillage = iris$Species, addEllipses = T, label='var')

Vous pouvez factoriser le code de la partie précédente sous forme de fonction qui prendra en paramètre la matrice de données quantitatives et un paramètre indiquant si la normalisation doit être effectuée.

Pour que summary ( et les fonctions screeplot biplot, …) fonctionne comme sur le résultat de princomp, il faut que cette fonction renvoie une variable ayant le même type de structures de données (et de même structure) :

class(acp)

## [1] "princomp"

names(acp)

## [1] "sdev"     "loadings" "center"   "scale"    "n.obs"    "scores"   "call"

pca = function(df, normalize=F) {
  ...
  # formattage du résultat comme princomp
  results = list( 
    ... ,
    call = sys.call()
    )
  class(results) = "princomp"
  return(results)
}

3 MDS : Multi-dimensional scaling

L’ACP s’utilise sur une matrice individus-variables (quantitatives).

Il arrive que l’on ne dispose pas de cette matrice mais seulement de la distance/similitude entre les individus.

On a alors une matrice dont les individus sont à la fois les lignes et les colonnes, et les valeurs correspondent à la distance/similitude entre une paire d’individus.

Simulons cette situation avec le jeu de données iris. On utilise dist qui par défaut utilise la distance euclidienne :

m = dist(iris[,-5])
as.matrix(m)[1:10, 1:10] %>% round(3)

##        1     2     3     4     5     6     7     8     9    10
## 1  0.000 0.539 0.510 0.648 0.141 0.616 0.520 0.173 0.922 0.469
## 2  0.539 0.000 0.300 0.332 0.608 1.091 0.510 0.424 0.510 0.173
## 3  0.510 0.300 0.000 0.245 0.510 1.086 0.265 0.412 0.436 0.316
## 4  0.648 0.332 0.245 0.000 0.648 1.166 0.332 0.500 0.300 0.316
## 5  0.141 0.608 0.510 0.648 0.000 0.616 0.458 0.224 0.922 0.529
## 6  0.616 1.091 1.086 1.166 0.616 0.000 0.995 0.700 1.459 1.010
## 7  0.520 0.510 0.265 0.332 0.458 0.995 0.000 0.424 0.548 0.480
## 8  0.173 0.424 0.412 0.500 0.224 0.700 0.424 0.000 0.787 0.332
## 9  0.922 0.510 0.436 0.300 0.922 1.459 0.548 0.787 0.000 0.557
## 10 0.469 0.173 0.316 0.316 0.529 1.010 0.480 0.332 0.557 0.000

MDS avec la commande cmdscale :

par(mfrow=c(1,2))
mds = cmdscale(m)
head(mds)

##           [,1]       [,2]
## [1,] -2.684126  0.3193972
## [2,] -2.714142 -0.1770012
## [3,] -2.888991 -0.1449494
## [4,] -2.745343 -0.3182990
## [5,] -2.728717  0.3267545
## [6,] -2.280860  0.7413304

plot(mds, col=iris$Species, pch=19, main='MDS')
plot(acp.iris$scores[,1], -acp.iris$scores[,2], col=iris$Species, pch=19, main='ACP')

Visualisation

library(gridExtra)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

mds.p = tibble(Dim1=mds[,1], Dim2= mds[,2], Species = iris$Species) %>%
  ggplot(aes(x=Dim1, y=Dim2, color = Species, shape=Species)) +
  geom_point(show.legend=F) +
  theme_light() +
  ggtitle('MDS')

acp.p = tibble(PC1=acp.iris$scores[,1], PC2= -acp.iris$scores[,2], Species = iris$Species) %>%
  ggplot(aes(x=PC1, y=PC2, color = Species, shape=Species)) +
  geom_point() +
  theme_light() +
  ggtitle('ACP')
grid.arrange(mds.p, acp.p, ncol=2)

4 Pour aller plus loin

Pour la réduction de dimensions et la visualisation, on pourra s’intéresser à d’autres méthodes selon les données et leur nature

qualitatives : Analyse factorielle des correspondances (AFC), et analyse des correspondances multiples (ACM)
UMAP : données volumimneuses, pseudo-distances locales

M1 BBS - TP Maths - Analyse en Composantes Principales