Guide d’utilisation de R pour la chimie des otolithes

Par William Fortin, Lola Cousseau

Avant toute chose, il faut preparer l’espace de travail.

rm(list=ls())

library("FactoMineR"); library("factoextra"); library("corrplot"); library("tidyr"); library("ggplot2"); library("ggpubr"); library("ggthemes"); library("RColorBrewer"); library("dplyr"); library("vegan");library("corrplot"); library("dendextend"); library("randomForest"); library("e1071"); library("nlme")
  • rm(list=ls()) : efface tous les objets de l’espace de travail
  • library(“bibliothèque”) : charge les bibliothèques qui seront utilisées tout au long du guide

1. Analyse de composante principale

L’ACP est une méthode de la famille de l’analyse des données et plus généralement de la statistique multivariée, qui consiste à transformer des variables liées entre elles (dites « corrélées » en statistique) en nouvelles variables décorrélées les unes des autres. Ces nouvelles variables sont nommées « composantes principales », ou axes principaux.

En effet, si notre jeu de données possède plus de trois variables (éléments dans notre cas), il pourrait être très difficile de visualiser les données dans une “hyper-espace” multidimensionnelle. L’ACP permet donc de réduire les dimensions d’une donnée multivariée à deux ou trois composantes principales, qui peuvent être visualisées graphiquement, en perdant le moins d’information possible.

Avant de pouvoir effectuer le test, il faut ouvrir le fichier contenant les données.

dataPCA <- read.csv("/Users/WillzyX/Desktop/Maitrise/Atlas_FSL_2019_log_R.csv", header=T, sep=',',dec=".", stringsAsFactors = T)[,-1] %>%
  drop_na(SrCa, BaCa, MgCa, NaCa)
options(max.print=60)

On peut ensuite préparer les données et les variables.

x<-subset(dataPCA, select = c(SrCa, BaCa, MgCa, NaCa))
xnum<-subset(dataPCA,select = c(SrCa, BaCa, MgCa, NaCa))

res.pca <- PCA(xnum)

var <- get_pca_var(res.pca)
ind <- get_pca_ind(res.pca)

La suite d’opérations est facultative, elle permet d’obtenir plus de détails sur l’ACP.

print(res.pca) # Liste contenant les détails de ce qui a été effectué pendant l'ACP
## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 923 individuals, described by 4 variables
## *The results are available in the following objects:
## 
##    name               description                          
## 1  "$eig"             "eigenvalues"                        
## 2  "$var"             "results for the variables"          
## 3  "$var$coord"       "coord. for the variables"           
## 4  "$var$cor"         "correlations variables - dimensions"
## 5  "$var$cos2"        "cos2 for the variables"             
## 6  "$var$contrib"     "contributions of the variables"     
## 7  "$ind"             "results for the individuals"        
## 8  "$ind$coord"       "coord. for the individuals"         
## 9  "$ind$cos2"        "cos2 for the individuals"           
## 10 "$ind$contrib"     "contributions of the individuals"   
## 11 "$call"            "summary statistics"                 
## 12 "$call$centre"     "mean of the variables"              
## 13 "$call$ecart.type" "standard error of the variables"    
## 14 "$call$row.w"      "weights for the individuals"        
## 15 "$call$col.w"      "weights for the variables"
eig.val <- get_eigenvalue(res.pca)
eig.val # % de la variance expliquée par chaque dimension
##       eigenvalue variance.percent cumulative.variance.percent
## Dim.1  1.9673606         49.18401                    49.18401
## Dim.2  0.9289273         23.22318                    72.40720
## Dim.3  0.6986069         17.46517                    89.87237
## Dim.4  0.4051053         10.12763                   100.00000
fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 50)) # Graph du % de la variance expliquée par chaque dimension

fviz_contrib(res.pca, choice = "var", axes = 1, top = 10) # Graph du % de la variance expliquée par la PC1

fviz_contrib(res.pca, choice = "var", axes = 2, top = 10) # Graph du % de la variance expliquée par la PC2

head(var$cos2) # Qualité de représentation des variables sur le graph de l’ACP à l'aide du cosinus carré
##          Dim.1        Dim.2      Dim.3      Dim.4
## SrCa 0.5989794 0.0034595607 0.27920838 0.11835266
## BaCa 0.7424321 0.0005253381 0.01200325 0.24503930
## MgCa 0.2960269 0.5187214126 0.16023048 0.02502121
## NaCa 0.3299222 0.4062209464 0.24716475 0.01669210
fviz_cos2(res.pca, choice = "var", axes = 1:2) # Graph du cos2 des variables sur les Dim.1 et Dim.2

