Analyse en Composantes Principales

Marie-Pierre Etienne

ENSAI - CREST

https://marieetienne.github.io/MAF/

2024-12-19

Et on commence avec les manchots

Rappel des données sur les manchots

Données disponibles dans le package palmerpenguins mises à disposition par le Dr. Kristen Gorman et the Palmer Station, Antarctica LTER.

On a mesuré les caractéristiques morphologiques de divers manchots :

Les 6 premières lignes (parmi 333 )

# A tibble: 333 × 4
  bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
           <dbl>         <dbl>             <int>       <int>
1           39.1          18.7               181        3750
2           39.5          17.4               186        3800
3           40.3          18                 195        3250
4           36.7          19.3               193        3450
5           39.3          20.6               190        3650
6           38.9          17.8               181        3625
# ℹ 327 more rows
  • bill_length_mm : la longueur du bec,

  • bill_depth_mm : l’épaisseur du bec,

  • flipper_length_mm : la longueur de la nageoire,

  • body_mass_g : le poids du corps.

Comment visualiser au mieux ces données pour faire apparaître les liens entre variables et identifier des resemblances entre individus ?

Voir c’est comprendre : comment représenter l’information contenue dans ce tableau ?

Idée 1 : on représente les liens des variables 2 à 2

Voir c’est comprendre : comment représenter l’information contenue dans ce tableau ?

Idée 1 : dangereuse

Voir c’est comprendre : comment représenter l’information contenue dans ce tableau ?

Objectifs

  • Représenter sans perdre trop d’information,
  • Quantifier la perte d’information,
  • Comprendre quelles sont les informations redondantes (variables liées),
  • Idéalement des individus éloignés dans le nuage initial, restent éloignés dans la représentation

Dans toute la suite, sans perte de généralité, on va considérer que les variables sont centrées.

Si on savait faire une Analyse en Composantes Principales

library(FactoMineR)
dta <- penguins %>% select(bill_length_mm, bill_depth_mm, flipper_length_mm, body_mass_g) 
penguins.PCA <- PCA(X = dta, ncp = 4, graph = FALSE)

L’ACP permet de trouver la projection en 2 dimensions qui préserve au mieux l’information contenue dans le nuage de points (avec peu de déformation).

plot(penguins.PCA, choix = "ind")

Reste à comprendre

  • Comment quantifier l’information contenue dans un nuage de points
  • Comment est construite cette meilleure projection
  • Quels sont les individus bien représentés ?

L’ACP définit de nouvelles variables, les composantes principales qui sont une combinaison linéaire des anciennes variables. On peut représenter les variables dans le premier plan de l’ACP.

plot(penguins.PCA, choix = "var")

Reste à comprendre

  • Comment interpréter cette représentation ?
  • Quelles sont les variables bien représentées ?
  • Que peut-on déduire de ce graphique sur le lien entre les variables.
  • Comment lier graphique des individus et graphique des variables

Les données

Mise en forme des données

L’ACP s’applique à des données quantitatives sous forme de tableau. Ce tableau est composé de n individus (lignes) et de p variables (colonnes).

Dans la suite on va noter \({\boldsymbol{X}}\) la matrice des données, \({\boldsymbol{X}}\in {\mathbb{R}}^{n\times p}\).

\(x_{\class{alea}{i}\class{orange}{j}}\) est la valeur de la variable \(\class{orange}{j}\) mesurée pour l’individu \(\class{alea}{i}.\)

\[{\boldsymbol{X}}= \overset{\color{orange}{\begin{matrix}var_1& \ldots \ & \ \ldots \ &\ \ldots\ & \ var_p \end{matrix}}}{\begin{pmatrix} x_{\class{alea}{1}}^{\class{orange}{1}} & \ldots & \ldots &\ldots & x_{\class{alea}{1}}^{\class{orange}{p}}\\ \vdots & & & &\vdots \\ & & x_{\class{alea}{i}}^{\class{orange}{j}} & & \\ \vdots & & & & \vdots \\ x_{\class{alea}{n}}^{\class{orange}{1}} & & & & \ldots x_{\class{alea}{n}}^{\class{orange}{p}}\\ \end{pmatrix}}\]

Deux points de vue complémentaires

Le nuage des individus \(C^n\)

On peut considérer qu’un individu \(i\) est un vecteur \(\class{alea}{\boldsymbol{x}_{i}}\) dans un espace de dimension \(p\). Par convention tous les vecteurs sont des vecteurs colonnes, donc on peut écrire

\[{\boldsymbol{X}}=\begin{pmatrix} \class{alea}{\boldsymbol{x}_{1}}^\top\\ \vdots \\ \class{alea}{\boldsymbol{x}_{n}}^\top\\ \end{pmatrix},\]

L’ensemble des \(n\) vecteurs forme le nuage des individus (ce qu’on représente classiquement).

Le nuage des variables \(C^p\)

On peut considérer qu’une variable \(j\) est un vecteur \(\class{orange}{\boldsymbol{x}^{j}}\) dans un espace de dimension \(n\) et on peut écrire

