---
title: "Résumé minimal du script AFC_profor.R pour les AFC et CAH sur la propriété forestière en France au début des années 2000, avec des données brutes"
author: "Vincent GODARD"
date: "02/01/2022 - V2"
output: 
  html_document: 
    toc: yes
    keep_md: yes
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(message = TRUE)
```

Analyse factorielle des correspondances (AFC) avec R et FactoMineR sur les données relatives à la propriété forestière  (sources : https://agreste.agriculture.gouv.fr/agreste-web/download/publication/publie/ChdAgr196/ChdAgr196.pdf)

Pour aller plus loin :

- http://factominer.free.fr/factomethods/analyse-factorielle-des-correspondances.html

- http://www.sthda.com/french/articles/38-methodes-des-composantes-principales-dans-r-guide-pratique/74-afc-analyse-factorielle-des-correspondances-avec-r-l-essentiel/#r-code-to-comput-ca

## Résumé des principales étapes du code R

- Calcul du khi2 ;
- Module AFC ;
- Valeurs propres et vecteurs propres ;
- Coordonnées ;
- Contributions ;
- Qualités (COS2) ;
- Aide à l’interprétation ;
- Module CAH ;
- Partitions ;
- Description des classes ;
- Edition des rapports.


## 1. Import des données

Avec mes fichiers exemples, préférez les paramètres nord-américains ("." séparateur de décimales, "," séparateur de miliers).

Salle C210, identifiant => "geographie2" !

```{r}
setwd("D:/3VG/UP8/enseignement/add/proforest/R/AFC")
### setwd("D:/Users/geographie2/vgodard/ADD/proforest/R/AFC")
### respecter ce cheminement si "geographie2"
### getwd()

# Remove all objects
rm(list = ls() )
```

## 2. Lecture des données

```{r}
propforafc <- read.table("PropFor19AFC.csv",
                      header=TRUE,
                      sep=";",
                      dec=".",### Si param. régio. non US, remplacer le "." par une ","
                      row.names=1,
                      check.names=FALSE,
                      fileEncoding="latin1",
                      stringsAsFactors = TRUE) # si stringsAsFactors = FALSE on perd les modalités
                      

```

## 3. Khi-deux


Pour une (petite) table de contingence, il est possible de tester la significativité du lien de dépendance à l'aide d'un test du khi-deux.


```{r message = FALSE}

# Calcul de la statistique du Chi-square (Khi-deux) sur les données numériques
chisq <- chisq.test(propforafc[,1:3])
chisq


```

Au risque 5% (communément admis), on rejette l’hypothèse nulle d’indépendance entre les lignes et les colonnes. 

Sur les grosses tables, il est pratiquement toujours significatif (p-value < 0.05)



## 4. Module AFC

Installer et exécuter préalablement les librairies "FactoMineR", "factoextra", puis les charger

```{r message = FALSE}
### Avec la librairie FactoMineR et factoextra
## install.packages(c("FactoMineR", "factoextra")) ## si pas déjà installés !

library("FactoMineR")
library("factoextra")


```


L'AFC avec tous les éléments actifs (90 lignes et 3 variables numériques et 2 qualitatives)

*-* CA est le nom de l'AFC dans Factominer 

*-* CA fonctionne a minima avec :

+ x : le data frame (la table de contingence) ;

+ ncp : nb de dimensions gardé dans le résultat final ;

+ graph : le graphe par défaut (pour ne pas afficher le graphe graph = FALSE) ;

*-* Si le nom du département est en première colonne :

+ le tableau "utile" va de la deuxième colonne à la 5ème [,2:6] ;

+ il n'y a pas de variables quantitatives supplémentaires ici ;

+ les variables qualitatives (supplémentaires) sont les dernières (quali.sup=c(5:6)) et doivent être courtes pour rendre lisible les graphiques et tableaux.


```{r}
res.propforafc.ca <- CA(propforafc[,1:5], # indiquer ttes les colonnes utiles
                         col.sup = NULL,
                         row.sup = NULL,
                         quali.sup= 4:5, # à numéroter dans le [1:5] cf. supra
                         axes = 1:2,
                         ncp = 3,
                         graph = TRUE)

```

Le graphique avec la fonction Biplot (cf. infra) apporte un net gain de lisibilité.



## 5. Les sorties de la fonction CA

Si ce n'est déjà fait, installer le package "factoextra", puis le charger.

```{r message = FALSE}
## install.packages("factoextra")

library("factoextra")
```

Si on veut connaître tous les résultats possibles de la fonction CA dans Factominer et/ou factoextra :

```{r}

### Listage des résultats