head(var$contrib) # Contributions des éléments à chacune des composantes principales en %
##         Dim.1      Dim.2     Dim.3     Dim.4
## SrCa 30.44584  0.3724254 39.966453 29.215285
## BaCa 37.73747  0.0565532  1.718169 60.487808
## MgCa 15.04691 55.8409077 22.935716  6.176471
## NaCa 16.76979 43.7301137 35.379663  4.120436
res.desc <- dimdesc(res.pca, axes = c(1,2), proba = 0.05)
res.desc$Dim.1 # Donne la corrélation des variables à la PC1
## $quanti
##      correlation       p.value
## BaCa   0.8616450 1.584300e-273
## SrCa   0.7739376 6.150812e-185
## MgCa   0.5440835  3.037277e-72
## NaCa  -0.5743885  3.892287e-82
## 
## attr(,"class")
## [1] "condes" "list"
res.desc$Dim.2 # Donne la corrélation des variables à la PC2
## $quanti
##      correlation       p.value
## MgCa   0.7202232 2.020209e-148
## NaCa   0.6373546 2.338192e-106
## 
## attr(,"class")
## [1] "condes" "list"
head(ind$coord) # Coordonnéess des individus
##      Dim.1      Dim.2     Dim.3       Dim.4
## 1 2.540642 -0.9172929 -2.949217  0.08801606
## 2 1.488525  0.2601334 -1.470317 -0.03599542
## 3 1.195112 -0.4185731 -1.848146 -0.19801623
## 4 1.671327  0.4554859 -1.735733  0.43790007
## 5 1.772030 -1.2086361 -2.283531  0.42419899
## 6 2.262230  0.8686664 -1.893387  0.02157735
head(ind$cos2) # Qualité des individus
##       Dim.1      Dim.2     Dim.3        Dim.4
## 1 0.4033807 0.05258284 0.5435524 0.0004841187
## 2 0.4983033 0.01521856 0.4861868 0.0002913908
## 3 0.2823632 0.03463647 0.6752487 0.0077516243
## 4 0.4501510 0.03343376 0.4855133 0.0309019065
## 5 0.3141552 0.14614810 0.5216938 0.0180028526
## 6 0.5411162 0.07978532 0.3790493 0.0000492281
head(ind$contrib) # Contributions des individus
##        Dim.1       Dim.2     Dim.3        Dim.4
## 1 0.35546873 0.098136972 1.3488971 0.0020718307
## 2 0.12201877 0.007892395 0.3352642 0.0003465172
## 3 0.07865589 0.020434270 0.5297099 0.0104865343
## 4 0.15382867 0.024197277 0.4672307 0.0512838260
## 5 0.17292450 0.170375728 0.8086845 0.0481248806
## 6 0.28183047 0.088008107 0.5559609 0.0001245164

On peut ensuite procéder à la création de la figure de l’ACP. Il peut être préférable de passer à la section 2. Regroupement hiérarchique avant d’effectuer cette figure si vous souhaitez séparer les stations par cluster comme dans l’exemple.

ind.p <- fviz_pca_biplot(res.pca, geom.ind = "point",
                pointshape = 21,
                pointsize = 2.5,
                fill.ind = dataPCA$Clusters,
                col.ind = "black",
                palette = c("#d73027","#fc8d59","#fee090","#ffffbf","#e0f3f8","#91bfdb","#4575b4"),
                arrowsize = 0.8,
                col.var = "black",
                addEllipses = FALSE,
                axes.linetype = "blank",
                mean.point=FALSE)
ggpubr::ggpar(ind.p,
              title = "",
              caption = "Source: factoextra",
              xlab = "PC1 (49.2%)", ylab = "PC2 (23.2%)",
              col.axis="black",
              col.lab="black",
              cex.axis = 1.5,
              cex.lab=1.5,
              legend.title = "Cluster", legend.position = "top",
              ggtheme = theme_classic(base_size = 12, base_family = "Times"))

La dernière étape cosiste à créer un .csv avec les résultats de l’ACP. Vous pouvez copier ces données et les coller à la suite de vos résultats dans votre fichier original.

write.infile(res.pca, "pca.csv", sep = ",")