\[{\boldsymbol{X}}=\begin{pmatrix} \class{orange}{\boldsymbol{x}^{1}} & \ldots & \class{orange}{\boldsymbol{x}^{p}} \end{pmatrix},\]

L’ensemble des \(p\) vecteurs forme le nuage des variables.

Bien sûr les deux nuages sont intimement liés

Mesurer l’information contenue dans un nuage de points

L’inertie par rapport à un point

Definition 1

L’inertie (terme hérité de la physique) par rapport à un point \(a \in {\mathbb{R}}^{p}\) est définie par
\[I_{\boldsymbol{a}} = \frac{1}{n} \sum_{i=1}^n \lVert\boldsymbol{x}_{i} - \boldsymbol{a}\rVert^2 = \frac{1}{n} \sum_{i=1}^n\sum_{k=1}^p (x_i^k -a^k)^2\]

L’inertie par rapport à un point quantifie l’information perdue lorsque le jeu de données est résumé par ce point.

Quel est selon vous le point \(\boldsymbol{a}\) qui minimise l’inertie ?

\[argmin_{\boldsymbol{a}} I_{\boldsymbol{a}} = ??\]

BINGO !!!

Il s’agit du point \(\boldsymbol{g}=\boldsymbol{x_{\bullet}}=({x}_{\bullet}^1, \ldots, {x}_{\bullet}^p )^\top,\) avec \({x}_{\bullet}^k=\frac{1}{n} \sum_{i=1}^n x_{i}^k\), la valeur moyenne de la variable \(k\) sur l’ensemble des individus.

Puisque nous avons supposé que nos variables étaient centrées, \(G=O\) l’origine sur repère.

\[\forall \boldsymbol{a}\in\mathbb{R}^p,\; I_\boldsymbol{a}=I_\boldsymbol{g} + \lVert\boldsymbol{a}- \boldsymbol{g}\rVert^2.\]

Conséquence

Le meilleur résumé du nuage de points se résumant à un point est le barycentre du nuage de points, c’est à dire la moyenne sur chacune des dimensions.

Remarque

\[I_{\boldsymbol{g}} = I = \frac{1}{n} \sum_{i=1}^n \left( \sum_{k=1}^p \left( x_{i}^k - x_{\bullet}^k \right)^2 \right) = \sum_{k=1}^p \left( \frac{1}{n} \sum_{i=1}^n \left( x_{i}^k - x_{\bullet}^k \right)^2 \right) = \sum_{k=1}^p Var(\boldsymbol{x}^{k})\] Si les variables sont réduites

\[I = p\]

L’inertie par rapport à un axe

Definition 2

L’inertie par rapport à l’axe \(\Delta\) est définie par \[I_{\boldsymbol{\Delta}} = \frac{1}{n} \sum_{i=1}^n \lVert\boldsymbol{x}_{i} - \boldsymbol{p_i}^{\Delta}\rVert^2 = \frac{1}{n} \sum_{i=1}^n\sum_{k=1}^p (x_i^k -p^{k,\Delta}_i)^2, \quad \mbox{où }\ \boldsymbol{p_i}^{\Delta} \mbox{projeté orthogonal de } \boldsymbol{x_i} \mbox{ sur }\Delta.\]

L’inertie par rapport à l’axe \(\Delta\) quantifie l’information perdue lorsque l’on résume le nuage de points à son projeté sur \(\Delta\).

Identifier \(\Delta\) tel que \(I_{\Delta}\) soit minimale

Ce qui revient à l’axe \(\Delta\) qui assure que la projection sur \(\Delta\) déforme le moins possible le nuage de points.

Dans le cadre de l’ACP, on veut construire une ensemble d’axes orthogonaux (une nouvelle base de \({\mathbb{R}}^p\)), de façon à ce que

  • la projection sur le premier axe, soit le meilleur résumé du nuage de points en une dimension,
  • le deuxième axe est orthogonal au premier et représente “le deuxième meileur choix” …. dans un sens que l’on va définir
  • le troisième …

  • Produit scalaire : Soient deux élements \(\boldsymbol{x}\) et \(\boldsymbol{y}\) de \(\mathbb{R}^p\), le produit scalaire est noté \(<\boldsymbol{x},\boldsymbol{y}>\) et \(<\boldsymbol{x},\boldsymbol{y}>=\sum_{i=1}^p x_i y_i = \boldsymbol{x}^\top \boldsymbol{y}.\)

  • Orthogonalité Soient deux élements \(\boldsymbol{x}\) et \(\boldsymbol{y}\) de \(\mathbb{R}^p\), on dit que \(\boldsymbol{x}\) et \(\boldsymbol{y}\) sont orthogonaux (noté \(\boldsymbol{x}\perp \boldsymbol{y}\)) si \(<\boldsymbol{x},\boldsymbol{y}> = \boldsymbol{x}^\top \boldsymbol{y}=0\).

  • Projection orthogonale On note \(\boldsymbol{a}_u\) le projeté de \(\boldsymbol{a}\) sur la droite \(\Delta\) définit par son vecteur directeur unitaire \(\boldsymbol{u}\) et passant par l’origine \[\boldsymbol{a}_u = < \boldsymbol{a}, \boldsymbol{u} > \boldsymbol{u} = (\boldsymbol{a}^\top \boldsymbol{u})\, \boldsymbol{u}.\] Le vecteur \(\boldsymbol{a_u}\) est orthogonal à \(\boldsymbol{a}-\boldsymbol{a_u}\): \[< \boldsymbol{a_u}, \boldsymbol{a}-\boldsymbol{a_u}> = \left ((\boldsymbol{a}^\top \boldsymbol{u}) \boldsymbol{u}\right)^\top \left ( \boldsymbol{a} - (\boldsymbol{a}^\top \boldsymbol{u})\, \boldsymbol{u} \right)=(\boldsymbol{a}^\top \boldsymbol{u}) (\boldsymbol{u}^\top \boldsymbol{a}) - (\boldsymbol{a}^\top \boldsymbol{u})^2 \boldsymbol{u}^\top \boldsymbol{u} =0 \]

