Atelier 9: Analyses multivariées

# Atelier 9: Analyses multivariées
## Série d’ateliers R
### Centre de la Science de la Biodiversité du Québec

---

# À propos de cet atelier
[![badge](https://img.shields.io/static/v1?style=for-the-badge&label=repo&message=dev&color=6f42c1&logo=github)](https://github.com/QCBSRworkshops/workshop09)
[![badge](https://img.shields.io/static/v1?style=for-the-badge&label=wiki&message=09&logo=wikipedia)](https://wiki.qcbs.ca/r_atelier9)
[![badge](https://img.shields.io/static/v1?style=for-the-badge&label=Diapos&message=09&color=red&logo=html5)](https://qcbsrworkshops.github.io/workshop09/workshop09-fr/workshop09-fr.html)
[![badge](https://img.shields.io/static/v1?style=for-the-badge&label=Diapos&message=09&color=red&logo=adobe-acrobat-reader)](https://qcbsrworkshops.github.io/workshop09/workshop09-fr/workshop09-fr.pdf)
[![badge](https://img.shields.io/static/v1?style=for-the-badge&label=script&message=09&color=2a50b8&logo=r)](https://qcbsrworkshops.github.io/workshop09/workshop09-fr/workshop09-fr.R)

---

# Packages requis

* [ape](https://cran.r-project.org/package=ape)
* [gclus](https://cran.r-project.org/package=gclus)
* [vegan](https://cran.r-project.org/package=vegan)

```R
install.packages(c('ape', 'gclus', 'vegan'))
```

---
# Objectifs d'apprentissage

##### Utiliser R pour faire des ordinations sans contraintes

#### Faire un dendrogramme avec R

---

# 1. Introduction
## Qu'est-ce que l'ordination?

---
# Une Dimension

Que se passe-t-il si nous voulons nous intéresser à la réponse de différentes espèces d'algues?

---
# Deux Dimensions

---
# Trois Dimensions

---
# 4,5,6, ou plus de Dimensions

---
# Ordination en espace réduit

---
# Ordination en espace réduit

- L'algèbre matricielle est complexe et difficile à comprendre

- Une compréhension générale est suffisante pour utiliser efficacement les méthodes d'ordination

---
# Méthodes pour la recherche scientifique

--
- **Questions / Hypothèses**
--

- **Design expérimental**

--
- **Collecte de données**
--

- **Transformation / Distance**
--

- **Analyses**
--

- **Rédaction**
--

- **Communication**

---
class: inverse, center, middle
# 2. Exploration des données

---
# Données de poissons de la rivière Doubs

Données de Verneaux (1973) :
- caractérisation des communautés de poissons
- 27 espèces
- 30 sites
- 11 variables environnementales

]

.pull.right[
![:scale 50%](images/DoubsRiver.png)
]

---
# Données de poissons de la rivière Doubs

Chargement des données espèces (`Doubs.Spe.csv`)

```r
spe <- read.csv("data/doubsspe.csv", row.names = 1)
spe <- spe[-8,] # supprimer le site vide
```

Chargement des données environnementales (`Doubs.Env.csv`)

```r
env <- read.csv("data/doubsenv.csv", row.names = 1)
env <- env[-8,] # remove site with no data
```

---
# Exporation des données

Explorer le contenu des données espèces :

```r
names(spe) # noms des objets
dim(spe) # dimensions
str(spe) # structure des objets
summary(spe) # résumé statistique
head(spe) # 6 premières lignes
```

```
#   CHA TRU VAI LOC OMB BLA HOT TOX VAN CHE BAR SPI GOU BRO PER BOU PSO ROT CAR
# 1   0   3   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
# 2   0   5   4   3   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
# 3   0   5   5   5   0   0   0   0   0   0   0   0   0   1   0   0   0   0   0
# 4   0   4   5   5   0   0   0   0   0   1   0   0   1   2   2   0   0   0   0
# 5   0   2   3   2   0   0   0   0   5   2   0   0   2   4   4   0   0   2   0
# 6   0   3   4   5   0   0   0   0   1   2   0   0   1   1   1   0   0   0   0
#   TAN BCO PCH GRE GAR BBO ABL ANG
# 1   0   0   0   0   0   0   0   0
# 2   0   0   0   0   0   0   0   0
# 3   0   0   0   0   0   0   0   0
# 4   1   0   0   0   0   0   0   0
# 5   3   0   0   0   5   0   0   0
# 6   2   0   0   0   1   0   0   0
```

---
# Fréquences des espèces

Observer la distribution de fréquence des espèces :

```r
ab <- table(unlist(spe))
barplot(ab, las = 1, col = grey(5:0/5),
 xlab = "Abondance des classes", ylab = "Fréquence")
```

---
# Fréquences des espèces

Combien de zéros?

```r
sum(spe == 0)
# [1] 408
```

Quelle proportion de zéros?

```r
sum(spe == 0)/(nrow(spe)*ncol(spe))
# [1] 0.5210728
```

---
# Richesse totale en espèce

Observer le nombre d'espèces présentes dans chaque site :

```r
site.pre <- rowSums(spe > 0)
barplot(site.pre, main = "Richesse spécifique",
 xlab = "Sites", ylab = "Nombre d'espèces",
 col = "grey ", las = 1)
```

---
# Comprenez vos données!

- Y-a-t-il beaucoup de zéros?

- Que veulent-ils dire?

.alert[Une mesure de 0 (e.g 0mg/L, 0°C) n'est pas équivalent à un 0 représentant une absence d'observation.]

---
# Avant de transformer vos données de composition des communautés...

--
- abondances/comptes/présence-absence relatives?

--
- distributions asymmétriques ?

--
- beaucoup d'espèces rares?

--
- surabondance d'espèces dominantes?

--
- problème de double Zéro?

---
# Transformer les données de composition des communautés

---
# Transformer les données de composition des communautés

## Exemples

Transformer des comptes en présence - absence

```r
library(vegan)
spec.pa <- decostand(spe, method = "pa")
```

Réduire le poids des espèces rares

```r
spec.hel <- decostand(spe, method = "hellinger")
spec.chi <- decostand(spe, method = "chi.square")
```

Réduire le poids des espèces abondantes

```r
spe.pa <- decostand(spe, method = "log")
```

---
# Données sur l'environnement

```r
names(env) # Names of objects
dim(env) # dimensions
str(env) # structure of objects
summary(env) # summary statistics
head(env) # first 6 rows
```

```r
head(env) # first 6 rows
#    das alt  pen  deb  pH dur  pho  nit  amm  oxy dbo
# 1  0.3 934 48.0 0.84 7.9  45 0.01 0.20 0.00 12.2 2.7
# 2  2.2 932  3.0 1.00 8.0  40 0.02 0.20 0.10 10.3 1.9
# 3 10.2 914  3.7 1.80 8.3  52 0.05 0.22 0.05 10.5 3.5
# 4 18.5 854  3.2 2.53 8.0  72 0.10 0.21 0.00 11.0 1.3
# 5 21.5 849  2.3 2.64 8.1  84 0.38 0.52 0.20  8.0 6.2
# 6 32.4 846  3.2 2.86 7.9  60 0.20 0.15 0.00 10.2 5.3
```
Explorer la colinéarité en visualisant les corrélations entre les variables

```r
pairs(env, main = "Bivariate Plots of the Environmental Data")
```

---
# Données sur l'environnement

![](images/EnvDat1.png)

---
# Standardisation

Standardiser les variables environnementales est indispensable car il est impossible de comparer des variables d'unités différentes :

```r
## ?decostand
env.z <- decostand(env, method = "standardize")
```
Cette fonction centre-réduit les données pour permettre la fiabilité des analyses :

```r
apply(env.z, 2, mean)
#           das           alt           pen           deb            pH 
# -7.959539e-17 -4.795165e-17  2.494600e-17 -7.323225e-17 -1.730430e-15 
#           dur           pho           nit           amm           oxy 
# -2.028505e-16  4.445790e-17  2.875893e-17  2.754434e-17 -4.038167e-16 
#           dbo 
#  9.829975e-17
apply(env.z, 2, sd)
# das alt pen deb  pH dur pho nit amm oxy dbo 
#   1   1   1   1   1   1   1   1   1   1   1
```

---
class: inverse, center, middle
# 3. Similarité / Dissimilarité

---
# Mesure d'association

L'algébre matricielle est au coeur de plusieurs méthodes d'analyses multivariées

- Explorer différentes mesures de distance entre objets permet de mieux comprendre le fonctionnement de l'ordination

---
# Au-délà de la 1ère dimension

- Les jeux de données écologiques correspondent souvent à de grandes matrices

- L'ordination calcule les relations entre espèces, ou entre objets

- Ces relations peuvent être simplifiées par des mesures de dissimilarités
]

![:scale 40%](images/PCAMatrix.png)

![:scale 40%](images/distMes.png)

![:scale 40%](images/distMat.png)
]

---
# Similarité / Dissimilarité

- Utile pour comprendre vos données
- Certains types d'ordination ou de groupement nécessitent des mesures appropriées

![](images/similarity.png)

---
# Mesures de distance des communautés

--

---
# Comparaison des sites de la rivière Doubs

La fonction `vegdist()` comprend les mesures de distances communes :

```r
?vegdist
```
Comment la composition des communautés diffère-t-elle entre les 30 sites de la rivière Doubs?

```r
spe.db.pa <- vegdist(spe, method = "bray")
```

---
# Comparaison des sites

---
# Comparaison des sites

---
# Visualisation d'une matrice de distances

---
# Défi #1 ![:cube]()

Discuter avec votre voisin:

.center[**Comment savoir si deux objets caractérisés par des données multidimensionnelles sont similaires?**]

- Faites une liste de vos suggestions

---
# Et qu'en est-il de l'ordination?

Avec des méthodes d'ordination, nous ordonnons vos objets (sites) en fonction de leur similarité

- Plus les sites sont similaires, plus ils sont proches dans l'espace d'ordination (plus petites distances)

- En écologie, on calcule habituellement la similarité entre sites en fonction de leur composition en espèces ou de leur conditions environnementales.

---
# Analyse schématique des analyses multivariées

---
# Groupement

- Permet de mettre en lumière des structures dans les données en partitionnant les objets

- Les résultats sont représentés sous forme de dendrogramme (arbre)

- Pas une méthode statistique!

---
# Aperçu de 3 méthodes hiérarchiques

- Groupement agglomératif à liens simples

- Groupement agglomératif à liens complets

- Groupement de Ward

- Les éléments de petits ensembles se regroupent en groupes plus vastes de rang supérieur
   - (e.g. espèces, genres, familles, ordres...)

---
# Groupement hiérarchique

À partir d'une matrice de distances, on classe les objets en ordre croissant

![](images/Hierachic1.png)

---
# Groupement à liens simples

![:scale 50%](images/singleClust1.png)
--

]

- Les deux objets les plus proches se regroupent

- Ensuite les deux objets les plus proches suivants

- et ainsi de suite.

![](images/singleClust2.png)

]

---
# Groupement à liens complet

![:scale 50%](images/compleClust1.png)

]

- Les deux objets les plus proches se regroupent

- Ensuite les groupes se lient à la distance à laquelle les objets qu'ils contiennent sont tous liés

![](images/compleClust2.png)
]

---
# Comparaison

Créer une matrice de distance à partir des données de la rivière Doubs transformées Hellinger et faire le groupement à liens simples :

```r
spe.dhe1 <- vegdist(spec.hel, method = "euclidean")
spe.dhe1.single <- hclust(spe.dhe1, method = "single")
plot(spe.dhe1.single)
```

---
# Comparaison

![](images/comparison.png)

**Liens simples :**

Les objets ont tendance à s'enchaîner (e.g. 19,29,30,26)
]

**Liens complets :**
Les groupes sont plus distincts
]

---
# Groupement de Ward

- Utilise la méthode des moindres carrés pour lier les objets
  - les groupes fusionnent de façon à minimiser la variance intragroupe
  - à chaque étape, la paire de groupes à fusionner est celle qui résulte à la plus petite augmentation de la somme des carrés des écarts intra-groupes

---
# Groupement de Ward

Faire le groupement de Ward et dessiner le dendrogramme en utilisant la racine carrée des distances :

```r
spe.dhel.ward <- hclust(spe.dhe1, method = "ward.D2")
spe.dhel.ward$height <- sqrt(spe.dhel.ward$height)
plot(spe.dhel.ward, hang = -1) # hang = -1 aligns objects at the same level
```

---
# Groupement de Ward

<img src="workshop09-fr_files/figure-html/unnamed-chunk-21-1.png" width="648" style="display: block; margin: auto;" />
Les objets ont tendance à former des groupes plus sphériques et homogènes

---
# Comment choisir la bonne méthode ?

- Dépend de votre objectif
  - démontrer des gradients? des contrastes?
- Si plus d'une méthode semble adéquate, comparer les dendrogrammes
- Encore une fois : ceci **n'est pas** une méthode statistique
 Mais! il est possible de:
  - déterminer le nombre de groupe optimal
  - faire des tests statistiques sur les résultats
  - combiner le groupement à l'ordination pour distinguer des groupes de sites
---
class: inverse, center, middle
# 4. Ordination non contrainte

---
# Définitions

- **Variance:** mesure de la dispersion d'une variable `$y_j$` de sa moyenne
--

- **Co-variance:** mesure de co-dispersion des variables `$y_j$` et `$y_i$` de leur moyenne

--
- **Corrélation:** mesure de la force du lien entre 2 variables : `$r_{ij} = (d_{ij} / d_j x d_k)$`

--
- **Valeurs propres:** proportion de variance (dispersion) représentée par un axe d'ordination

--
- **Orthogonalité:** angle droit entre 2 axes ou 2 flèches, ce qui veut dire qu'ils sont indépendants = non corrélés

--
- **Score:** position d'un point sur un axe. Tous les scores d'un point donnent ses coordonnées dans l'espace multidimensionnel. Ils peuvent être utilisés pour d'autres analyses (e.g combinaison linéaire de variables mesurées)

--
- **Dispersion** (inertie): Mesure de la variabilité totale du diagramme de dispersion de l'espace multidimensionnel en fonction de son centre de gravité

---
# Ordination non contrainte

- Évalue la relation **dans** un ensemble de variables (espèces ou variables environnementales, et non **parmi** les ensembles, i.e analyse sous contraintes)

- Trouve les composants clés de la variation entre échantillons, sites, espèces, etc...

- Réduit le nombre de dimensions dans les données multivariées sans perte d'informations considérables

- Créer de nouvelles variables pour des analyses ultérieures (comme la régression)

---
# 4.1. Analyse en Composantes Principales (ACP ou PCA)

- Préserve, en 2 dimensions, le maximum de variation des données
- Il en résulte des variables synthétiques orthogonales entre elles (et donc non corrélées)

---
# PCA - Ce qu'il vous faut

- Un jeu de données correspondant à des variables réponses (eg. composition de communautés) OU à des variables explicatives (e.g variables environnementales)

**PAS LES DEUX!**

.pull-left[
- Échantillons correspondant à des mesures du même jeu de variables
- Généralement, un jeu de données plus long que large est préféré
]

---
# PCA - Principes

|Site|Species 1| Species 2|
|---|------|------|
|A|7|3|
|B|4|3|
|C|12|10|
|D|23|11|
|E|13|13|
|F|15|16|
|G|18|14|

---
# PCA - Principes

.small[
.alert[En 2D, les sites seraient disposés de cette façon... Notez la dispersion dans le diagramme de dispersion]]

---
# PCA - Principes

.small[
.alert[La première composantes principale est celle qui maximise la variation observée... la meilleure droite entre les sites]]

---
# PCA - Principes

.center[
![:scale 70%](images/PCA2.png)]
.small[
.alert[La seconde composante principale est ajoutée perpendiculairement au premier axe]]

---
# PCA - Principes

.small[Le graphique final subit une rotation afin que les deux axes correspondent aux composantes principales (et non plus aux espèces)]
---
# PCA - Cas multidimensionnel

- **PC1** --> axe qui maximise la variance des points projetés perpendiculairement sur les axes.
- **PC2** --> doit être orthogonal à PC1, mamis sa direction maximise la variance des points projetés.
- **PC3** --> et ainsi de suite : orthogonale à PC1 et PC2...

.alert[Quand il y a plus de deux dimensions, la PCA produit un nouvel espace dans lequel tous les axes sont orthogonaux (i.e. la corrélation entre les axes =0) et où les axes sont ordonnés selon le pourcentage de variation des données brutes qu'ils représentent (valeur propre).]

---
# PCA - Essayons sur les données Poissons!

- La PCA (tout comme la RDA) est implémentée par la fonction `rda()` de la librairie vegan

- Effectuer une PCA sur les abondances de poissons transformées Hellinger

```r
spe.h.pca <- rda(spec.hel)

summary(spe.h.pca)
# 
# Call:
# rda(X = spec.hel) 
# 
# Partitioning of variance:
#               Inertia Proportion
# Total          0.5025          1
# Unconstrained  0.5025          1
# 
# Eigenvalues, and their contribution to the variance 
# 
# Importance of components:
#                          PC1     PC2     PC3     PC4     PC5     PC6     PC7
# Eigenvalue            0.2580 0.06424 0.04632 0.03850 0.02197 0.01675 0.01472
# Proportion Explained  0.5133 0.12784 0.09218 0.07662 0.04371 0.03334 0.02930
# Cumulative Proportion 0.5133 0.64118 0.73337 0.80999 0.85370 0.88704 0.91634
#                           PC8      PC9     PC10     PC11     PC12     PC13
# Eigenvalue            0.01156 0.006936 0.006019 0.004412 0.002982 0.002713
# Proportion Explained  0.02300 0.013803 0.011978 0.008781 0.005935 0.005399
# Cumulative Proportion 0.93934 0.953144 0.965123 0.973903 0.979838 0.985237
#                           PC14     PC15     PC16      PC17      PC18      PC19
# Eigenvalue            0.001835 0.001455 0.001118 0.0008309 0.0005415 0.0004755
# Proportion Explained  0.003651 0.002895 0.002225 0.0016535 0.0010776 0.0009463
# Cumulative Proportion 0.988888 0.991783 0.994008 0.9956612 0.9967389 0.9976852
#                            PC20      PC21      PC22      PC23      PC24
# Eigenvalue            0.0003680 0.0002765 0.0002253 0.0001429 7.618e-05
# Proportion Explained  0.0007324 0.0005503 0.0004483 0.0002845 1.516e-04
# Cumulative Proportion 0.9984176 0.9989678 0.9994161 0.9997006 9.999e-01
#                           PC25      PC26      PC27
# Eigenvalue            4.99e-05 1.526e-05 9.118e-06
# Proportion Explained  9.93e-05 3.036e-05 1.814e-05
# Cumulative Proportion 1.00e+00 1.000e+00 1.000e+00
# 
# Scaling 2 for species and site scores
# * Species are scaled proportional to eigenvalues
# * Sites are unscaled: weighted dispersion equal on all dimensions
# * General scaling constant of scores:  1.93676 
# 
# 
# Species scores
# 
#          PC1      PC2       PC3        PC4        PC5       PC6
# CHA  0.17336  0.08295 -0.064963  0.2539861 -0.0285801  0.019057
# TRU  0.64860  0.01162 -0.261994 -0.1606020 -0.0745819 -0.088616
# VAI  0.51810  0.14773  0.165304  0.0241017  0.1012928  0.104748
# LOC  0.38606  0.16615  0.242995 -0.0275216  0.1258011  0.048299
# OMB  0.16893  0.06274 -0.096143  0.2426514  0.0140574  0.062117
# BLA  0.07786  0.14644 -0.031402  0.2339394 -0.1032338 -0.040810
# HOT -0.18491  0.04901 -0.045107  0.0199377  0.0687305  0.009650
# TOX -0.14644  0.17834 -0.010937  0.0649955 -0.0006229 -0.106955
# VAN -0.11436  0.15673  0.142223 -0.0127266 -0.1989404  0.013897
# CHE -0.09682 -0.15449  0.242943  0.1124210  0.0233830 -0.039996
# BAR -0.19826  0.21211 -0.053980  0.0969899  0.0067098 -0.035442
# SPI -0.17689  0.16250 -0.033112  0.0397113  0.0323159 -0.072908
# GOU -0.23138  0.09782  0.064144 -0.0013887 -0.1503303  0.130575
# BRO -0.15129  0.12804  0.040303 -0.1203826 -0.1006077  0.066242
# PER -0.15719  0.18144  0.057029 -0.0940032 -0.0412984 -0.060409
# BOU -0.22853  0.13870 -0.062197 -0.0125024  0.0798647 -0.006907
# PSO -0.22790  0.08231 -0.065797  0.0172143  0.0611434 -0.001407
# ROT -0.19221  0.03090 -0.006264 -0.0739133 -0.0731548  0.074581
# CAR -0.18699  0.13388 -0.050804  0.0001803  0.0403961 -0.031005
# TAN -0.19169  0.15719  0.114415 -0.0818330  0.0142624 -0.072024
# BCO -0.20174  0.08807 -0.067086 -0.0529106  0.0737228  0.037312
# PCH -0.14717  0.05829 -0.067311 -0.0458414  0.0501013  0.031605
# GRE -0.30155 -0.01785 -0.084333 -0.0181797  0.0226500  0.126639
# GAR -0.35245 -0.14076  0.168014  0.0185946  0.0213462 -0.129788
# BBO -0.24317  0.03679 -0.082731 -0.0384489  0.0939828  0.063369
# ABL -0.42536 -0.26155 -0.054190  0.1021959 -0.0078085  0.044540
# ANG -0.20631  0.11889 -0.062079 -0.0175733  0.0718743 -0.001956
# 
# 
# Site scores (weighted sums of species scores)
# 
#          PC1      PC2      PC3      PC4       PC5       PC6
# 1   0.367401 -0.39935 -1.08857 -0.63304 -0.512027 -0.858378
# 2   0.503582 -0.05683 -0.19259 -0.43441  0.389533  0.069451
# 3   0.461709  0.02262 -0.06522 -0.49798  0.309425  0.270577
# 4   0.298336  0.15130  0.26748 -0.53196  0.003088  0.184821
# 5  -0.002222  0.07631  0.54769 -0.50936 -0.780261 -0.169353
# 6   0.212816  0.08345  0.55091 -0.42210 -0.139518 -0.104278
# 7   0.438055 -0.06114  0.15590 -0.31150  0.158686  0.036565
# 9   0.040794 -0.44269  0.89022  0.09609  0.641193 -0.646943
# 10  0.298011 -0.01094  0.56837 -0.10013 -0.088124  0.515072
# 11  0.467609 -0.12622 -0.15505  0.29459  0.325464  0.200912
# 12  0.476845 -0.07691 -0.16329  0.29384  0.360112  0.194576
# 13  0.483620  0.06649 -0.44723  0.53734  0.048587  0.182565
# 14  0.371728  0.16555 -0.21939  0.62130 -0.183604  0.364847
# 15  0.277048  0.23525  0.08928  0.61773 -0.475769  0.124107
# 16  0.077024  0.47455  0.17116  0.34361 -0.570434 -0.572740
# 17 -0.053860  0.42290  0.02810  0.42376 -0.059203 -0.586419
# 18 -0.135418  0.37780  0.03233  0.39706 -0.007199 -0.347064
# 19 -0.269281  0.30751  0.18022  0.09354  0.178657 -0.016299
# 20 -0.378830  0.19764  0.04939 -0.03438  0.157660 -0.056696
# 21 -0.409369  0.22888 -0.08401 -0.12823  0.152787  0.096105
# 22 -0.443679  0.17698 -0.13708 -0.13152  0.103294  0.030004
# 23 -0.242292 -1.11711  0.15254  0.40512  0.045573 -0.576778
# 24 -0.358333 -0.83372 -0.17314  0.27200  0.181192  0.347231
# 25 -0.325288 -0.61983  0.10487  0.01059 -1.034438  0.750325
# 26 -0.441703  0.02111 -0.13742 -0.14346  0.200775  0.244356
# 27 -0.444529  0.12735 -0.15915 -0.14112  0.179240  0.123487
# 28 -0.446407  0.12774 -0.18830 -0.15467  0.239617  0.117101
# 29 -0.355788  0.28044 -0.28006 -0.02003  0.110181  0.079568
# 30 -0.467578  0.20086 -0.29797 -0.21269  0.065512  0.003276
```

---
# Fonction `rda()`

- RDA en 2 étapes :

- régressions multiples
  - PCA sur les valeurs régressées

- Si on donne seulement un tableau à la fonction `rda()`, la fonction roule une PCA sans faire les régressions

rda(Y) ou rda(X) ![:faic](arrow-right) PCA ]]