print(res.propforafc.ca)
```


### 5.1. Les *eigenvalues* ou valeurs propres


Contrairement à l'ACP, une *eigenvalue* n'est jamais supérieure à 1 en AFC.

```{r}

eig.val <- get_eigenvalue(res.propforafc.ca)
eig.val


```

Il est aussi possible d'apprécier les ruptures dans la succession des *eigenvalues*. Sur cet exemple à "deux" dimensions, on pourrait s'en passer.

```{r}

 # Visualisation des eigenvalues
fviz_eig(res.propforafc.ca, addlabels = TRUE, ylim = c(0, 80)) # calibrer le ylim avec la lecture du tableau précédent


```

Il n'y a pas de règles pour choisir le nombre d'axes. Eventuellement, retenir les axes qui ont une *eigenvalues* d'une valeur supérieure à :

- 1/(nrow(propforafc)-1 = 1/90 soit environ 1.4%)
ou

- 1/(ncol(propforafc[, 2:4])-1 = 1/2 soit environ 50%)
Sources : http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/113-ca-correspondence-analysis-in-r-essentials/#r-code-to-compute-ca (According to M. T. Bendixen 1995)



```{r}
 # Visualisation des eigenvalues
sp <- fviz_eig(res.propforafc.ca,addlabels = TRUE, ylim = c(0, 80)) # calibrer le ylim avec la lecture du tableau précédent
         
sp + geom_hline (yintercept = 50, linetype = 2, color = "red") # tracer le yintercept avec la valeur de la *eigenvalues* moyenne 
        

```

Liste des paramètres les plus utiles :

```{r}

get_ca_row(res.propforafc.ca)

```


### 5.2. Coordonnées des variables et des individus sur les axes

Interprétation des coordonnées des variables sur les axes factoriels (les colonnes : col).

```{r}

res.propforafc.ca$col$coord ### [toutes les lignes ; colonnes 1 à 3, cf résultats *eigenvalue*]

```

Valeur des coordonnées des variables qualitatives supplémentaires (les categorical var.).

Ici, coordonnées de chaque modalité de chaque variable sur les Dim 1 et 2 !

```{r}

res.propforafc.ca$quali.sup$coord  ## coord. for the supplementary categorical var.

```


Interprétation des coordonnées des individus sur les axes factoriels

```{r}

res.propforafc.ca$row$coord ### [ligne 1 à 90 ; colonnes 1 à 3, mais 2 dimensions !]

```


### 5.3. Contribution des variables et des individus sur les axes

Interprétation des contributions des variables sur les axes factoriels.

Rappel :

+ La somme en colonnes est égale à 100 ;

+ Une variable est considérée comme contributive si elle est supérieure à 100 / nb de variables (ici 30)

+ Un individu est considéré comme contributif s'il est supérieur à 100 / nb d'individus (ici 90)



```{r}
res.propforafc.ca$col$contrib ### [toutes les lignes ; colonnes toutes, cf résultats *eigenvalue*]

```

Interprétation des contributions des individus sur les axes factoriels.

```{r}
res.propforafc.ca$row$contrib ### [ligne toutes ; colonnes toutes] faire un copier-coller pour rechercher les plus contributifs dans un tableur.

```


### 5.4. Qualité de la représentation (COS2) des variables et des individus sur les axes

Interprétation des COS2 des variables sur les axes factoriels.

```{r}
res.propforafc.ca$col$cos2 ### [toutes les lignes ; colonnes toutes, cf résultats *eigenvalue*] pour rechercher les mieux représentées dans un tableur.

```

Interprétation du degré d’association entre les colonnes et un axe particulier des variables qualitatives supplémentaires.

Rappel : La somme en lignes est égale à 1 !


```{r}
res.propforafc.ca$quali.sup$cos2 ##  "cos2 for supplementary categorical var." [toutes les lignes ; colonnes toutes, cf résultats *eigenvalue*] en faire un copier-coller pour rechercher les mieux représentées dans un tableur.
```

Interprétation des qualités de représentation des individus sur les axes factoriels.

```{r}
res.propforafc.ca$row$cos2 ### [ligne toutes ; colonnes toutes] en faire un copier-coller pour rechercher les mieux représentés dans un tableur.
```


## 6. Aide à l'interprétation de l'AFC

### 6.1. Graphiques de la contribution


-  des variables

```{r}

          fviz_ca_col(res.propforafc.ca,
                  col.col = "contrib",
                  axes = 1:2,
                  gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))
                  


```

- des individus

```{r}
          fviz_ca_row(res.propforafc.ca,
                      col.row = "contrib",
                      axes = 1:2,
                      gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))