L’inertie par rapport à un sous espace affine

  • Espace orthogonal Soit \(E\) un sous-espace vectoriel de \(\mathbb{R}^p\). On definit \(E^\perp=\{\boldsymbol{x}\in \mathbb{R}^p:\; \forall \boldsymbol{y}\in E,\; \boldsymbol{x}\perp \boldsymbol{y}\}\).

Dans \({\mathbb{R}}^2\), on considère le sous espace vectoriel \(E\), tel que \[ E =\left \lbrace \boldsymbol{x}=\begin{pmatrix} x_1 & x_2 \end{pmatrix}^\top \in {\mathbb{R}}^2 : x_1 = x_2 \right \rbrace \] \(E\) est la droite d’équation \(y=x\).

Qui est \(E^\perp\) ?

Soit \(\boldsymbol{y}\in E^\perp,\) alors

\[\begin{align} 0 & = \boldsymbol{x}^\top \boldsymbol{y}\\ & = \begin{pmatrix} x_1 & x_2 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \end{pmatrix}\\ & = x_1 y_1 + x_2 y_2 \\ & = x_1 (y_1 + y_2) \quad (\boldsymbol{x}\in E \Rightarrow x_1 =x_2) \\ \end{align}\]

\[E^\perp = \left \lbrace \boldsymbol{y}=\begin{pmatrix} y_1 & y_2 \end{pmatrix}^\top \in {\mathbb{R}}^2 : y_1 + y_2 = 0 \right \rbrace \]

Definition 3

L’inertie par rapport à un sous espace vectoriel \(E\) est définie par \[I_{E} = \frac{1}{n} \sum_{i=1}^n \lVert\boldsymbol{x}_{i} - \boldsymbol{p_i}^{E}\rVert^2 = \frac{1}{n} \sum_{i=1}^n\sum_{k=1}^p (x_i^k -p^{k,E}_i)^2, \quad \mbox{où }\ \boldsymbol{p_i}^{E} = argmin_{\boldsymbol{y}\in E} d(\boldsymbol{y}, \boldsymbol{x}_{i}), \mbox{i.e. projeté orthogonal de } \boldsymbol{x_i} \mbox{ sur }E.\]

C’est l’inertie perdue lorsque l’on résume le nuage de points à sa projection sur le sous espace \(E\).

Décomposition de l’inertie

Soit \(E\) un sous espace vectoriel de \({\mathbb{R}}^p\), alors \[I = I_E + I_{E^\perp}.\]

Preuve

Un dessin vaut mieux qu’un long discours et Pythagore est ton ami !

\[d(\class{alea}{\boldsymbol{x}_{i}}, O)^2 = d(\boldsymbol{x}_{i}, \class{bleu}{\boldsymbol{x}_{i}^E})^2 + d(\class{alea}{\boldsymbol{x}_{i}}, \class{rouge}{\boldsymbol{x}_{i}^{E^\perp}})^2.\]

\[\begin{align} I & = \frac{1}{n} \sum_{i=1}^n d(\class{alea}{\boldsymbol{x}_{i}}, O)^2 \\ & = \frac{1}{n} \sum_{i=1}^n d(\class{alea}{\boldsymbol{x}_{i}}, \class{bleu}{\boldsymbol{x}_{i}^E})^2 + \frac{1}{n} \sum_{i=1}^n d(\class{alea}{\boldsymbol{x}_{i}}, \class{rouge}{\boldsymbol{x}_{i}^{E^\perp}})^2 \\ & = \class{bleu}{I_E} + \class{rouge}{I_{E^\perp}} \end{align}\]