Pour plus de détails sur les fondements de l’ACP, lire l’article d’où le code a été tiré : http://www.sthda.com/french/articles/38-methodes-des-composantes-principales-dans-r-guide-pratique/73-acp-analyse-en-composantes-principales-avec-r-l-essentiel/

2. Regroupement hiérarchique

La notion de regroupement hiérarchique recouvre différentes méthodes de clustering, et se catégorise en deux grandes familles : les méthodes ’ascendantes’ et les méthodes ’descendantes’.

L’approche ascendante, aussi appelée clustering agglomératif : on considère tout d’abord que chaque point est un cluster. Il y a donc autant de clusters que de points. Ensuite, on cherche les deux clusters les plus proches, et on les agglomère en un seul cluster. On répète cette étape jusqu’à ce que tous les points soient regroupés en un seul grand cluster. La fonction utilisée dans ce guide utilise l’approche ascendante.

L’approche descendante, aussi appelée clustering divisif : c’est l’inverse. On part d’un grand cluster contenant tous les points, puis on le divise successivement jusqu’à obtenir autant de clusters que de points.

On obtient donc une arborescence qui a un cluster tout en haut, et qui se divise petit à petit jusqu’à avoir autant de clusters que de points. On appelle cette arborescence un dendrogramme.Si vous avez beaucoup de points, le bas du dendrogramme risque d’être illisible, c’est pourquoi on ne représente parfois que le haut de l’arborescence.

Avant de pouvoir effectuer le test, il faut ouvrir le fichier contenant les données.

dataHC <- read.csv("C:/Users/WillzyX/Desktop/Maitrise/Atlas_FSL_2019_Stations_log_PCA_R.csv", header=T, sep=',',dec=".", stringsAsFactors = T)[,-1]

On peut ensuite préparer les données et les variables.

chimie <- as.matrix(dataHC[,c(2:3)])
rownames(chimie) = dataHC$Stations
head(chimie)
##                 PC1     PC2
## Batiscan    -2.0251  0.7026
## Becancour   -1.0975 -0.5721
## Chateauguay -0.0985  0.8750
## DuChene     -0.5585  0.0741
## DuLoup      -0.9694  0.3117
## Etchemin    -2.6240  1.5013
chimie <- na.omit(chimie)
chimie <- scale(chimie)

On peut ensuite procéder à la création du dendrogramme.

dend <- chimie %>%
  dist %>%
  hclust(method = "average") %>%
  as.dendrogram %>%
  set("branches_lwd",2) %>%
  set("labels_cex", 1.2) %>%
  set ("branches_k_col",
      value = c("#d73027","#fc8d59","#fee090","#ffffbf","#e0f3f8","#91bfdb","#4575b4"), k = 7) %>% 
  plot(main = "", ylim = c(-0.4,4))

Pour plus de détails sur les fondements du regroupement hiérarchique, lire l’article d’où le code a été tiré : http://www.sthda.com/english/wiki/beautiful-dendrogram-visualizations-in-r-5-must-known-methods-unsupervised-machine-learning

3. Random Forest

RF est l’un des algorithmes d’apprentissage supervisés les plus populaires et les plus performants, permettant de réaliser à la fois des régressions et des tâches de classification. Pour comprendre cet algorithme, il faut d’abord comprendre les fondements d’un arbre de décision. Un arbre de décision est un outil d’aide à la décision, représentant un ensemble de choix sous la forme graphique d’un arbre. Chaque nœud interne représente un test sur un attribut, chaque branche désigne les différents résultats possibles, et les « feuilles » de l’arbre, aussi nommés nœuds externes, représentent les étiquettes de classes qui peuvent être atteintes en fonction de décisions prises à chaque étape.

