02-ACP.Rmd

# L'ACP

## Principe de l'ACP

### Type de données

On dispose d'un tableau de $n$ lignes et $p$ colonnes actives qui sont uniquement des variables quantitatives (au sens large, des binaires codées 0/1 ou des variables ordinales passent). 

<center>![](images/DF1.png)</center>

### Objectifs

Une ACP peut permettre de : 

- Résumer l'information.
- Identifier les corrélations entre variables actives.
- Identifier les proximités entre les individus.

### Méthode

L'ACP va déterminer, dans un espace à $p$ dimensions (sous l'hypothèse qu'aucune des variables n'est une combinaison linéaire des autres), l'axe le long duquel les points représentant les individus sont les plus "étalés". Cet axe $F_1$ est une combinaison linéaire des variables de départ.

Dans l'exemple graphique ci-dessous et en 2D, la dispersion des points est essentiellement le long de l'axe bleu. L'axe rouge est orthogonal au bleu. Dans le repère $(O, x, y)$, la variance est à peu près équivalente sur $x$ et sur $y$. A l'inverse, dans le repère constitué des axes de couleur, la variance est portée essentiellement sur le bleu tandis que le rouge ne porte qu'une faible partie de l'inertie du nuage de point. En gros, l'axe bleu est la première composante principale, généralement notée $PC_1$ ou $F_1$, et le rouge $PC_2$ ou $F_2$.

```{r acp_1, echo = FALSE}
library (ggExtra)
set.seed (123)
x <- rnorm (300, 50, 10)
y <- x + rnorm (300, 50, 3)
mod <- lm (y ~ x)
df <- data.frame (x,y)
p <- ggplot(df, aes(x, y)) +
  geom_point () +
  geom_smooth (method = 'lm', formula = y ~ x, se = FALSE, lwd = 1.2, col = 'blue') +
  geom_abline (slope = -1 / mod$coefficients[2], intercept = 155, col = 'red', lwd = 1.2) +
  coord_fixed (ratio = 1)
ggExtra::ggMarginal (p, type = "histogram")
```

Ce changement de repère peut être visualisé ainsi :

```{r acp_2, echo = FALSE}
library (grid)
df <- data.frame (F1 = mod$fitted, F2 = mod$residuals)
p <- ggplot (data = df, aes (x = F1, y = F2)) +
  geom_point () +
  geom_hline (yintercept = 0, col = 'blue', lwd = 1.2) +
  geom_vline (xintercept = mean (mod$fitted, na.rm = T), col = 'red', lwd = 1.2) +
  coord_fixed (ratio = 1) +
  geom_point (aes(x=df [255,1], y=df [255,2]), colour="purple") +
  geom_segment (aes(x = mean (mod$fitted, na.rm = T), y = 0,
                    xend = df [255,1], yend = df [255,2]),
                arrow = arrow(length = unit(0.5, "cm")), color = "purple", lwd = 1) +
  geom_segment (aes(x = df [255,1], y = 0,
                    xend = df [255,1], yend = df [255,2]),
                color = "purple", linetype = "dashed", lwd = 1) +
  geom_segment (aes(x = mean (mod$fitted, na.rm = T), y = df [255,2],
                    xend = df [255,1], yend = df [255,2]),
                color = "purple", linetype = "dashed", lwd = 1) +
  annotate (geom = "text", x = 108, y = 6.5, label = "d1", color = "purple") +
  annotate (geom = "text", x = 117, y = 3, label = "d2", color = "purple") +
  annotate (geom = "text", x = 107, y = 1.5, label = "d", color = "purple", cex = 5)
  
ggExtra::ggMarginal (p, type = "histogram")
```

Dans ce nouveau repère, l'essentiel de la dispersion des points peut être "résumé" par leur position sur $F_1$.

La *variance* du nuage de points caractérise sa dispersion. C'est la moyenne, sur l'ensemble des points, des carrés des distances au centre d'inertie du nuage. D'après le théorème de Pythagore, elle peut être décomposée en composantes orthogonales car $d^2=d_1^2+d_2^2$.