Principe de la décomposition à venir

  1. Identifier le vecteur \(\boldsymbol{u_1}\) de \({\mathbb{R}}^p\), tel que \(I_{E_1}\) soit minimale avec \(E_1=\left \lbrace \lambda \boldsymbol{u_1}, \lambda \in {\mathbb{R}}\right\rbrace,\)

  2. Dans \(E_1^\perp\), identifier le vecteur \(\boldsymbol{u_2}\) de \({\mathbb{R}}^p\), tel que \(I_{E_2}\) soit minimale avec \(E_2=\left \lbrace \lambda \boldsymbol{u_2}, \lambda \in {\mathbb{R}}\right\rbrace,\)

  3. Dans \(\left ( E_1 \oplus E_2 \right)^\perp\), identifier le vecteur \(\boldsymbol{u_3}\) de \({\mathbb{R}}^p\), tel que ….

Finalement on aura la décomposition

\[I= I_{E_1} + I_{E_2} + \ldots + I_{E_p}, \quad \mbox{avec } I_{E_1} \leq I_{E_2} \leq \ldots \leq I_{E_p}\]

Construction de la base des composantes principales

Se donner des intuitions

Se donner des intuitions

Identifier \(\boldsymbol{u_1}\) - Formalisation

  • On cherche \(\boldsymbol{u_1}\) tel que \(I_{u_1}\) soit minimale mais \(I_{u_1} + I_{u_1^\perp}=I\). Minimiser \(I_{u_1}\) revient à maximiser \(I_{u_1^\perp}\), i.e l’inertie du nuage projeté.

    • Le premier individu se projette en \(\class{alea}{x_{1}}^\top \boldsymbol{u_1}\) sur \(\boldsymbol{u_1}\)
    • Le deuxième individu se projette en \(\class{alea}{x_{2}}^\top \boldsymbol{u_1}\) sur \(\boldsymbol{u_1}\)
  • On rappelle que \[{\boldsymbol{X}}= \begin{pmatrix} \class{alea}{x_{1}}^\top \\ \vdots \\ \class{alea}{x_{n}}^\top \end{pmatrix}.\] La projection de l’individu \(i\) est obtenue en multipliant le vecteur ligne \(\class{alea}{x_{1}}^\top\) par \(\boldsymbol{u_1}\). Le vecteur des coordonnées sur \(\boldsymbol{u_1}\) pour chaque individu est

\[\begin{pmatrix} \class{alea}{x_{1}}^\top \boldsymbol{u_1} \\ \vdots \\ \class{alea}{x_{n}}^\top \boldsymbol{u_1} \end{pmatrix} = X \boldsymbol{u_1}.\]

Maximiser \(I_{u_1^\perp}\)

sous la contrainte \(\lVert\boldsymbol{u_1}\rVert\) \[\begin{align} I_{u_1^\perp} & = \frac{1}{n}\lVert X \boldsymbol{u_1}\rVert^2 \\ & = \frac{1}{n} (X \boldsymbol{u_1})^\top (X \boldsymbol{u_1})\\ & = \boldsymbol{u_1}^\top \left (\frac{1}{n} X^\top X \right) \boldsymbol{u_1} \end{align}\]

  • \(\frac{1}{n} X^\top X\) est la matrice de covariance, on écrit souvent \(\frac{1}{n} X^\top X\) sous la forme \(X^\top W X\) (\(W\) est une matrice diagonale, dont chaque terme diagonal vaut \(1/n\), \(W_{ii}\) est le poids attribué à l’observation \(i\)).

  • \(X^\top W X\) est symétrique et définie positive donc diagonalisable.

  • Il existe une matrice \(P\) unitaire \(P^\top P = I_p,\) et \(D\) une matrice diagonale telle que \[X^\top W X = P D P^\top.\]

  • Les termes diagonaux de \(D\) sont les valeurs propres, toutes positives, de \(X^\top W X\) rangées par odre décroissant, \(\lambda_1 > \lambda_1 \ldots >\lambda_p\) et \(U\) est la matrice des vecteurs propres associés.

\(\boldsymbol{b}\) est une base orthonormale, la base formée par les vecteurs propres de \(X^\top W X\).

Identifier les directions principales

Chercher \(\boldsymbol{u_1}\) tel que \(I_{u_1^\perp}\) soit maximale

En se plaçant dans la nouvelle base, posons \(\boldsymbol{v_1} = P^\top \boldsymbol{u_1},\) \(\lVert\boldsymbol{v_1}\rVert=1\):

Maximiser \(\boldsymbol{v_1}^\top D \boldsymbol{v_1} = \sum_{k=1}^p v_{1p}^2 \lambda_p\), sous la contrainte \(\sum_{k=1}^p v_{1p}^2=1\)

\[\boldsymbol{v_1} = (1, 0, \ldots, 0)^\top, \quad \boldsymbol{u_1} = P \boldsymbol{v_1} = \begin{pmatrix} P_{11} \\ \vdots \\P_{p1} \end{pmatrix}.\]

\(\boldsymbol{u_1}\) est le premier vecteur propre de \(X^\top W X\). (la matrice de covariance).

Chercher \(\boldsymbol{u_2}\)

On cherche \(\boldsymbol{u_2}\perp \boldsymbol{u_1}\), tel que \(\boldsymbol{u_2}^\top X^\top W X \boldsymbol{u_2}\) soit maximale.