Avec RF, de nombreux arbres binaires sont construits à partir d’un sous-ensemble aléatoire des données en utilisant un rééchantillonnage bootstrap avec remplacement des individus. Ceux-ci sont appelés individus « in-bag ». Des individus « out-of-bag » sont ensuite utilisés pour mesurer la précision de prédiction de cet arbre, puisqu’ils ne sont pas utilisés pour construire un arbre spécifique. Pour un arbre donné, à chaque nœud, un sous-ensemble aléatoire de prédicteurs est recherché pour en trouver un qui maximise l’homogénéité au sein du groupe lorsqu’il est utilisé comme critère pour la division des nœuds. Chaque nouveau groupe est ensuite scindé à nouveau de manière récursive. Le processus récursif s’arrête lorsqu’une nouvelle division n’entraîne aucun gain supplémentaire d’homogénéité. Ainsi, deux niveaux de randomisation se produisent dans le processus RF : un dans la sélection initiale des individus (signatures de poissons) pour la construction de chaque arbre, et un dans la sélection des variables (éléments chimiques) utilisées pour la division des groupes à chaque nœud. RF construit une collection d’arbres de classification à partir de laquelle un consensus est établi : des individus « in-bag » sont descendus dans chaque arbre pour obtenir une seule prédiction d’origine, chaque arbre “votent” donc pour une classe. Les prédictions sur les arbres composant une forêt sont combinées selon des règles de majorité, RF choisira donc la classification ayant obtenu le plus de votes pour obtenir une prédiction finale de la forêt. Ainsi, on peut mesurer la précision de notre modèle, sous forme de pourcentage de bonne classification, en regardant la proportion d’individus « out-of-bag » qui ont été correctement classés.

Il faut d’abord aller chercher les données et les filtrer.

dataRF <- read.csv("/Users/Willzyx/Desktop/Maitrise/Atlas_FSL_2019_log_7Cluster_Moy_PCA_R.csv", header=T,
                 sep=',',dec=".", stringsAsFactors = T)[, -1]

trib<-filter(dataRF, Clusters == "A" | Clusters == "B" | Clusters ==  "C" |
               Clusters == "D" | Clusters == "E" | Clusters ==  "F" |
               Clusters == "G") %>% droplevels()

chim<-dataRF %>% drop_na(BaCa, MgCa, NaCa, SrCa)

Il faut ensuite déterminer les variables qui seront utilisées :

NUMLOOPS = 10
NTREES = 1000
PROPTRAINING = 0.75

confusion = NULL

Il faut créer la matrice d’affectation avant de pouvoir utiliser l’algorithme Random Forest.

stations = sort(as.character(unique(chim$Clusters)))
nstations = length(stations)
assignments = data.frame(matrix(0L, nrow=nstations, ncol=nstations))
names(assignments) = stations
rownames(assignments) = stations

On peut maintenant passer à travers les boucles de Random Forest.

for (i in 1:NUMLOOPS) {

  sample.ind = sample(2,nrow(chim),replace = T,prob = c(PROPTRAINING, 1 - PROPTRAINING))
  
  chim.train = chim[sample.ind==1,] 
  chim.test = chim[sample.ind==2,]
  
  chim_RF <- randomForest(x=chim.train[, c(2:5)],
                          y=chim.train[,1],
                          ntree=NTREES,
                          proximity = T,
                          oob.prox=TRUE)
  
  if (is.null(confusion)) {
    confusion = chim_RF$confusion
  } else {
    confusion = confusion + chim_RF$confusion
  }
  
  chim.test$predicted <- predict(chim_RF, chim.test[, c(2:5)])
  
  predicted_stations = as.data.frame(cbind(as.character(chim.test$Clusters),
                             as.character(chim.test$predicted)))
  
  print(head(predicted_stations))
  
# Le texte qui suit permet de compter combien de prédicitons ont été attribuées à chaque stations (ou clusters). La ligne précédente permet de visualiser comment l'algorithme procède (Valeur réelle en V1 et valeur prévue en V2).
  
  for (i in 1:nrow(predicted_stations)) {
      s = predicted_stations[i, 1]
      p = predicted_stations[i, 2]

      assignments[s, p] = assignments[s, p] + 1
  }
}
##   V1 V2
## 1  F  F
## 2  F  F
## 3  F  B
## 4  F  F
## 5  F  F
## 6  F  F
##   V1 V2
## 1  F  F
## 2  F  F
## 3  F  F
## 4  F  A
## 5  F  F
## 6  F  F
##   V1 V2
## 1  F  F
## 2  F  B
## 3  F  A
## 4  F  F
## 5  F  F
## 6  F  F
##   V1 V2
## 1  F  B
## 2  F  F
## 3  F  F
## 4  F  F
## 5  F  F
## 6  F  F
##   V1 V2
## 1  F  F
## 2  F  F
## 3  F  F
## 4  F  F
## 5  F  F
## 6  F  F
##   V1 V2
## 1  F  F
## 2  F  F
## 3  F  F
## 4  F  F
## 5  F  F
## 6  F  F
##   V1 V2
## 1  F  F
## 2  F  F
## 3  F  F
## 4  F  F
## 5  F  F
## 6  F  F
##   V1 V2
## 1  F  F
## 2  F  F
## 3  F  B
## 4  F  F
## 5  F  F
## 6  F  F
##   V1 V2
## 1  F  F
## 2  F  F
## 3  F  F
## 4  F  F
## 5  F  F
## 6  F  F
##   V1 V2
## 1  F  F
## 2  F  F
## 3  F  F
## 4  F  F
## 5  F  F
## 6  F  F