```


### 6.2. Graphiques de la qualité


- des variables (ici pas très intéressant car 100% sur les deux seuls axes)

-- il semblerait que dans "fviz_ca", il ne soit pas possible de représenter les "quali.sup"


```{r}

          fviz_ca_col(res.propforafc.ca,
                  col.col = "cos2",
                  axes = 1:2,
                  gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))
                  


```

- des filtrages sont possibles

-- en fonction du Cos2  (ici pas très intéressant car 100% sur les deux seuls axes)

```{r}
          ## Lignes avec un cos2 > 0.8
          fviz_ca_row(res.propforafc.ca,
                      select.row = list(cos2 = 0.8),
                      col.row = "cos2",
                      axes = 1:2,
                      gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))
```

-- en fonction des contributions (intéressant pour les individus les plus contributifs)

```{r}
          ## Top 5 des contributions lignes et colonnes
          fviz_ca_biplot(res.propforafc.ca,
                      select.row = list(contrib = 5), 
                      #select.col = list(contrib = 5), # si on veut toutes les variables
                      arrow = c (TRUE, TRUE), # pour mieux voir la proximité avec les axes
                      axes = 1:2,
                      repel = TRUE, # pour éviter les chevauchements dans biplot
                      )
```

Il est possible de faire des Biplots asymétriques, où les coordonnées sont pondérées par le poids de la valeur propre (en fait la racine carrée de la valeur propre *eigenvalue*).
- Voir : Biplot asymétrique in : http://www.sthda.com/french/articles/38-methodes-des-composantes-principales-dans-r-guide-pratique/74-afc-analyse-factorielle-des-correspondances-avec-r-l-essentiel/#r-code-to-comput-ca


### 6.3. Aide à la description et à l'interprétation


Identification des variables les plus remarquables (moins efficace que dans PCA car **pas significatives** en l'absence de test) par composante. Simple tri ds coordonnées.

```{r}
### fonction dimdesc() [in FactoMineR]

res.desc <- dimdesc(res.propforafc.ca, axes = c(1,2))


```
- Aide sur l'axe 1

-- Pour les lignes (pour nous, les individus)

```{r}
### Description de la dimension 1 pour 4 lignes
head(res.desc[[1]]$row, 4)


```
-- Pour les colonnes (pour nous, les variables)

```{r}
### Description de la dimension 1 pour 4 colonnes

head(res.desc[[1]]$col, 4)

```

- Aide sur l'axe 2

-- Pour les lignes (pour nous, les individus)

```{r}
### Description de la dimension 2 pour toutes les lignes
res.desc[[2]]$row


```
-- Pour les colonnes (pour nous, les variables)

```{r}
### Description de la dimension 2 pour toutes les colonnes

res.desc[[2]]$col

```

### 6.4. Export des résultats de l'AFC

- Exporter des pdf ou des png

-- En premier, création des graphes

```{r message = FALSE}
# Scree plot
scree.plot <- fviz_eig(res.propforafc.ca)

# Biplot of row and column variables
biplot.ca <- fviz_ca_biplot(res.propforafc.ca)


```

-- En second, export des graphes

```{r message = FALSE}

library(ggpubr)

## Export d'un pdf (un graphe par page)
ggexport(plotlist = list(scree.plot, biplot.ca), 
         filename = "res.propforafc.ca.pdf")

## Export d'un png par graphe
ggexport(plotlist = list(scree.plot, biplot.ca), 
         filename = "res.propforafc.ca.png")


```
- Exporter les résultats en txt ou csv

A l'aide de write.table() [package FactoMineR] ou write.infile.

```{r message = FALSE}
# Export into a TXT file de la qualité de représentation des variables
write.table(res.propforafc.ca$col$cos2, "res.propforafc.ca$col$cos2_v2.txt", sep = "\t")

# Export into a CSV file de la qualité de représentation des variables
write.table(res.propforafc.ca$col$cos2, file="res.propforafc.ca$col$cos2_v2.csv", col.names=TRUE, sep = ";")

# Export into a CSV file de la qualité de représentation des variables qualitatives suplémentaires
write.table(res.propforafc.ca$quali.sup$cos2, file="res.propforafc.ca$quali.sup$cos2_v2.csv", col.names=TRUE, sep = ";")

# Export into a CSV file des contributions des individus
write.table(res.propforafc.ca$row$contrib, file="res.propforafc.ca$row$contrib_v2.csv", col.names=TRUE, sep = ";")


```           

## 7. Classification ascendante hiérarchique (CAH)

Si ce n'est déjà fait, installer le package "FactoMineR", puis le charger.


```{r message = FALSE}

library(FactoMineR)


```

### 7.1. Détermination du nombre de classes de la CAH

Ne se calcule que sur les variables numériques actives

```{r message = FALSE}