Posons \(\boldsymbol{v_2} = P^\top \boldsymbol{u_2},\) \(\lVert\boldsymbol{v_2}\rVert=1\): Maximiser \(\boldsymbol{v_2}^\top D \boldsymbol{v_2} = \sum_{k=1}^p v_{2p}^2 \lambda_p\), sous la contrainte \(\sum_{k=1}^p v_{1p}^2=1\) et \(v_{1p}=0\) pour assurer l’orthogonalité avec \(\boldsymbol{v_1}\)

\[\boldsymbol{v_2} = (0, 1, \ldots, 0)^\top, \quad \boldsymbol{u_2} = P \boldsymbol{v_2} = \begin{pmatrix} P_{12} \\ \vdots \\U_{p2} \end{pmatrix}.\]

\(\boldsymbol{u_2}\) est le deuxième vecteur propre de de \(X^\top W X\).

Mesurer l’information préservée par projection

Sur un axe principal

  • Lorsque l’on résume le nuage à sa projection sur un axe factoriel \(\boldsymbol{u_k}\), on ne garde que l’information associée à cet axe soit \(I_{\boldsymbol{u_k}^\perp}=\lambda_k\).

  • La quantité d’information préservée par projection sur l’axe \(k\) est donc \(\lambda_k\).

  • La part d’information portée par l’axe \(k\) est donnée par \(\displaystyle\frac{\lambda_k}{\sum_{k=1}^p\lambda_k}\)

Sur un axe principal

  • L’inertie du nuage projeté sur plan \(k_1, k_2\) est donnée par \[\lambda_{k_1} +\lambda_{k_2}\]

  • La part d’information portée par le plan \(k_1, k_2\) est donnée par \[\frac{\lambda_{k_1} +\lambda_{k_2}}{\sum_{k=1}^p\lambda_k}\]

Attention aux unités - Illustration sur les données de manchots

library(FactoMineR)

penguins.PCA <- PCA(X = penguins, 
                    scale.unit = FALSE, ## should we scale the data ?
                    ncp = 4, ## number of components : as many as initial variables
                    graph = FALSE, ## no graph right now
                    quali.sup = c(1, 2, 7,8)) ## to indicate the column number of qualitative variable
penguins.PCA$eig

penguins.PCA$svd$V ## what is named P in the lecture
         eigenvalue percentage of variance cumulative percentage of variance
comp 1 6.465859e+05           9.998932e+01                          99.98932
comp 2 5.058345e+01           7.822325e-03                          99.99715
comp 3 1.610647e+01           2.490736e-03                          99.99964
comp 4 2.353778e+00           3.639929e-04                         100.00000
             [,1]        [,2]        [,3]          [,4]
[1,]  0.004003162  0.31927773  0.94126475 -0.1098470736
[2,] -0.001154327 -0.08684753  0.14449479  0.9856862728
[3,]  0.015194547  0.94354238 -0.30518986  0.1278908060
[4,]  0.999875876 -0.01571702  0.00103611 -0.0003657482

La dispersion du nuage projeté sur le premier axe représente \(6.45 10^5\) en terme d’inertie soit \(99.98%\) de l’inertie totale.

Trop beau pour être vrai ?

Réfléchissons aux données, ce qui est crucial dans toutes analyses de données.

penguins %>% 
  summarise(across(where(is.numeric), list(mean = mean, sd= sd)))  |> 
  pivot_longer(cols = everything(),  names_to = c("Var", "Resume"), 
      names_pattern = "(.*)_(mean|sd)") |> 
  pivot_wider(names_from = Resume,values_from = value)
# A tibble: 4 × 3
  Var                 mean     sd
  <chr>              <dbl>  <dbl>
1 bill_length_mm      44.0   5.47
2 bill_depth_mm       17.2   1.97
3 flipper_length_mm  201.   14.0 
4 body_mass_g       4207.  805.  

Les variables ont des ordres de grandeur variable.

La contribution à la distance entre individu de la variable Massest plus importante que celle de la variable bill_length_mm.

Est-ce souhaitable ?

Attention aux unités - Illustration sur les données de manchots normalisées

penguins.PCA <- PCA(X = penguins, 
                    scale.unit = TRUE, ## should we scale the data,
                    ncp = 4, ## number of components : as many as initial variables
                    graph = FALSE, ## no graph right now
                    quali.sup = c(1, 2, 7,8)) ## to indicate the column number of qualitative variable
penguins.PCA$eig

penguins.PCA$svd$V ## what is named P in the lecture
       eigenvalue percentage of variance cumulative percentage of variance
comp 1  2.7453557              68.633893                          68.63389
comp 2  0.7781172              19.452929                          88.08682
comp 3  0.3686425               9.216063                          97.30289
comp 4  0.1078846               2.697115                         100.00000
           [,1]       [,2]       [,3]       [,4]
[1,]  0.4537532 0.60019490 -0.6424951 -0.1451695
[2,] -0.3990472 0.79616951  0.4258004  0.1599044
[3,]  0.5768250 0.00578817  0.2360952  0.7819837
[4,]  0.5496747 0.07646366  0.5917374 -0.5846861