Les prochaines étapes permettent de prendre des décisions lors des tests. Ne pas effectuer cette partie lorsque vous ferez des boucles avec plus de 1 test, puisque les variables vont se référer à la valeur qu’ils avaient lors de la dernière boucle et non à la moyenne de toutes les boucles.

plot(chim_RF)

print(chim_RF)
## 
## Call:
##  randomForest(x = chim.train[, c(2:5)], y = chim.train[, 1], ntree = NTREES,      proximity = T, oob.prox = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 1000
## No. of variables tried at each split: 2
## 
##         OOB estimate of  error rate: 33.24%
## Confusion matrix:
##   A   B   C   D  E  F  G class.error
## A 7   0   0   0  0  7  0   0.5000000
## B 1 102  16  11  0  8 17   0.3419355
## C 0  12 118  23  1  0  3   0.2484076
## D 0  10  23 119  2  1 12   0.2874251
## E 0   3   1   8 19  2  1   0.4411765
## F 2   7   0   0  1 41  4   0.2545455
## G 0  29   7  10  3  2 50   0.5049505
varImpPlot(chim_RF, main =  "Importance relative des variables") 

chim_RF$importance
##      MeanDecreaseGini
## NaCa         128.6153
## MgCa         162.2333
## SrCa         131.8600
## BaCa         126.1519
MDSplot(chim_RF, chim.train$Clusters)

print(confusion / NUMLOOPS) # Performance de la reclassification. La classe.erro est juste, le problème est que c'est le nombre d'échantillons correspndant aux données d'entrainement qui sont dans chaque assignation.
##      A    B     C     D    E    F    G class.error
## A 10.1  0.0   0.0   0.0  0.0  4.6  0.0   0.3160827
## B  0.7 99.3  14.0   9.5  1.0  6.7 18.5   0.3374927
## C  0.0 13.2 123.1  25.0  0.7  0.0  4.6   0.2613744
## D  0.0 10.3  21.5 122.8  3.3  0.9  9.0   0.2682312
## E  0.0  2.3   1.0   6.9 19.0  2.0  0.4   0.3987355
## F  1.2  4.0   0.0   0.2  1.2 46.4  2.4   0.1628658
## G  0.0 24.5   6.4   9.9  1.5  1.5 59.1   0.4263419
print(assignments / NUMLOOPS) # Assignation donnees test. C'est la qu'on peut voir la performance du modele. Comme c'est un procede aleatoire il serait bien de générer une boucle (ex. 1000 fois) et utiliser la moyenne pour avoir une valeur plus fiable.
##     A    B    C    D   E    F    G
## A 3.6  0.1  0.0  0.0 0.0  1.6  0.0
## B 0.3 30.1  5.2  2.4 0.3  2.1  5.9
## C 0.0  5.4 43.5 10.2 0.2  0.0  1.1
## D 0.0  3.3  6.8 41.8 1.7  0.4  3.2
## E 0.0  0.2  0.0  2.3 7.1  0.4  0.4
## F 0.8  1.6  0.0  0.0 0.4 16.8  1.0
## G 0.0  8.4  1.4  4.2 0.1  0.7 19.3
print(apply(confusion / NUMLOOPS, 1, sum)) # Nombre d'échantillons moyen par station du groupe d'entrainement
##         A         B         C         D         E         F         G 
##  15.01608 150.03749 166.86137 168.06823  31.99874  55.56287 103.32634
print(apply(assignments / NUMLOOPS, 1, sum)) # Nombre d'échantillons moyen par station du groupe 'test'
##    A    B    C    D    E    F    G 
##  5.3 46.3 60.4 57.2 10.4 20.6 34.1
print(table(chim$Clusters) * (1 - PROPTRAINING)) # Devrait tendre vers ces chiffres (nombre attendu avec PROPTRAINING): 
## 
##     A     B     C     D     E     F     G 
##  5.00 49.00 56.75 56.25 10.50 19.00 34.25
print(sum(apply(assignments / NUMLOOPS, 1, sum))) # Nombre d'échantillons totaux moyens du groupe 'test'
## [1] 234.3