propforafc.num <- subset(propforafc[,1:3]) ## (extraction dans le tableau de départ) 

head(propforafc.num)


```

Si ce n'est déjà fait, installer le package "factoextra", puis le charger.


```{r message = FALSE}

library("factoextra")


```

Essai avec la méthode "silhouette" pour déterminer le nombre optimal de clusters

(On lira avec intérêt à ce sujet : https://www.datanovia.com/en/lessons/determining-the-optimal-number-of-clusters-3-must-know-methods/)


```{r message = FALSE}

fviz_nbclust(propforafc.num, kmeans, method = "silhouette") ## tester les autres méthodes pour voir...


```

### 7.2. Partition de la CAH

Utilisation des sauts d'inertie calculés à la section 7.1 pour indiquer le nombre de classes retenu.

Avec FactoMiner, c'est "HCPC" la Classification ascendante hiérarchique (CAH)


```{r message = FALSE}

res.propforafc.ca.hcpc <- HCPC(res.propforafc.ca, graph = TRUE, nb.clust = 2) ## nb.clust = 2 en fonction du résultat du 7.1


```

### 7.3. Exploitation des résultats de la CAH

Affichage du nom des 5 colonnes résultats d’HCPC


```{r message = FALSE}

names(res.propforafc.ca.hcpc) # Nommage des 5 colonnes résultats d’HCPC


```

Résumé des 5 colonnes résultats d’HCPC


```{r message = FALSE}

summary(res.propforafc.ca.hcpc)


```

Listage des résultats


```{r message = FALSE}

#print(res.propforafc.ca.hcpc) # Intéressant à exporter pas à afficher (sauf contrôle)


```

Affectation des départements à chacune des classes (pour la cartographier les partitions par exemple)

```{r message = FALSE}

res.propforafc.ca.hcpc$data.clust # affectation individus => cluster


```

Description des sur et sous représentations des variables dans les classes (pour établir la légende de la cartographie des résulats par exemple)

+ variables qantitatives

+ variables qualitatives supplémentaires


```{r message = FALSE}

res.propforafc.ca.hcpc$desc.var # description des classes pour qualifier les postes de légende


```


### 7.4. Export des résultats de la CAH

A l'aide de write.table() [package FactoMineR] ou write.infile (le seul des deux qui exporte "$desc.var") !

```{r message = FALSE}

# Export into a CSV file de l'affectation des individus aux clusters
write.table(res.propforafc.ca.hcpc$data.clust, file="res.propforafc.ca.hcpc$data.clust_v2.csv", col.names=TRUE, sep = ";")

# Export into a CSV file de la description des variables quantitatives et qualitatives suplémentaires
write.infile(res.propforafc.ca.hcpc$desc.var, file="res.propforafc.ca.hcpc$desc.var_v2.csv", sep = ";")


```     


La lecture des caractéristiques des classes se présente comme suit :


- Pour les variables quantitatives [on parle de fréquence en AFC (frequency)]

-- glob% => % de la fréquence dans l'échantillon - GLOBAL : représente la part (en pourcentage) des modes de gestion (ForDom, ForColec, ForPriv) dans l'échantillon global (càd la répartition des modes de gestion pour toute la France) ;

-- intern% => % de la fréquence dans la classe - FRE / CLA : représente la part (en pourcentage) de chaque modes de gestion (ForDom, ForColec, ForPriv) dans la classe.

On crée une variable "écart au global" (dans une feuille de calcul par exemple) pour connaître la sur ou sous représentation de cette fréquence dans la classe par rapport à la répartition France entière.  


- Pour les variables qualitatives suplémentaires [on parle de catégorie ou modalité (category)]:

-- Global => % de la modalité dans l'échantillon - GLOBAL : représente la part (en pourcentage) des départements de la régions par rapport à l'ensemble des départements.

-- Mod/Cla => représente le pourcentage des départements de la région dans cette classe possédant la modalité.

-- Cla/Mod => représente le pourcentage des départements de la région présentant la modalité et appartenant à la classe. 

On crée une variable "écart au global" (dans une feuille de calcul par exemple) pour connaître la sur ou sous représentation de cette modalité dans la classe par rapport à la répartition France entière. 



## 8. Le couteau suisse

Le package FactoInvestigate décrit et interprète automatiquement les résultats de votre analyse factorielle (ACP, AFC ou ACM) en choisissant les graphes les plus appropriés pour le rapport (sources :  http://factominer.free.fr/reporting/index_fr.html, mais ce package semble inaccessible à la date de rédaction !)

```{r message = FALSE}
#install.packages("FactoInvestigate")

library(FactoInvestigate)

Investigate(res.propforafc.ca)
```