On a centré réduit les variables pour les rendre comparables et leur donner la même importance.

Dans ce cas l’inertie totale est le nombre de variables (voir remarque Section 3.1.1)

Lorsque l’on veut s’abstraire des unités, on centre te on réduit les données pour réaliser une ACP normée.

Choix de la distance

L’ACP normée (ACP sur les données centrées et réduites) peut-être envisagée comme une ACP avec une autre distance que la distance euclidienne.

Une autre manière de présenter la même approche est de dire qu’on change la métrique utilisée

Definition (métrique)

Une métrique \(M\) permet de définir une distance entre les observations (ainsi, M est une matrice de \({\mathbb{R}}^p\times{\mathbb{R}}^p\) symétrique, définie et positive).

Calcul de distance

Soient \(\boldsymbol{x}_1 \in {\mathbb{R}}^p\) et \(\boldsymbol{x}_é\in{\mathbb{R}}^p\), la distance entre \(\boldsymbol{x}_1\) et \(\boldsymbol{x}_2\)pour la métrique \(M\) est \[d_M(\boldsymbol{x}_1, \boldsymbol{x}_2)^2= \lVert\boldsymbol{x}_1 - \boldsymbol{x}_2\rVert_M^2 = (\boldsymbol{x}_1 - \boldsymbol{x}_2)^\top M (\boldsymbol{x}_1 - \boldsymbol{x}_2).\]

Si \(M= I_p\),

\[d_M(\boldsymbol{x}_1, \boldsymbol{x}_2)^2= \lVert\boldsymbol{x}_1 - \boldsymbol{x}_2\rVert_M^2 = (\boldsymbol{x}_1 - \boldsymbol{x}_2)^\top I_p (\boldsymbol{x}_1 - \boldsymbol{x}_2) =\sum_{k=1}^p (x_{1k}-x_{2k})^2.\]

Si \(M\) est diagonale est \(M_{ii}= 1/ s_i^2\)\(s_i^2\) est la variance de la variable \(i\),

\[d_M(\boldsymbol{x}_1, \boldsymbol{x}_2)^2= (\boldsymbol{x}_1 - \boldsymbol{x}_2)^\top M (\boldsymbol{x}_1 - \boldsymbol{x}_2) =\sum_{k=1}^p \frac{(x_{1k}-x_{2k})^2}{s_k^2}=\sum_{k=1}^p \left ( \frac{x_{1k}}{s_k}-\frac{x_{2k}}{s_k} \right)^2.\]

Dans une ACP normée on s’intéresse donc aux valeurs propres de la matrice \(X^\top M W X,\) il s’agit de la matrice de corrélation.

Etude du nuage des individus

Représentation des individus - Composantes principales

Definition (Composantes principales)

Les coordonnées des \(n\) observations sur les \(p\) axes factoriels sont groupées dans la matrice \(C\): \(c_{ik}\) est la coordonnée de l’observation \(i\) sur l’axe \(\Delta_{\boldsymbol{u_k}}\). La \(k\)-ème composante principale est la colonne \(k\) de la matrice \(C\), notée \(C^k\). \[ C = [\class{orange}{C^1} | \ldots | \class{orange}{C^p}] \text{ avec } \class{orange}{C^k}=\begin{pmatrix} c_{1k}\\ \vdots\\ c_{nk} \end{pmatrix}. \]

On a les expressions suivantes \[ C^k=Xu_k \mbox{ et } C=XP, \]\(P\) désigne la matrice de changement base, i.e la matrice contenant les vecteurs propres.

  • Les \(C^k\) sont des combinaisons linéaires des variables de départ.
  • Elles sont centrées, de variance \(\lambda_k=\sum_{i=1}^n p_i c_{ik}^2\) et non-corrélées 2 à 2.
  • Les \(C^k\) sont vecteurs propres de \(XX^\top\) de valeur propre \(\lambda_k\).

Représentation des manchots

X <- penguins |> 
  mutate(year = as.factor(year))|>  ## year willnot be considered as numerical variable
  select(where(is.numeric)) |>      ## select all numeric columns
  mutate_all(list(scale2)) |>          ## and scale them
  as.matrix()

n <- nrow(X) # number of individuals
d <- ncol(X)

penguins_svd <- svd( 1/n * t(X)%*% X ) # to get P and D from t(X)WX

## eigenvalues
penguins_eigenvalue <- penguins_svd$d

P <- penguins_svd$v
cp <- X %*% P
head(cp[, 1:2])
cat("\n")
head(penguins.PCA$ind$coord[,1:2])

#plot(penguins.PCA, choix = "ind", axes = c(1,2),invisible=c('ind.sup'),label =c('ind') )
library(factoextra)
fviz_pca_ind(penguins.PCA, label="none", addEllipses=FALSE)

         [,1]        [,2]
[1,] 1.853593 -0.03206938
[2,] 1.316254  0.44352677
[3,] 1.376605  0.16123048
[4,] 1.885288  0.01235124
[5,] 1.919981 -0.81759813
[6,] 1.773020  0.36622296

      Dim.1       Dim.2