$V={\frac  1n}\sum _{{i=1}}^{n}(d_{i})^{2}\\
  ={\frac  1n}\sum _{{i=1}}^{n}[(d_{1i})^{2}+(d_{2i})^{2}]\\
  ={\frac  1n}\sum _{{i=1}}^{n}(d_{1i})^{2}+{\frac  1n}\sum _{{i=1}}^{n}(d_{2i})^{2}\\
  =V_1+V_2$

Lors du changement de repère, on a conservé la variance totale du nuage de points, mais au lieu d'être répartie également entre les variables $x$ et $y$, elle est concentrée sur $F_1$.

Pour un nombre de variables $p>2$, on généralise. $F_1$ est l'axe qui porte le plus d'inertie, puis $F_2$ est, dans le sous-espace orthogonal à $F_1$, celui qui porte le plus d'inertie, et ainsi de suite. Le nombre total d'axes est $p$ (sauf à avoir des variables qui sont des combinaisons linéaires les unes des autres, ce qui intervient par exemple quand la somme des colonnes fait 100%).

Faire une ACP revient donc à effectuer un changement de repère. Dans notre *dataframe* de départ, chaque colonne peut être vue comme une coordonnée de chacun des individus. Ce n'est pas un repère orthonormé car les variables sont plus ou moins corrélées les unes aux autres. Quand on représente le nuage de points dans le repère $(O, F_1, F_2)$, les axes sont indépendants (corrélation nulle).

En d'autres termes, l'ACP consiste à construire un nouvel espace vectoriel orthonormé (les variables construites sont non-corrélées 2 à 2) de même dimension que l'espace de départ, mais où l'inertie sera concentrée sur les premiers *axes factoriels*. Mathématiquement, cela conduit à diagonaliser la matrice de variance-covariance ; les valeurs propres correspondent à la part de l'inertie totale portée par chaque axe factoriel.

### Centrer - réduire ?