1 -1.853593  0.03206938
2 -1.316254 -0.44352677
3 -1.376605 -0.16123048
4 -1.885288 -0.01235124
5 -1.919981  0.81759813
6 -1.773020 -0.36622296

Qualité de la représentation d’un individu

  • Même avec ces précautions, on déforme le nuage de points
  • Certains sont plus déformés que d’autre, il faut les identifier

Un individu est d’autant moins déformé qu’il est proche du plan sur lequel on le projette, ce qui est donné par l’angle entre le vecteur initial et le vteur projeté

\[\begin{align} <\class{alea}{x_i}, \boldsymbol{u_1}> &= \lVert\class{alea}{x_i}\rVert \lVert\boldsymbol{u_1}\rVert cos(\theta_i) \\ & = \class{alea}{x_i}^\top \boldsymbol{u_1}, \quad \left( \text{ou } \class{alea}{x_i}^\top M \boldsymbol{u_1} \text{avec une autre métrique} \right) \\ \end{align}\]

\[\cos^2(\theta_i) = \frac{(\class{alea}{x_i}^\top \boldsymbol{u_1})^2}{\lVert\class{alea}{x_i}\rVert^2}\]

Interprétation: plus \(\cos^2(\theta_i)\) est proche de \(1\), mieux l’individu est représenté.

Contribution d’un individu à la définition de l’axe

  • Les individus sont plus ou moins important dans la définition des axes

Quels sont les points qui vont le plus contribuer à la détermination des axes ?

La contribution relative de l’observation \(\class{alea}{x_i}\) à la construction de l’axe principal \(k\) est la part d’inertie de l’axe \(k\) expliquée par \(\class{alea}{x_i}\) donnée par

\[Ctr(i, k) = \frac{1}{n} \frac{c_{ik}^2}{\lambda_k}.\] En effet \(\lambda_k = \frac{1}{n} \sum_{i=1}^n c_{ik}^2\) (ie. la variance du nuage projeté).

Remarque : si les individus ont des poids \(w_i\) alors \[Ctr(i, k) = \frac{w_i c_{ik}^2}{\lambda_k}.\]

Etude des variables

Qualité de la représentation d’une variable sur un plan factoriel

  • Une variable est bien représentée sur l’axe \(k\) si l’angle entre cette variable et la composante \(C^k\) est petit, son cosinus est proche de 1 ou -1.

Or \(<X^j, C^k> = \lVert X^j\rVert\lVert C^k\rVert cos(\theta_{jk})\)

\[\begin{align} <X^j, C^k> & = <(CP^\top)^j, C^k> \\ & = \underbrace{\left((CP^\top)^j\right)^\top}{= (P C^\top)_j} C^k \\ \end{align}\]

Ainsi \(<X^j, C^k>\) est le \(j^\text{ème}\) terme du vecteur \(PC^\top C^k\), or \(C^\top C\) est diagonale avec \(C_k^k = n lambda_k\), donc \(C^\top C^k=(0, \dots, n\lambda_k, 0, \ldots)\) et

\[\cos(\theta_{jk}) = \frac{n \lambda_k \boldsymbol{u_k^j}}{\sqrt{\lVert X^j\rVert^2 \, \lVert C_k\rVert^2}}= \sqrt{\lambda_k} \boldsymbol{u_k^j}, \quad \cos^2(\theta_{jk}) = \lambda_k{(\boldsymbol{u_k^j})^2}\]

Ajouter des informations supplémentaires

Il est possible d’ajouter des informations supplémentaires

Des individus supplémentaires qui ne seront pas utiliser pour cuitruire les axes

Si on veut ajouter un nouvel individu \(x_{new}\)

  • on calcule \(\boldsymbol{x}_{new} - \boldsymbol{g}\), le point centré (et potentiellement on divise par l’écart type des variables calculés sur le jeu de données initial)

  • On projette sur les axes principaux : \((\boldsymbol{x}_{new} - \boldsymbol{g} )^\top \boldsymbol{u_k}\)

Une variable supplémentaire

Si on souhaite ajouter une nouvelle variable \(\boldsymbol{x}^{new}= (x^{new}_1, \ldots, x^{new}_n)^\top \in {\mathbb{R}}^n\)

  • On soustrait sa moyenne \(\tilde{\boldsymbol{x}}^{new} = \boldsymbol{x}^{new} - m^{new}\), \(m^{new}=\frac{1}{n} \sum_{i=1}^n x^{new}_i\) (ou \(m^{new}= \sum_{i=1}^n w_i x^{new}_i\) si toutes les observations n’ont pas le même poids)

  • On projette \(\tilde{\boldsymbol{x}}^{new}\) sur \(\boldsymbol{u_k},\)

\[Cor(\tilde{\boldsymbol{x}}^{new}, u_k) = \frac{(\tilde{\boldsymbol{x}}^{new})^\top W \boldsymbol{u_k} }{\lVert\tilde{\boldsymbol{x}}^{new}\rVert}\]