Si l'on a des variables de dispersion (variance) très différentes, ou d'unités différentes, en général il fait centrer (retrancher la moyenne) et réduire (diviser par l'écart-type) chacune des variables avant d'effectuer l'ACP. Dans ce cas chacune des variables a la même importance (variance de 1) $\Rightarrow$ l'inertie du nuage vaut $p$.

C'est l'option par défaut dans la fonction `PCA` de FactoMineR (`scale.unit = TRUE`), qui réalise donc une ACP dite "normée".

**Note :** Les espaces des variables et des individus ne sont pas les mêmes en ACP $\Rightarrow$ on ne peut pas les représenter simultanément. Dans les autres analyses, on le peut.

## L'ACP avec FactoMiner

### Ressources

- Une [explication en vidéo par l'auteur du package FactoMiner](https://www.youtube.com/watch?v=1QPRsg3Bxok)
- Un [tuto associé](http://factominer.free.fr/course/doc/AnaDo_ACP_Facto_Decathlon_Markdown.pdf) en pdf.
- Des [éléments plus généraux sur l'ACP](http://www.math.univ-toulouse.fr/~besse/Wikistat/pdf/st-m-explo-acp).

### Exemple

#### Données utilisées

- Données communales téléchargées sur [Géoïdd](http://geoidd.developpement-durable.gouv.fr/geoclip_stats_o3/index.php?profil=FR#s=2012;v=map1;i=ocs.p_surf_c1;l=fr) et converties en format *.csv* (**attention aux valeurs N/A ou autres dans le fichier Excel qui peuvent parasiter l'importation de données**).
- Les individus statistiques sont les communes (lignes) et les variables (colonnes) différents indicateurs décrivant ces communes.
- Objectif : voir quels indicateurs différencient le plus les communes, et quelles sont les communes qui s'écartent le plus de la moyenne.
- Utilisation de variables qualitatives supplémentaires pour compléter la description (présence d'un agenda21 et d'un PPRN).


#### Importation et exploration rapide

Les données sont dans le fichier *ACP.csv*.

Le fait d'affecter les code géographiques comme `row.names` du dataframe permet de mieux identifier les individus sur les graphiques (et les sorties) par la suite.

ATTENTION, certaines fonctions du package `dplyr` du tidyverse, comme `mutate`, suppriment les noms des lignes.

```{r acp_3}
dat <- geoidd %>%
  select (-code, -communes)
row.names (dat) <- geoidd$code

dat$agd21 <- as.factor (paste ("Agenda21 ", (dat$part_agenda21 > 0) + 0))
dat <- dat %>%
  select (- part_agenda21) %>%
  na.omit()
```

Matrice de corrélations sur les variables quantitatives :

```{r acp_4}
dat %>% select (-PPRN, -agd21) %>% cor() %>% round (digits = 2)
```

Tableau croisé sur les variables qualitatives :

```{r acp_5}
table (dat$agd21, dat$PPRN)
```

On remarque dans un premier temps que les corrélations ne sont pas très élevées entre les variables. La corrélation la plus élevée est celle entre la densité de population et la part de terres artificialisées.
Le lien entre les variables *PPRN* et *Agenda21* est peu évident avec le tableau de fréquences.

On peut utiliser les fonctions vues dans le module 3 pour explorer les données :

```{r acp_6, cache=TRUE}
select (dat, -PPRN, -agd21) %>% 
  ggpairs()
```

#### Réaliser l'ACP

Pour voir l'aide : `?PCA`.

On voit que par défaut les variables sont centrées - réduites et que les valeurs manquantes sont remplacées par la moyenne de la colonne à laquelle elles appartiennent. Il est important de se questionner sur la pertinence de ces options "par défaut".

Réalisation de l'ACP :

```{r acp_7}
acp <- PCA (X = dat, quali.sup = c (5,7), graph = FALSE)
```

La fonction `PCA` retourne un objet de type liste (de liste), qui contient toutes les informations nécessaires à l'interprétation des résultats et leur utilisation : on y retrouve notamment, pour les individus **et** pour les variables, sur chacun des axes :

- Les coordonnées factorielles (*coord*).
- La qualité de représentation (*cos2*) sur chaque axe.
- La contribution à la formation de l'axe (c'est à dire la part de variance de l'axe portée par l'individu / la variable).

Pour accéder aux résultats de l'ACP, on peut voir les objets contenus dans cette liste :

```{r acp_8}
names (acp)
```

Chacun de ces sous-objets (par exemple `acp$ind`) contient à son tour des éléments :

```{r acp_9}
names (acp$ind)
```

Pour tout visualiser d'un coup, on peut utiliser la fonction `str` :

```{r acp_10, eval = FALSE}
str (acp)
```

L'élément `eig` donne les valeurs propres de la matrice de variance-covariance, autrement dit la part d'inertie portée par chacun des axes factoriels. Pour accéder aux coordonnées factorielles des individus, on apppelle l'élément *coord* de l'élément *ind* de l'objet *acp*.

```{r acp_11}
head (acp$ind$coord, n = 10) %>% round (digits = 2) %>% datatable ()
```

Contributions des individus à chacun des axes :

```{r acp_12}
acp$ind$contrib %>% round (digits = 5) %>% datatable ()
```

Recherche des individus qui contribuent le plus aux axes. Pas mal de fonctions du tidyverse font disparaître les noms des lignes $\Rightarrow$ on les stocke provisoirement dans une variable avec la fonction `rownames_to_column()` puis on les ré-attribue avec la fonction réciproque `column_to_rownames()`

```{r acp_13}
acp$ind$contrib %>%
  round (digits = 5) %>%
  as.data.frame () %>% 
  rownames_to_column (var = 'Commune') %>% 
  mutate (somme = rowSums (select(., starts_with ("Dim")))) %>% 
  filter (somme > 0.2) %>% 
  column_to_rownames (var = 'Commune') %>% 
  datatable ()
```

En ordonnant le tableau selon les différentes composantes principales, on peut identifier les communes qui contribuent le plus à chacun des axes.

On peut également utiliser les fonctions génériques `summary` et `plot` (ou `plot.PCA` ; cf. `?plot.PCA`) pour aborder ce nouvel objet.

```{r acp_14}
plot (acp)
```

#### Nombre d'axes à retenir

On accède aux valeurs propres par l'objet `acp$eig` :

```{r acp_15}
acp$eig
```

On peut représenter graphiquement les valeurs propres (« éboulis » des valeurs propres).

```{r acp_16}
eig <- as.data.frame (acp$eig)
mm <- mean (eig$`percentage of variance`) / 100

ggplot (eig, aes (x = 1:nrow(eig), weight = `percentage of variance`)) +
  geom_bar (fill = "grey20") + 
  ggtitle ("Valeurs propres") + 
  theme (axis.title.x = element_blank(), axis.title.y = element_blank()) +
  scale_y_continuous (labels = function(x) paste0(x, "%")) +
  geom_hline (yintercept = 20, colour="red", linetype="dashed") +
  coord_flip ()
```

Pour choisir le nombre d'axes à conserver, on peut utiliser plusieurs critères :

- Critère de l'inertie moyenne : on retient les axes qui représentent plus d'inertie que la moyenne (=inertie totale/nombre de variables). Cette inertie vaut **1** dans le cas de l'ACP normée (soit 20% dans notre cas). Cela revient à se dire qu'on ne prend que les nouvelles variables qui portent plus d'inertie que les variables initiales.
- Critère du coude : on retient les premiers axes jusqu'à oberver un "décrochage" dans l'éboulis des valeurs propres. On peut l'objectiver par le calcul des différences secondes entre valeurs propres.

```{r acp_17}
acp$eig %>% diff () %>% diff ()
```

Ici, le critère de l'inertie moyenne incite à retenir les 2 premiers axes alors que le critère du coude en retiendrait 3 (la dérivée seconde change de signe entre la 3e et 4e valeur propre). Généralement, on prend le critère le plus parcimonieux pour s'épargner du travail ensuite !

#### Interpréter les axes

- Rappel : les axes factoriels sont des combinaisons linéaires des variables initiales. On regarde donc l'importance de ces variables initiales dans les différents axes pour leur *donner du sens*.
- On doit regarder trois grandeurs pour interpréter (dans cet ordre) : 
    + Contribution
    + Qualité de représentation
    + Coordonnée

On visualise le **cercle des corrélations**. En ACP, qualité et contribution dépendent directement de la coordonnée : il suffit que celle-ci soit élevée pour que la qualité de représentation soit bonne et la contribution forte. 
Pour obtenir ce graphique, on utilise la fonction `plot.PCA`, avec l'argument `choix = "var"`.

```{r acp_18}
plot.PCA (acp, choix = "var", col.var = "blue")
```

$\Rightarrow$ on regarde avant tout les variables proches du cercle unité.

A partir de ce cercle (sur le premier plan factoriel), on voit que :

- Les variables *part_artif* et *Densite* sont corrélées positivement, et toutes deux négativement à *part_proprio*. 
- La variable *ind_vieill* est indépendante (~orthogonale) aux trois variables précédentes.
- L'axe 1 porte (~synthétise) 43,1% de l'inertie totale. Il est principalement formé par les variables *part_artif*, *Densite* et *part_proprio*.
- L'axe 2 porte 21,8% de l'inertie. Il est formé par les variables *ind_vieill* et, dans une moindre mesure, *tx_emploi*.

On peut interpréter de façon plus fine grâce aux chiffres donnés par la fonction `summary`.

```{r acp_19 summary}
summary (acp)
```

La fonction `dimdesc` donne les variables les plus liées à chacun des axes.

```{r acp_20}
dimdesc (acp)
```

On peut également regarder la projection des individus dans le plan. Ce n'est pas toujours très informatif, ça peut permettre de repérer :
- Des groupes d'individus quand il existe des discontinuités.
- Des valeurs extrêmes qui peuvent contribuer exagérément aux axes (donc qui nuisent à la représentation du reste des individus).

```{r acp_21, cache=TRUE}
plot.PCA (acp, choix = "ind", col.ind = "lightblue")
```

N'oubliez pas de consulter l'aide des fonctions présentées pour plus d'options.

#### Les variables supplémentaires

Dans toutes les analyses factorielles, on distingue les **variables actives**, qui participent de la création du nouvel espace vectoriel, des **variables supplémentaires**, que l'on projette *a posteriori* sur cet espace, mais qui n'interviennent pas dans la construction des variables synthétiques. 

- Elles peuvent être qualitatives ou quantitatives.
- Elles sont projetée sur l'espace des variables (variable continue) ou des individus (qualitative). 
- Elles ne rentrent pas dans la matrice qui sert à la définition des axes. 
- Elles sont utiles pour voir la corrélation d'une variable (une classe d'individus par exemple) avec toutes les variables simultanément

Dans la fonction PCA, on désigne les variables supplémentaires avec les arguments `quali.sup` ou `quanti.sup` ; attention, on les désigne par le numéro de colonne de la variable.

$\Rightarrow$ Pas lien entre le PPRN et les autres variables : elles se projettent au centre du nuage de points. Les communes disposant d'un PPRN ou d'un Agenda21 n'ont donc pas de caractéristiques particulières.

On peut utiliser le paramètre `select` pour mettre en évidence certains individus ou variables en fonction de leur contribution ou qualité de représentation. Par exemple `select = "contrib 20"` permet de n'afficher les étiquettes que des individus qui contribuent le plus aux deux axes représent

```{r acp_22, cache=TRUE}
plot.PCA (acp, choix = "ind", select = "contrib 20")
plot.PCA (acp, choix = "ind", select = "cos2 20")
```

### Repérer un effet taille

L'effet taille se recontre assez fréquemment quand on réalise une ACP : il se manifeste par :

- Toutes les variables sont de même signe sur le premier axe factoriel (donc elles sont toutes corrélées positivement entre elles) et celui-ci contient une très grand part de l'inertie.
- Dans ce cas, l'axe 1 constitue un **gradient** : il permet de classer les individus du plus "petit" au plus "grand", sur toutes les variables simultanément.

```{r acp_23}
taille <- read.csv2 ("data/Effet_taille.csv", header = TRUE,
                     encoding = "latin1")
t <- PCA (na.omit (taille[,-c(1,2)]), graph = FALSE)
plot.PCA (t, choix = c ("var"), col.var = "blue")
```

Ce type de configuration peut décevoir quand on ne s'y attend pas. Cependant, il peut y avoir un intérêt car l'axe $F_1$ synthétise, de manière objective et multidimensionnelle, le lien entre les variables. Ici, revenus, diplômes et statut de cadre sont très liés. Si l'on voulait étudier si ces variables sont corrélées à d'autres caractéristiques des individus, plutôt que de faire 3 analyses semblables, ou de choisir une de ces 3 variables, on pourrait prendre les coordonnées des individus sur $F_1$ (`t$ind$coord[, 1`) comme variable de synthèse.

## Exercice

A partir du `dataframe` *iris* inclus dans R :

- Explorer rapidement le jeu de données. Quel est le type des variables ?
- Réalisez une ACP sur ce jeu de données.
- Interprétez les résultats obtenus.
- Inclure la variable "Species" dans l'analyse en tant que variable qualitative supplémentaire.

```{r exercice02acp, collapse = T, cache=TRUE}
data ("iris")
ggpairs (iris)
acp.iris <- select (iris, -Species) %>% PCA()
summary (acp.iris)
dimdesc (acp.iris)
acp.iris <- PCA(iris, quali.sup = 5)
```