Retournons à nos manchots

penguins.PCA <- PCA(X = penguins, 
                    scale.unit = TRUE, ## should we scale the data,
                    ncp = 4, ## number of components : as many as initial variables
                    quali.sup = c(1, 2, 7,8)) ## to indicate the column number of qualitative variable

 # library(Factoshiny)
 # PCAshiny(penguins)

Comprendre les données du Doubs

Présentation des données

L’exemple des caractéristiques du Doubs

On a mesuré les caractéristiques physico chimiques sur 30 sites différents le long de la rivière Doubs.

Les 6 premières lignes (parmi 30) du jeu de données doubs.env

# A tibble: 30 × 11
    das   alt   pen   deb    pH   dur   pho   nit   amm   oxy   dbo
  <dbl> <int> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
1   0.3   934  48    0.84   7.9    45  0.01  0.2   0     12.2   2.7
2   2.2   932   3    1      8      40  0.02  0.2   0.1   10.3   1.9
3  10.2   914   3.7  1.8    8.3    52  0.05  0.22  0.05  10.5   3.5
4  18.5   854   3.2  2.53   8      72  0.1   0.21  0     11     1.3
5  21.5   849   2.3  2.64   8.1    84  0.38  0.52  0.2    8     6.2
6  32.4   846   3.2  2.86   7.9    60  0.2   0.15  0     10.2   5.3
# ℹ 24 more rows
  • das : distance à la source ( \(km\) ),
  • alt : altitude ( \(m\) )
  • pen : la pente (dénivelé pour 1000m)
  • deb : le débit () \(m^3.s^{-1}\) )
  • pH : le pH de l’eau,
  • dur : la concentration en calcium ( \(mg.L^{-1}\) ),
  • pho : concentration en phosphate ( \(mg.L^{-1}\) ),
  • nit : concentration en nitrate ( \(mg.L^{-1}\) ),
  • amn : concentration en ammonium ( \(mg.L^{-1}\) ),
  • oxy : concentration d’oxygène dissous ( \(mg.L^{-1}\) ),
  • dbo : Demande biologique en oxygène ( \(mg.L^{-1}\) ).

Comment visualiser au mieux ces données pour faire apparaître les liens entre variables et identifier des ressemblances entre individus ?

Analyse par ACP

  • Poids des individus ?
  • Quelle distance entre les variables ?
  • Lancer l’analyse
    • quelles sont les composantes intéressantes (règle empririque du choix du nombre d’axes) ?
    • quelles informations sont visibles sur le premier plan factoriel ?
    • Que peut on dire des liens entre les variables ?
    • A votre avis quels sites contribuent le plus à la formation de l’axe 1 ?
    • Que peut on dire du pH ?

L’ACP comme une méthode de compression

Ne garder que l’information pertinente

On a utilisé l’ACP comme une technique de visualisation, mais on peut sen servir comme une technique de compression (avec perte). Il suffit de ne garder que les première composantes principales.

L’exemple sur une image

Point Technique

Matrice d’inertie

Revisiter l’inertie

L’inertie du nuage de points est définie comme \[I = \sum_{i=1}^n w_i \lVert\class{alea}{x_i} -\boldsymbol{g}\rVert^2_M\] On a considéré des variables centrées donc \(\boldsymbol{g} = \boldsymbol{0}\).

  • Or

\[\begin{align} I = & \sum_{i=1}^n w_i x_i^\top M x_i = tr(\sum_{i=1}^n w_i x_i^\top M x_i ) \\ = & tr(\sum_{i=1}^n M x_i w_i x_i^\top ) \quad \text{(commutatitvité de la trace)} \\ = & tr(M \sum_{i=1}^n x_i w_i x_i^\top )\\ = & tr(M X^\top W X)\\ = & tr( X^\top W X M)\\ = & tr( V M)\\ \end{align}\]

Exemple sur les cas présentés

  • Dans le cas de l’ACP, la matrice \(V\) est la matrice de covariance et \(M = Id\),

  • Dans le cas de l’ACP normée , la matrice \(V\) est la matrice de covariance et \(M = Diag(s_k^{-2})\),

Bilan

Ce que vous pensez devoir retenir

  • l’ACP permet de mettre en valeur l’information principale contenu dans un jeu de données

    • les relations entre variables
    • les proximités entre individus
  • Cette méthode permet d’obtenir une projection sur queluqes axes indépendants et bien choisis, au sens ou ils résument le mieux l’information

  • Techniquement l’ACP consiste à diagonaliser la matrice de covariance

  • Les axes sont définis par les vecteurs propres de la matrice
    \[ X^\top W X\]

  • Ils sont rangés par ordre d’importance décroissante, au sens où le premier est celui qui préserve le plus l’inertie du nuage projeté (i.e a la plus grande valeur propre)

  • Il faut adapter la métrique au probleme

  • l’information portée par un axe est mesurée par la valeur propre associé à l’axe

  • il faut savoir interpréter le nuage des individus, éclairés par les liens entre axes principaux et anciennes variables