Atelier 9: Analyses multivariées

.title[
# Atelier 9: Analyses multivariées
]
.subtitle[
## Série d’ateliers R
]
.author[
### 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=Présentation&message=09&color=BF616A)](https://r.qcbs.ca/workshop09/pres-fr/workshop09-pres-fr.html)    [![badge](https://img.shields.io/static/v1?style=for-the-badge&label=Script&message=09&color=D08770&logo=r)](https://r.qcbs.ca/workshop09/book-fr/workshop09-script-fr.R) 
[![badge](https://img.shields.io/static/v1?style=for-the-badge&label=Livre&message=09&color=EBCB8B)](https://r.qcbs.ca/workshop09/book-fr/index.html)
[![badge](https://img.shields.io/static/v1?style=for-the-badge&label=Site&message=r.qcbs.ca&color=A3BE8C)](https://r.qcbs.ca/fr/workshops/r-workshop-09/)
[![badge](https://img.shields.io/static/v1?style=for-the-badge&label=GitHub&message=09&color=B48EAD&logo=github)](https://github.com/QCBSRworkshops/workshop09)

---

en modifiant et en améliorant son contenu dans le cadre du Prix d'apprentissage et de développement (*Le*arning *a*nd *D*evelopment *A*ward)
]

**2022** - **2021** - **2020**

Pedro Henrique P. Braga

Katherine Hébert

Mi Lin

Linley Sherin

**2019** - **2018** - **2017**

Gabriel Muñoz

Marie Hélène-Brice

Pedro Henrique P. Braga

]
]

**2016** - **2015** - **2014**

Bérenger Bourgeois

Xavier Giroux-Bougard

Amanda Winegardner

Emmanuelle Chrétien

Monica Granados

]
]

.center[
__Si vous voulez contribuer aussi__, visitez [r.qcbs.ca/fr/contributing](https://r.qcbs.ca/fr/contributing/) et n'hésitez pas à [nous contacter](mailto:csbq.qcbs.r@gmail.com)!
]

---

# Matériel requis

Cet atelier requiert les dernières versions de [RStudio](https://rstudio.com/products/rstudio/download/#download) et de [R](https://cran.rstudio.com/).

.pull-left[
Vous devez également avoir installé ces paquets:
* [ape](https://cran.r-project.org/package=ape)
* [ade4](https://cran.r-project.org/package=ade4)
* [codep](https://cran.r-project.org/package=codep)
* [gclus](https://cran.r-project.org/package=gclus)
* [vegan](https://cran.r-project.org/package=vegan)
* [GGally](https://cran.r-project.org/package=GGally)
* [PlaneGeometry](https://cran.r-project.org/package=PlaneGeometry)
* [remotes](https://cran.r-project.org/package=remotes)

]

Pour les installer à partir de CRAN, roulez:

```r
install.packages(c("ape",
                   "ade4",
                   "codep",
                   "gclus",
                   "vegan",
                   "GGally",
                   "PlaneGeometry",
                   "remotes"))
```
]

.pull-left2[
Au cours de cet atelier, il y aura une série de **défis** que vous pouvez reconnaître par ce cube de Rubik.

**Pendant ces défis, n'hésitez pas à collaborer!**
]

]
]

---

# Objectifs d'apprentissage

1. Apprendre les bases de l'analyse multivariée pour révéler des modèles dans les données sur la composition des communautés.

2. Utiliser `R` pour effectuer une ordination sans contrainte.

3. Apprendre les coefficients de similarité et de dissimilarité et les transformations des données.

4. Utiliser `R` pour créer des dendrogrammes.

5. Apprendre les méthodes suivantes :

* Analyse de regroupement
  * Analyse en composantes principales (ACP)
  * Analyse en coordonnées principales (ACoP)
  * Échelle multidimensionnelle non métrique (PMnM).

---

# 1. Préambule

---

# Récapitulation : Analyses univariées

Nous avons appris une multitude d'analyses qui nous ont permis d'interpréter des données écologiques en décrivant les effets d'_une ou plusieurs_ variables sur _une_ variable de réponse.

1. Modèles linéaires généraux
  1. `lm()` ;
  2. `anova()` ;
  3. `t.test()` ;
  4. `lmer()`.

2. Modèles linéaires généralisés
  1. `glm()` et `glmer()` avec plusieurs fonctions de liaison `family()`.

3. Modèles additifs généralisés
  1. `gam()`.
]

1. Quels sont les effets des précipitations et de la température sur la richesse des espèces ?
2. Comment l'abondance des microbes change-t-elle entre les hôtes ?
3. Les poissons cooccurrents deviennent-ils plus agressifs après avoir été induits à la peur ?
]

---

# Statistiques multivariées

Cependant, on peut être intéressé à faire des inférences à partir de données écologiques contenant _plus d'une_ variable dépendante.

Cet intérêt peut être motivé par la validation d'hypothèses et la modélisation, mais peut aussi être entièrement exploratoire.

Par exemple, notre question de recherche pourrait être la suivante

1. Comment la _composition bactérienne_ des feuilles d'érable change-t-elle le long du gradient d'altitude ?

2. Quelle est la _dissimilarité compositionnelle_ des communautés de chauves-souris ?

3. Quel est le rapport entre les communautés locales d'araignées et leur _composition_ ?

Dans tous ces cas, le résultat est composé de plusieurs variables, _e.g._ généralement une matrice échantillon par espèce ou échantillon par environnement.

---
# Statistiques multivariées

Nous allons maintenant plonger dans les **statistiques multivariées**, un ensemble d'outils qui nous permettront de répondre à des questions nécessitant l'observation ou l'analyse simultanée de plus d'une variable de résultat.

---

# Contexte des méthodes multivariées

1. Mesures et matrices d'association (ou de dissimilarité)

2. Analyse de classification (ou _cluster_)

3. Ordination sans contrainte

4. Ordination contrainte (ou canonique)

---

# Algèbre matricielle : un résumé très bref

L'algèbre matricielle est bien adaptée à l'écologie, car la plupart des _ensembles de données_ avec lesquels nous travaillons sont dans un format _matriciel_.

.pull-left[
Les tableaux de données écologiques sont obtenus sous forme d'observations d'objets ou d'unités d'échantillonnage, et sont souvent enregistrés comme ceci :

| Objets | `$y_1$` | `$y_2$` | `$\dots$` | `$y_n$` |
| :-------------: |:-------------:| :-----:|:-----:|:-----:|
| `$x_1$` | `$y_{1,1}$` | `$y_{1,2}$` | `$\dots$` | `$y_{1,n}$` |
| `$x_2$` | `$y_{2,1}$` | `$y_{2,2}$` | `$\dots$` | `$y_{2,n}$` |
| `$\vdots$` | `$\vdots$` | `$\vdots$` | `$\ddots$` | `$\vdots$` |
| `$x_m$` | `$y_{m,1}$` | `$y_{m,2}$` | `$\dots$` | `$y_{m,n}$` |
 
où `$x_m$` est l'unité d'échantillonnage `$m$` ; et `$y_n$` est le descripteur écologique qui peut être, _e.g._, l'espèce présente ou une variable chimique.
]

Le même tableau de données écologiques peut être représenté en _notation matricielle_ de la manière suivante :
`$$Y = [y_{m,n}] =
\begin{bmatrix}
y_{1,1} & y_{1,2} & \cdots & y_{1,n} \\
y_{2,1} & y_{2,2} & \cdots & y_{2,n} \\
\vdots  & \vdots  & \ddots & \vdots  \\
y_{m,1} & y_{m,2} & \cdots & y_{m,n}
\end{bmatrix}$$`

où les lettres minuscules indiquent les _éléments_, et les lettres en indice indiquent la _position de ces éléments_ dans la matrice (et dans le tableau !).
]

---

# Algèbre matricielle : un résumé très bref

L'algèbre matricielle est bien adaptée à l'écologie, car la plupart des _ensembles de données_ avec lesquels nous travaillons sont dans un format _matriciel_.

.center[_une matrice de lignes_]
`$$\begin{bmatrix}
y_{1,1} & y_{1,2} & \cdots & y_{1,n} \\
\end{bmatrix}$$`

.center[_une matrice à colonnes_]
`$$\begin{bmatrix}
y_{1,1} \\ y_{2,2} \\ \vdots \\ y_{m,2}
\end{bmatrix}$$`

]

.pull-right[
Le même tableau de données écologiques peut être représenté en _notation matricielle_ de la manière suivante :
`$$Y = [y_{m,n}] =
\begin{bmatrix}
y_{1,1} & y_{1,2} & \cdots & y_{1,n} \\
y_{2,1} & y_{2,2} & \cdots & y_{2,n} \\
\vdots  & \vdots  & \ddots & \vdots  \\
y_{m,1} & y_{m,2} & \cdots & y_{m,n}
\end{bmatrix}$$`

où les lettres minuscules indiquent les _éléments_, et les lettres en indice indiquent la _position de ces éléments_ dans la matrice (et dans le tableau !).
]

---
# Algèbre matricielle : un résumé très bref

###### Matrices d'association

Deux matrices importantes peuvent être dérivées de la matrice des données écologiques : la _**matrice d'association entre objets**_ et la _**matrice d'association entre descripteurs**_.

<div class="math">
\[ Y =
\begin{array}{cc}
\begin{array}{ccc}
x_1 \rightarrow\\
x_2 \rightarrow\\
\vdots \\
x_m \rightarrow\\
\end{array}
&
\begin{bmatrix}
y_{1,1} & y_{1,2} & \cdots & y_{1,n} \\
y_{2,1} & y_{2,2} & \cdots & y_{2,n} \\
\vdots & \vdots & \ddots & \vdots \\
y_{m,1} & y_{m,2} & \cdots & y_{m,n}
\end{bmatrix}
\end{array}
\]
</div>

on peut examiner la relation entre les deux premiers objets :

<div class="math">
\[x_1 \rightarrow \begin{bmatrix}
y_{1,1} & y_{1,2} & \cdots & y_{1,n} \\
\end{bmatrix}
\]
</div>

<div class="math">
\[x_2 \rightarrow
\begin{bmatrix}
y_{2,1} & y_{2,2} & \cdots & y_{2,n} \\
\end{bmatrix}
\]
</div>

et obtenir $a_{1,2}$.

]

Nous pouvons remplir la matrice d'association `$A_{n,n}$` avec les relations entre tous les objets de `$Y$` :

<div class="math">
\[A_{n,n} =
\begin{bmatrix}
a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\
a_{2,1} & a_{2,2} & \cdots & a_{2,n} \\
\vdots & \vdots & \ddots & \vdots \\
a_{n,1} & a_{n,2} & \cdots & a_{n,n}
\end{bmatrix}\]
</div>

Parce que $A_{n,n}$ a le même nombre de lignes et de colonnes, elle est notée matrice carrée.

Par conséquent, elle possède $n^2$ éléments.

]

---
# Algèbre matricielle : un résumé très bref

###### Matrices d'association

Deux matrices importantes peuvent être dérivées de la matrice des données écologiques : la _**matrice d'association entre objets**_ et la _**matrice d'association entre descripteurs**_.

Ces matrices sont la base des analyses **_mode Q_** et **_mode R_** en écologie.

.pull-left[
`$$Y_{m,n} =
\begin{bmatrix}
y_{1,1} & y_{1,2} & \cdots & y_{1,n} \\
y_{2,1} & y_{2,2} & \cdots & y_{2,n} \\
\vdots  & \vdots  & \ddots & \vdots  \\
y_{m,1} & y_{m,2} & \cdots & y_{m,n}
\end{bmatrix}$$`
]

`$$A_{m,m} =
\begin{bmatrix}
a_{1,1} & a_{1,2} & \cdots & a_{1,m} \\
a_{2,1} & a_{2,2} & \cdots & a_{2,m} \\
\vdots  & \vdots  & \ddots & \vdots  \\
a_{m,1} & a_{m,2} & \cdots & a_{m,m}
\end{bmatrix}$$`

]

`$$A_{n,n} =
\begin{bmatrix}
a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\
a_{2,1} & a_{2,2} & \cdots & a_{2,n} \\
\vdots  & \vdots  & \ddots & \vdots  \\
a_{n,1} & a_{m,2} & \cdots & a_{n,n}
\end{bmatrix}$$`

]

`$\leftarrow$` Analyse **_mode R_** pour les descripteurs ou les espèces

]

---
## Nos projets pour les prochaines étapes...

Nous allons nous plonger dans les analyses en **R-mode** et **Q-mode**, et nous allons explorer :

1. Les coefficients d'association : dissimilarité et similarité
2. La transformation des données de composition
3. Analyses de regroupement
4. Analyses d'ordination dans l'espace réduit

---
# Communautés de poissons de la rivière Doubs

.pull-left3[
Verneaux (1973) a proposé d'utiliser les espèces de poissons pour caractériser les zones écologiques le long des cours d'eau européens.

Il a recueilli des données dans **30 sites** le long de la rivière Doubs près de la frontière franco-suisse, dans les montagnes du Jura. Il a montré que les communautés de poissons étaient des indicateurs biologiques de ces masses d'eau.
]

.pull.right3[
.center[![:scale 25%](images/DoubsRiver.png)]
]

1. L'abondance de 27 espèces de poissons dans les communautés ;
2. Les variables environnementales enregistrées sur chaque site ; et,
3. Les coordonnées géographiques de chaque site.
]

.pull.right3[
.center[![:scale 30%](https://upload.wikimedia.org/wikipedia/commons/c/cc/Doubs_Laissey.jpg)]

.xsmall[Verneaux, J. (1973) _Cours d'eau de Franche-Comté (Massif du Jura). Recherches écologiques sur le réseau hydrographique du Doubs_. Essai de biotypologie. Thèse d'état, Besançon. 1–257.]
]

---
# Communautés de poissons de la rivière Doubs

Vous pouvez télécharger les jeux de données sur [r.qcbs.ca/fr/workshops/r-workshop-09](http://r.qcbs.ca/fr/workshops/r-workshop-09/).

```r
spe <- read.csv("data/doubsspe.csv",
 row.names = 1)
env <- read.csv("data/doubsenv.csv",
 row.names = 1)
```

]

Leurs données peuvent également être récupérées dans le paquet `ade4` :

```r
library (ade4)
data (doubs)
spe <- doubs$fish
env <- doubs$env
```

Alternativement, à partir du paquet `codep` :

```r
library (codep)
data (Doubs)
spe <- Doubs.fish
env <- Doubs.env
```
]

---
# Communautés de poissons de la rivière Doubs

Jetons un coup d'œil aux données `spe` :

```r
head(spe)[, 1:8]
#   CHA TRU VAI LOC OMB BLA HOT TOX
# 1   0   3   0   0   0   0   0   0
# 2   0   5   4   3   0   0   0   0
# 3   0   5   5   5   0   0   0   0
# 4   0   4   5   5   0   0   0   0
# 5   0   2   3   2   0   0   0   0
# 6   0   3   4   5   0   0   0   0
```

```r
str(spe)
...
# 'data.frame':	30 obs. of  27 variables:
#  $ CHA: int  0 0 0 0 0 0 0 0 0 0 ...
#  $ TRU: int  3 5 5 4 2 3 5 0 0 1 ...
#  $ VAI: int  0 4 5 5 3 4 4 0 1 4 ...
#  $ LOC: int  0 3 5 5 2 5 5 0 3 4 ...
#  $ OMB: int  0 0 0 0 0 0 0 0 0 0 ...
#  $ BLA: int  0 0 0 0 0 0 0 0 0 0 ...
#  $ HOT: int  0 0 0 0 0 0 0 0 0 0 ...
...
```
]

```r
# Essayez-en quelques-uns !
names(spe)   # noms d'objets
dim(spe)     # dimensions
str(spe)     # structure des objets
summary(spe) # sommaire
head(spe)    # debut
```

```r
# Nous devrons enlever des sites avec 
# aucune spèce

# Obtenir les lignes avec qui ont des 
# valeurs differents de zero
row_sub <- apply(spe, 
 1, function(row) any(row !=0 ))

# Garder ces lignes dans `spe`
spe <- spe[row_sub, ]
```

]

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

Nous pouvons alors explorer les objets contenant nos données nouvellement chargées. Jetons un coup d'oeil aux données `env` :

```r
# Garder les lignes avec != 0 dans `env`
env <- env[row_sub,]

str(env)
# 'data.frame':	29 obs. of  11 variables:
#  $ das: num  0.3 2.2 10.2 18.5 21.5 ...
#  $ alt: int  934 932 914 854 849 846 841 752 617 483 ...
#  $ pen: num  48 3 3.7 3.2 2.3 3.2 6.6 1.2 9.9 4.1 ...
#  $ deb: num  0.84 1 1.8 2.53 2.64 2.86 4 4.8 10 19.9 ...
#  $ pH : num  7.9 8 8.3 8 8.1 7.9 8.1 8 7.7 8.1 ...
#  $ dur: int  45 40 52 72 84 60 88 90 82 96 ...
#  $ pho: num  0.01 0.02 0.05 0.1 0.38 0.2 0.07 0.3 0.06 0.3 ...
#  $ nit: num  0.2 0.2 0.22 0.21 0.52 0.15 0.15 0.82 0.75 1.6 ...
#  $ amm: num  0 0.1 0.05 0 0.2 0 0 0.12 0.01 0 ...
#  $ oxy: num  12.2 10.3 10.5 11 8 10.2 11.1 7.2 10 11.5 ...
#  $ dbo: num  2.7 1.9 3.5 1.3 6.2 5.3 2.2 5.2 4.3 2.7 ...
```

]

|Variable |Description|
|:--:|:--|
|das|Distance de la source [km] |
|alt|Altitude [m a.s.l.] |
|pen|Pente [par millier m.] |
|deb|Decharge min. moy. [m3s-1] |
|pH|pH de l'eau |
|dur|Conc. Ca (hardness) [mgL-1] |
|pho|Conc. K [mgL-1] |
|nit|Conc. N [mgL-1] |
|amn|Conc. NH₄⁺ [mgL-1] |
|oxy|Diss. d'oxygen [mgL-1] |
|dbo|Demande biol. d'oxygen [mgL-1] |

]

```r
names(env)
dim(env) # dimensions
summary(env) # sommaire
head(env)
```

]

---
# Mesures de (dis)similarité

La ressemblance des communautés est frequement évaluée sur la composition des espèces, sous la forme d'un tableau de données site par espèce `$Y_{m,n}$`.

On peut obtenir une matrice d'association `$A_{m,m}$` sous la forme de distances ou de dissimilarités par paires `$D_{m,m}$` (ou de similarités `$S_{m,m}$`), puis les analyser.

```r
dist(spe)
#            1         2         3         4         5         6         7
# 2   5.385165                                                            
# 3   7.416198  2.449490                                                  
# 4   7.874008  4.123106  3.000000                                        
# 5  10.816654 10.677078 10.862780  9.219544                              
# 6   7.348469  4.582576  4.123106  2.828427  8.185353                    
# 7   6.855655  2.449490  2.000000  3.605551 10.488088  3.605551          
# 9   7.810250  8.717798  9.380832  8.774964  9.380832  6.708204  8.485281
# 10  6.708204  5.099020  5.291503  5.000000  9.273618  3.605551  4.472136
# 11  4.472136  3.316625  5.000000  5.477226 10.246951  5.291503  4.795832
# 12  6.708204  3.162278  3.464102  4.582576 11.045361  4.795832  3.162278
# 13  7.071068  4.358899  5.196152  6.164414 11.532563  6.633250  5.385165
# 14  9.110434  6.324555  6.164414  6.557439 11.916375  7.000000  6.324555
# 15  9.899495  7.810250  7.681146  7.483315 10.816654  6.633250  6.855655
# 16 11.090537  9.899495  9.797959  9.327379 10.198039  8.426150  9.055385
# 17 10.630146  9.486833  9.486833  8.774964  9.797959  8.062258  9.055385
# 18  9.848858  9.591663  9.899495  8.660254  9.055385  7.810250  9.380832
# 19 11.704700 11.135529 11.045361 10.148892  9.380832  8.888194 10.770330
# 20 13.453624 14.212670 14.628739 13.453624 10.770330 12.206556 14.212670
# 21 14.456832 15.362291 15.811388 14.525839 11.489125 13.601471 15.556349
# 22 16.643317 17.663522 18.110770 16.763055 13.416408 15.779734 17.663522
# 23  3.872983  7.483315  9.055385  9.000000 10.583005  7.937254  8.485281
# 24  7.071068  9.539392 10.816654 10.583005 11.357817  9.486833 10.246951
# 25  5.291503  8.306624  9.643651  9.273618  9.746794  8.246211  9.110434
# 26 11.135529 12.609520 13.304135 12.409674 10.723805 11.224972 12.922848
# 27 15.362291 16.462078 16.881943 15.811388 13.076697 14.696938 16.583124
# 28 16.522712 17.549929 18.000000 16.941074 14.352700 15.905974 17.606817
# 29 18.761663 19.364917 19.621417 18.384776 15.524175 17.776389 19.364917
# 30 20.396078 21.377558 21.794495 20.542639 17.117243 19.849433 21.517435
#            9        10        11        12        13        14        15
# 2                                                                       
# 3                                                                       
# 4                                                                       
# 5                                                                       
# 6                                                                       
# 7                                                                       
# 9                                                                       
# 10  6.480741                                                            
# 11  7.549834  4.582576                                                  
# 12  8.717798  5.477226  3.872983                                        
# 13 10.049876  6.855655  4.000000  3.316625                              
# 14 10.583005  7.615773  6.244998  4.242641  3.316625                    
# 15  9.746794  6.855655  7.745967  6.403124  6.633250  5.000000          
# 16 10.862780  8.485281 10.049876  9.591663  9.746794  8.944272  5.916080
# 17  9.899495  8.246211  9.219544  9.273618  9.643651  9.380832  8.185353
# 18  8.602325  7.874008  8.774964  9.380832  9.949874  9.797959  8.660254
# 19  9.486833  9.273618 11.000000 11.313708 12.041595 11.832160 10.630146
# 20 11.489125 12.727922 13.527749 14.422205 15.000000 14.966630 13.674794
# 21 13.266499 14.142136 14.662878 15.684387 16.031220 16.062378 15.066519
# 22 15.033296 16.248077 16.941074 17.832555 18.303005 18.165902 16.763055
# 23  6.324555  6.633250  5.744563  8.366600  8.774964 10.392305 10.630146
# 24  7.549834  8.544004  8.124038 10.148892 10.583005 11.789826 11.747340
# 25  7.280110  7.000000  6.782330  9.110434  9.486833 10.723805 10.583005
# 26 10.148892 11.445523 11.747340 13.000000 13.490738 13.892444 13.190906
# 27 13.892444 15.264338 15.748016 16.703293 17.146428 17.117243 16.062378
# 28 14.899664 16.309506 16.822604 17.720045 18.193405 18.220867 17.175564
# 29 17.521415 18.303005 18.761663 19.416488 19.646883 19.313208 18.330303
# 30 19.467922 20.371549 20.736441 21.610183 21.863211 21.931712 21.023796
#           16        17        18        19        20        21        22
# 2                                                                       
# 3                                                                       
# 4                                                                       
# 5                                                                       
# 6                                                                       
# 7                                                                       
# 9                                                                       
# 10                                                                      
# 11                                                                      
# 12                                                                      
# 13                                                                      
# 14                                                                      
# 15                                                                      
# 16                                                                      
# 17  6.164414                                                            
# 18  7.615773  3.464102                                                  
# 19  9.165151  6.633250  6.000000                                        
# 20 12.649111  9.380832  7.615773  6.324555                              
# 21 13.928388 11.224972  9.380832  7.874008  3.741657                    
# 22 15.556349 13.416408 11.489125 10.583005  6.480741  4.000000          
# 23 11.489125 10.295630  9.055385 10.488088 11.916375 13.114877 15.362291
# 24 12.449900 10.723805  9.110434  9.746794 10.148892 11.269428 13.304135
# 25 11.180340 10.148892  8.660254  9.539392 10.630146 11.618950 13.964240
# 26 13.076697 10.816654  8.774964  7.810250  5.385165  5.385165  7.549834
# 27 15.329710 13.304135 11.532563 10.246951  6.244998  4.582576  5.000000
# 28 16.370706 14.352700 12.489996 11.489125  7.745967  6.000000  5.099020
# 29 17.058722 14.662878 13.076697 12.767145  9.000000  7.000000  5.567764
# 30 19.621417 17.117243 15.652476 15.000000 11.180340  9.000000  8.185353
#           23        24        25        26        27        28        29
# 2                                                                       
# 3                                                                       
# 4                                                                       
# 5                                                                       
# 6                                                                       
# 7                                                                       
# 9                                                                       
# 10                                                                      
# 11                                                                      
# 12                                                                      
# 13                                                                      
# 14                                                                      
# 15                                                                      
# 16                                                                      
# 17                                                                      
# 18                                                                      
# 19                                                                      
# 20                                                                      
# 21                                                                      
# 22                                                                      
# 23                                                                      
# 24  4.358899                                                            
# 25  3.000000  3.741657                                                  
# 26  9.433981  7.071068  8.000000                                        
# 27 14.035669 11.916375 12.649111  5.291503                              
# 28 15.231546 13.076697 13.964240  7.000000  3.872983                    
# 29 17.804494 15.874508 16.431677  9.899495  6.633250  5.567764          
# 30 19.416488 17.832555 18.055470 11.916375  9.165151  7.280110  6.480741
```

]

Par mesure de simplification, nous avons fait cela sans spécifier d'arguments.

Essayez `dist(spe)` de votre côté, notez ce que vous observez et voyez ce que les commandes ci-dessous vous montrent :

```r
class(dist(spe))
```

```r
str(dist(spe))
```

```r
as.matrix(dist(spe))
```

```r
dim(as.matrix(dist(spe)))
```

]

???
Si possible, faites cette partie avec les participants.

Montrez-leur comment est structuré un objet de classe dist. Soulignez qu'il s'agit d'une matrice triangulaire inférieure, où les distances entre les mêmes objets sont cachées.

---
# Standardisation

La standardisation des variables environnementales est cruciale car vous ne pouvez pas comparer les effets de variables ayant des unités différentes :

```r
## ?decostand
env.z <- decostand(env, method = "standardize")
```

Cela centre et réduit les variables pour rendre votre analyse en aval plus appropriée:

```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
```

---

# Mesures de (dis)similarité

Il existe trois groupes de coefficients de distance : _métriques_, _..._ , _..._ .

Le premier groupe est constitué de *métriques*, et ses coefficients satisfont aux propriétés suivantes :

.pull-left[
1. minimum 0 : si l'espèce `$a$` est égale à l'espèce `$b$`, alors `$D(a,b)=0$` ;
2. positivité : si `$a \neq b$`, alors `$D(a,b) > 0$` ;
3. la symétrie : `$D(a,b) = D(b,a)$` ;
4. inégalité triangulaire : `$D(a,b) + D(b,c) \geq D(a,c)$`. La somme de deux côtés d'un triangle tracé dans l'espace euclidien est égale ou supérieure au troisième côté.
]

Nous pouvons repérer toutes ces propriétés ci-dessous :

```r
as.matrix(dist(spe))[1:3, 1:3]
#          1        2        3
# 1 0.000000 5.385165 7.416198
# 2 5.385165 0.000000 2.449490
# 3 7.416198 2.449490 0.000000
```

]

???

Soulignez la nature triangulaire de la classe `dist`, le fait qu'elle peut être convertie en matrice, et les dimensions (en notant qu'elle a le même nombre d'espèces).

---
# Mesures de (dis)similarité

Il existe trois groupes de coefficients de distance : _métriques_, _semimétriques_, _..._ .

Le deuxième groupe est constitué de *semimétriques*, et ils violent la propriété d'_inégalité des triangles_ :

.pull-left[
1. minimum 0 : si l'espèce `$a$` est égale à l'espèce `$b$`, alors `$D(a,b)=0$` ;
2. positivité : si `$a \neq b$`, alors `$D(a,b) > 0$` ;
3. la symétrie : `$D(a,b) = D(b,a)$` ;
4. ~~inégalité triangulaire~~ : `${D(a,b) + D(b,c) \geq ou < D(a,c)}$`. La somme de deux côtés d'un triangle tracé dans l'espace euclidien n'est _pas_ égale ou supérieure au troisième côté.
]

]

???

Soulignez la nature triangulaire de la classe `dist`, le fait qu'elle peut être convertie en matrice, et les dimensions (en notant qu'elle a le même nombre d'espèces).

---
# Mesures de (dis)similarité

Il existe trois groupes de coefficients de distance : _métriques_, _semimétriques_, _nonmétriques_.

Le troisième groupe est constitué de coefficients *non métriques*, et ils violent la propriété de _positivité_ des distances métriques :

.pull-left[
1. minimum 0 : si l'espèce `$a$` est égale à l'espèce `$b$`, alors `$D(a,b)=0$` ;
2. ~~positivité:~~ si `$a \neq b$`, alors `$D(a,b) > ou < 0$` ;
3. la symétrie : `$D(a,b) = D(b,a)$` ;
4. ~~inégalité triangulaire~~ : `${D(a,b) + D(b,c) \geq ou < D(a,c)}$`. La somme de deux côtés d'un triangle tracé dans l'espace euclidien n'est _pas_ égale ou supérieure au troisième côté.
]

]

---

# Mesures de (dis)similarité : Distances euclidiennes

La mesure de distance métrique la plus courante est la _distance euclidienne_.

.pull-left4[
Elle est calculée à l'aide de la formule de Pythagore :
`$$D_{1} (x_1,x_2) = \sqrt{\sum_{j=1}^p(y_{1j} - y_{2j})^2}$$`

]

```r
spe.D.Euclid <- dist(x = spe,
 method = "euclidean")
```

Et, nous pouvons tester si une distance est euclidienne en utilisant :

```r
is.euclid(spe.D.Euclid)
# [1] TRUE
```

]

???

La figure sur le côté doit être terminée. Elle est censée représenter la distance euclidienne entre les sites x1 et x2 dans l'espace cartésien 2D. Les axes sont les descripteurs y1 et y2, respectivement. La flèche en bas représente la distance entre la position de y21 et celle de y11.

Vous pouvez l'utiliser pour vous rappeler le théorème de Pythagore que tout le monde apprend au lycée, où la longueur de l'hypoténuse et a et b désignent les deux longueurs des branches d'un triangle rectangle.

`method = "euclidean"` est le paramètre par défaut.

---

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

**Votre tour!** En utilisant la fonction `dist()`, calculez la matrice de distance euclidienne `$D_{hmm}$` pour les abondances des espèces par matrice de site `$Y_{hmm}$` ci-dessous :

```r
Y.hmm <- data.frame(
 y1 = c(0, 0, 1),
 y2 = c(4, 1, 0),
 y3 = c(8, 1, 0))
```
]

| Sites | `$y_1$` | `$y_2$` | `$y_3$` |
|:----: | :----:| :---: | :---: |
| `$s_1$` | 0 | 4 | 8 |
| `$s_2$` | 0 | 1 | 1 |
| `$s_3$` | 1 | 0 | 0 |

]

Après cela, examinez les chiffres, réfléchissez-y de manière critique et soyez prêts à discuter ! (_5 à 10 minutes_)

**Solution:**

```r
Y.hmm.DistEu <- dist(x = Y.hmm,
 method = "euclidean")
as.matrix(Y.hmm.DistEu)
```

]

**Résultat**:

```
#          1        2        3
# 1 0.000000 7.615773 9.000000
# 2 7.615773 0.000000 1.732051
# 3 9.000000 1.732051 0.000000
```

]

.center[
Suggestion : Examinez la composition et les distances entre les sites `$s_2$` et `$s_3$` et entre `$s_1$` et `$s_2$`.
]

???

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

**Solution:**

```r
Y.hmm
#   y1 y2 y3
# 1  0  4  8
# 2  0  1  1
# 3  1  0  0
```

]

**Résultat**:

```r
as.matrix(Y.hmm.DistEu)
#          1        2        3
# 1 0.000000 7.615773 9.000000
# 2 7.615773 0.000000 1.732051
# 3 9.000000 1.732051 0.000000
```

]

La distance euclidienne entre les sites `$s_2$` et `$s_3$`, qui n'ont aucune espèce en commun, est plus petite que la distance entre `$s_1$` et `$s_2$`, qui partagent les espèces `$y_2$` et `$y_3$` ( !).

D'un point de vue écologique, _cette évaluation de la relation entre les sites est problématique_.

Ce problème est comme le **paradoxe du double zéro**, _i.e._ les doubles zéros sont traités de la même manière que les doubles présences, de sorte que les doubles zéros réduisent la distance entre deux sites.

Les distances euclidiennes ( `$D_1$` ) ne devraient donc _pas_ être utilisées pour comparer des sites sur la base de l'abondance des espèces.

---

# Mesures de (dis)similarité : la distance des cordes

Orlóci (1967) a proposé la _distance de corde_ pour analyser la composition des communautés.

1\. Normaliser les données, _c.-à-d. _ mettre à l'échelle les vecteurs de site à la longueur 1 en divisant les abondances des espèces dans un échantillon donné par la somme à racine carrée des abondances carrées dans tous les échantillons, comme suit

`$$y'_{Uj}=y_{Uj}/\sum^s_{j=1}{y^2_{Uj}}$$`
]

`$$D_{3} (x_1,x_2) = \sqrt{\sum_{j=1}^p(y'_{1j} - y'_{2j})^2}$$`
]

---

# Mesures de (dis)similarité : la distance des cordes

Nous pouvons utiliser `vegan::vegdist()` pour cela :

```r
spe.D.Ch <- vegdist(spe,
 method = "chord")
as.matrix(spe.D.Ch)[1:3, 1:3]
# 1 2 3
# 1 0.0000000 0.7653669 0.9235374
# 2 0.7653669 0.0000000 0.2309609
# 3 0.9235374 0.2309609 0.0000000
```

Lorsque deux sites partagent les mêmes espèces dans les mêmes proportions du nombre d'individus, la valeur de `$D_3$` est de `$0$`, et lorsqu'aucune espèce n'est partagée, sa valeur est de `$\sqrt{2}$`.

---
# Mesures de (dis)similarité : la distance des cordes

Il existe de nombreuses autres distances métriques d'intérêt, que nous pouvons utiliser :

Que se passe-t-il si nous calculons les distances d'accord dans la même matrice site-espèce `$Y_{hmm}$` ?
.pull-left[

```r
Y.hmm
#   y1 y2 y3
# 1  0  4  8
# 2  0  1  1
# 3  1  0  0
```

```r
as.matrix(Y.hmm.DistEu)
#          1        2        3
# 1 0.000000 7.615773 9.000000
# 2 7.615773 0.000000 1.732051
# 3 9.000000 1.732051 0.000000
```

]

```r
Y.hmm.DistCh <- vegdist(Y.hmm,
 method = "chord")
```

```r
as.matrix(Y.hmm.DistCh)
#           1         2        3
# 1 0.0000000 0.3203645 1.414214
# 2 0.3203645 0.0000000 1.414214
# 3 1.4142136 1.4142136 0.000000
```
]

L'ajout d'un nombre quelconque de doubles zéros à une paire de sites ne modifie pas la valeur de `$D_3$`.

Par conséquent, les _distances de cordes_ peuvent être utilisées pour comparer des sites décrits par des abondances d'espèces !

---
# Mesures de (dis)similarité : Jaccard

Un autre coefficient d'association populaire est le _coefficient de similarité de Jaccard_ (1900).

Il n'est approprié que pour les **données binaires**, et son coefficient de distance est défini par la taille de l'intersection divisée par la taille de l'union des ensembles d'échantillons

`$$D_{7}(x_1,x_2)=1-{ {|x_1 \cap x_2|}\over{|x_1 \cup x_2|} } = 1-{ {|x_1 \cap x_2|}\over{|x_1| + |x_2| - |x_1 \cap x_2|} } = 1-\frac{a}{a+b+c}$$`
où,

- `$a$` est le nombre d'espèces partagées entre `$x_1$` et `$x_2$` qui sont codées `$1$` ;
- `$b$` est le nombre d'occurrences où l'on sait que `$x_1$` et `$x_2$` sont différents ;
- `$c$` est le nombre d'absences communes entre `$x_1$` et `$x_2$`, _c.-à-d._ toutes deux `$0$`.

| `$x_1,x_2$` | `$y_1$` | `$y_2$` | `$y_3$` | `$y_4$` | `$y_5$` |
|:----: | :----:| :---: | :---: | :---: | :---: |
| `$x_1$` | 0 | 1 | 0 | 1 | 0 |
| `$x_2$` | 0 | 1 | 1 | 1 | 1 |
]

]
.pull-right[
 
`$D_{7}(x_1,x_2) =$`

`$1-\frac{2}{2+2+1}=$`

`$0.6$`
]
]

---
# Mesures de (dis)similarité : Sørensen

Tous les paramètres du _coefficient de similarité de Jaccard_ ont le même poids.

`$$D_{7}(x_1,x_2)=1-\frac{a}{a+b+c}$$`

Cependant, vous pouvez envisager que la présence d'une espèce soit plus informative que son absence.

La distance correspondant au _coefficient de similarité de Sørensen_ (1948) donne du poids aux doubles présences :

`$$D_{13}(x_1,x_2)=1-\frac{2a}{2a+b+c}=\frac{b+c}{2a+b+c}$$`

où,
- `$a$` est le nombre d'espèces partagées entre `$x_1$` et `$x_2$` qui sont codées `$1$` ;
- `$b$` est le nombre d'occurrences où l'on sait que `$x_1$` et `$x_2$` sont différents ;
- `$c$` est le nombre d'absences communes entre `$x_1$` et `$x_2$`, _c.-à-d. _ toutes deux `$0$`.

---
# Mesures de (dis)similarité : Bray-Curtis

Le _coefficient de dissimilarité de Bray-Curtis_ est une version modifiée de l'indice de Sørensen et tient compte de l'abondance des espèces :

`$$D_{14}(x_1,x_2)=1 - \frac{2W}{A+B}$$`
où,
- `$W$` est la somme des plus faibles abondances de chaque espèce trouvées entre les sites `$x_1$` et `$x_2$` ;
- `$A$` est la somme de toutes les abondances dans `$x_1$` ; et,
- `$B$` est la somme de toutes les abondances dans `$x_2$`.
]

| `$x_1,x_2$` | `$y_1$` | `$y_2$` | `$y_3$` | `$y_4$` | `$y_5$` |
|:----: | :----:| :---: | :---: | :---: | :---: |
| `$x_1$` | _2_ | _1_ | _0_ | 5 | 2 |
| `$x_2$` | 5 | _1_ | 3 | _1_ | 1 |

Alors:
- `$W = 2 + 1 + 0 + 1 + 1 = 5$`
- `$A = 2 + 1 + 0 + 5 + 2 = 10$`
- `$B = 5 + 1 + 3 + 1 + 1 = 10$`

`$$D_{14}(x_1,x_2) = 1-\frac{2 \times 5}{10+10}=$$`
`$$D_{14}(x_1,x_2) = 0.5$$`
]

---
##### Mesures de (dis)similarité : Coefficients de Jaccard et de Sørensen

Dans `R`, vous pouvez utiliser la fonction `vegan::vegdist()` pour calculer les indices de Jaccard et de Sørensen :

```r
spe.D.Jac <- vegdist(spe,
 method = "jaccard",
 binary = TRUE)
```

]

```r
spe.D.Sor <- vegdist(spe,
 method = "bray",
 binary = TRUE)
```
]
 
> Comme la méthode de Jaccard et celle de Sørensen ne sont appropriées que pour les données de présence-absence, vous devez effectuer une transformation binaire des données d'abondance en utilisant `binary = TRUE` dans `vegdist()`.

##### Coefficient de dissimilarité de Bray-Curtis

Pour calculer le _coefficient de dissimilarité de Bray-Curtis_, qui peut tenir compte des abondances, vous devez définir `binary = FALSE`.

```r
spe.D.Bray <- vegdist(spe,
 method = "bray",
 binary = FALSE)
```

---
### Mesures de (dis)similarité : Distance de Mahalanobis

.pull-left[
Les variables sont souvent corrélées et l'utilisation d'une mesure de distance qui n'en tient pas compte peut conduire à des résultats erronés.

La distance de Mahalanobis tient compte de la corrélation entre les variables dans les données.

`$$D_{Mahal}^2 = (x - \mu)^T \Sigma^{-1} (x - \mu)$$`

`$D$` représente la distance de Mahalanobis, `$x$` est le point à partir duquel vous mesurez la distance, `$\mu$` est le vecteur moyen de la distribution que vous mesurez, et `$\Sigma$` est la matrice de covariance de la distribution.

La distance de Mahalanobis est également robuste aux changements dans la distribution des données (elle ne présuppose pas la normalité).
]

Dans `R`, nous pouvons calculer la distance de Mahalanobis comme suit :

```r
# Pour la reproductibilité
set.seed(123)
# Générer des données aléatoires à partir d'une distribution normale multivariée
data <- data.frame(x = rnorm(100, mean = 10, sd = 2),
 y = rnorm(100, mean = 5, sd = 1),
 z = rnorm(100, mean = 20, sd = 4))

# Calculer la distance de Mahalanobis
mahalanobis_dist <- mahalanobis(data, 
 center = colMeans(data), 
 cov = cov(data))
```

]

---
# Mesures de (dis)similarité : représentation

Nous pouvons créer des représentations graphiques des matrices d'association en utilisant la fonction `coldiss()` :

```r
coldiss(spe.D.Jac)
```

---
# Transformations pour les données de composition des communautés

Les communautés échantillonnées dans des conditions environnementales homogènes ou courtes peuvent avoir des compositions d'espèces avec peu de zéros, de sorte que les distances euclidiennes pourraient suffire à les caractériser.

Néanmoins, c'est rarement la réalité.

Les espèces peuvent être très fréquentes lorsque les conditions sont favorables, ou être absentes de nombreux sites. Parfois, cette asymétrie peut introduire des problèmes parasites dans nos analyses.

Nous pouvons alors être amenés à transformer nos données de composition pour les analyser de manière appropriée.

Dans `R`, nous pouvons compter sur `vegan::decostand()` pour de nombreux types de transformations.

Jetez un coup d'oeil à l'aide de cette fonction pour voir les options disponibles :

```
?decostand()
```

Voyons-en quelques-uns.

---
## Transformations : présence-absence

Nous pouvons changer l'argument `method` en `"pa"` pour transformer nos données d'abondance en données de présence-absence :

```r
spe[1:6, 1:6]
#   CHA TRU VAI LOC OMB BLA
# 1   0   3   0   0   0   0
# 2   0   5   4   3   0   0
# 3   0   5   5   5   0   0
# 4   0   4   5   5   0   0
# 5   0   2   3   2   0   0
# 6   0   3   4   5   0   0
```
]

Transformons les abondances `spe` en présence-absence :

```r
spe.pa <- decostand(spe, method = "pa")
spe.pa[1:6, 1:6]
# CHA TRU VAI LOC OMB BLA
# 1 0 1 0 0 0 0
# 2 0 1 1 1 0 0
# 3 0 1 1 1 0 0
# 4 0 1 1 1 0 0
# 5 0 1 1 1 0 0
# 6 0 1 1 1 0 0
```

---
## Transformations : profils d'espèces

Parfois, on souhaite éliminer les effets des unités très abondantes. Nous pouvons transformer les données en profils d'abondance relative des espèces grâce à l'équation suivante :

`$$y'_{ij} = \frac{y_{ij}}{y_{i+}}$$`

```r
spe[1:5, 1:6]
#   CHA TRU VAI LOC OMB BLA
# 1   0   3   0   0   0   0
# 2   0   5   4   3   0   0
# 3   0   5   5   5   0   0
# 4   0   4   5   5   0   0
# 5   0   2   3   2   0   0
```
]

Avec `decostand()`:

```r
spe.total <- decostand(spe,
 method = "total")
spe.total[1:5, 1:6]
# CHA TRU VAI LOC OMB BLA
# 1 0 1.00000000 0.00000000 0.00000000 0 0
# 2 0 0.41666667 0.33333333 0.25000000 0 0
# 3 0 0.31250000 0.31250000 0.31250000 0 0
# 4 0 0.19047619 0.23809524 0.23809524 0 0
# 5 0 0.05882353 0.08823529 0.05882353 0 0
```

???

yi+ est utilisé pour indiquer le nombre total de l'échantillon pour toutes les espèces j=1,...,m, pour le ième échantillon.

---
## Transformations : Hellinger

Nous pouvons prendre la racine carrée de la _transformation du profil de l'espèce_ et obtenir la _transformation de Hellinger_, qui a de très bonnes propriétés mathématiques et permet de réduire les effets des valeurs de `$y_{ij}$` qui sont extrêmement grandes.

`$$y'_{ij} = \sqrt{\frac{y_{ij}}{y_{i+}}}$$`

```r
spe[1:5, 1:6]
#   CHA TRU VAI LOC OMB BLA
# 1   0   3   0   0   0   0
# 2   0   5   4   3   0   0
# 3   0   5   5   5   0   0
# 4   0   4   5   5   0   0
# 5   0   2   3   2   0   0
```
]

Avec `decostand()`:

```r
spe.total <- decostand(spe,
 method = "hellinger")
spe.total[1:5, 1:6]
# CHA TRU VAI LOC OMB BLA
# 1 0 1.0000000 0.0000000 0.0000000 0 0
# 2 0 0.6454972 0.5773503 0.5000000 0 0
# 3 0 0.5590170 0.5590170 0.5590170 0 0
# 4 0 0.4364358 0.4879500 0.4879500 0 0
# 5 0 0.2425356 0.2970443 0.2425356 0 0
```

---
exclude: true
# Richesse totale en espèces

Visualisez le nombre d'espèces présentes sur chaque site :

```r
site.pre <- rowSums(spe > 0)
barplot(site.pre, main = "Species richness",
 xlab = "Sites",
 ylab = "Number of species",
 col = "grey ", las = 1)
```

---
exclude: true

# Visualisation des matrices de distance

```r
# the code for the coldiss() function is in the workshop script.
coldiss(spe.D.Bray)
```

---
class: inverse, middle, center

# _Clustering_ ou regroupement

---
# Regroupement

Une application des matrices d'association est le _regroupement hierarchique._

Le regroupement met en évidence les structures des données en partitionnant soit les objets, soit les descripteurs.

Un objectif des écologistes pourrait être de diviser un ensemble de sites en groupes en fonction de leurs conditions environnementales ou de la composition de leur communauté.

Ces résultats peuvent être représentés sous forme de dendrogrammes (diagrammes en forme d'arbre), qui décrivent la proximité des observations.

???
.center[![:scale 80%](images/cluster1_revised.jpg)]

---
# Méthodes de regroupement hiérarchique

Les algorithmes de regroupement s'appuient généralement sur cette série d'étapes :

1. Calculer une matrice d'association avec toutes les similitudes par paires entre tous les objets ;
2. Réunir les paires d'objets qui sont les plus similaires (ou dissemblables) ;
3. Recalculer la matrice de similarité pour ce groupe _par rapport à tous les objets restants_ ;
4. Répétez les étapes 2 et 3 jusqu'à ce que tous les objets soient réunis.

--

Parmi les nombreux algorithmes de clustering hiérarchique existants, nous allons explorer :

1. Le regroupement agglomératif à liaison unique ;
2. Liaison complète, regroupement agglomératif ;
3. Regroupement à variance minimale de Ward.

Dans `R`, nous utiliserons les fonctions `hclust()` et `agnes()` pour construire nos dendrogrammes.

???

Dans le cas d'un regroupement agglomératif par liaison simple (également appelé tri par
voisins), les objets les plus proches s'agglomèrent.
ce qui génère souvent de longues et fines grappes ou chaînes d'objets. À l'inverse, dans le clustering agglomératif à liens complets, un objet ne s'agglomère à un groupe que lorsqu'il est lié à l'élément le plus éloigné du groupe, qui le lie à son tour à tous les membres de ce groupe. Il formera de nombreux petits groupes séparés, et est plus approprié pour rechercher des contrastes, des discontinuités dans les données.
Le clustering à variance minimale de Ward diffère de ces deux méthodes en ce que
elle regroupe les objets en groupes en utilisant le critère des moindres carrés
(similaire aux modèles linéaires). Son dendogramme montre par défaut les distances au carré. Pour comparer avec d'autres méthodes, calculez d'abord la racine carrée des distances.

---

# Regroupement complet de liens

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

]

- Les objets sont répartis en petits groupes (1-2, 3-4, 5)

- Relier les petits groupes en utilisant la plus grande distance entre leurs éléments
* (1-3=0,15, 2-4=0,35, 2-3=0,6, choisir 0,6 pour relier le groupe 1-2 et 3-4)

![](images/compleClust2.png)
]

---
# Comparison

Créez une matrice de distance à partir des données sur les rivières du Doubs transformées par Hellinger et calculez le regroupement par liaison simple :

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

---
# Comparaison entre methodes

```r
spe.dhe1 <- vegdist(spe.hel, method = "euclidean")
spe.dhe1.complete <- hclust(spe.dhe1, method = "complete")
plot(spe.dhe1.single, main="Single linkage clustering", hang =-1)
plot(spe.dhe1.complete, main="Complete linkage clustering", hang=-1)
```

**Single linkage:**

Des chaînes d'objets se produisent  (e.g. 19, 29, 30, 26)
]

**Complete linkage:**
Les groupes contrastés sont formés d'objets se produisent
]

???
![](images/comparison.png)

---
# Méthode de variance minimale de Ward

- À chaque étape, la paire de clusters qui fusionne est celle qui entraîne l'augmentation minimale de la somme totale des carrés à l'intérieur du groupe.

- Les groupes générés par cette méthode ont tendance à être plus sphériques et à contenir un nombre similaire d'objets.
]

---
# Méthode de variance minimale de Ward

Calculez le regroupement à variance minimale de Ward et tracez 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 0
```

> L'algorithme de clustering de Ward ne peut être appliqué qu'avec des distances euclidiennes. Les distances de Bray-Curtis ne sont pas euclidiennes, elles doivent donc être racinées au carré avant.

---
# Questions importantes à aborder

- Le choix de la méthode "appropriée" dépend de votre objectif ;
- Voulez-vous mettre en évidence des gradients ou des contrastes ?
- Il ne s'agit *pas* d'une méthode statistique et vous devrez peut-être vous appuyer sur des méthodes complémentaires pour étayer vos hypothèses.

Cependant, le clustering nous permet de :
- Déterminer le nombre optimal de clusters interprétables ;
- Calculer les statistiques de clustering ;
- Combiner le clustering à l'ordination pour distinguer des groupes de sites.

---

#### Et maintenant ?

Alors que l'**analyse de regroupements** remonte les *discontinuités* dans un ensemble de données, l'**ordination** extrait les principales tendances sous la forme d'axes continus.

Nous allons examiner maintenant trois types de **méthodes d'ordination sans contrainte**...

--
 
.center[*Si personne ne s'exprime, choisissez une personne "volontaire" !*]

**Les ordinations non contraintes** évaluent les relations chez un seul ensemble de variables. *Aucune tentative* n'est faite pour définir la relation entre un ensemble de variables indépendantes et une ou plusieurs variables dépendantes.

L'interprétation des effets potentiels d'autres facteurs sur les patrons observés ne peut se faire qu'indirectement, car ces facteurs ne sont *pas* explicitement inclus dans les analyses !

Ici, nous allons explorer :
.small[
.center[
.pull-left2[
**A**nalyse des **C**omposants **P**rincipales

**A**nalyse des **Co**ordonnées **P**rincipales
]

---
### Mais d'abord, faisons un *récapitulatif*...

Nous avons déjà compris la signification de **variance** et **moyenne**, et comment les calculer :

Elles sont très utiles pour comprendre le *centre* et la *dispersion* d'une variable ou d'une dimension donnée.

Néanmoins, nous sommes souvent intéressés **par plus d'une dimension**, et nous voulons mesurer dans quelle mesure chaque dimension varie de la moyenne *par rapport aux autres*.

La **covariance** est une telle mesure, qui peut décrire comment **deux dimensions co-varient** :

---

### Mais d'abord, faisons un *récapitulatif*...

Intuitivement, nous pouvons mesurer la **covariance** entre plus de deux variables. Disons, entre les variables `$x$`, `$y$`, et `$z$` :

`$$cov_{z,y}=\frac{\sum_{i=1}^{N}(z_{i}-\bar{z})(y_{i}-\bar{y})}{N-1}$$`
Nous pouvons représenter ces calculs dans une **matrice de covariance** :

`$${C(x, y, z)} = \left[ \begin{array}{ccc}
cov_{x,x} & cov_{y,x} & cov_{z,x} \\
cov_{x,y} & cov_{y,y} & cov_{z,y} \\
cov_{x,z} & cov_{y,z} & cov_{z,z}
\end{array} \right]$$`

**QUIZ**: *Que sont les diagonales ?* *Et, que se passe-t-il si les variables sont indépendantes ?*

---

##### Toujours en train de *récapituler*...

Si, `$cov_{x,y}=\frac{\sum_{i=1}^{N}(x_{i}-\bar{x})(y_{i}-\bar{y})}{N-1}$`, alors:

`$$cov_{x,x}=\frac{\sum_{i=1}^{N}(x_{i}-\bar{x})(x_{i}-\bar{x})}{N} = \frac{\sum_{i=1}^{n}(x_i - \bar{x})^2} {N} = var_x$$`
--

Afin que :

`$${C(x, y, z)} = \left[ \begin{array}{ccc}
cov_{x,x} & cov_{y,x} & cov_{z,x} \\
cov_{x,y} & cov_{y,y} & cov_{z,y} \\
cov_{x,z} & cov_{y,z} & cov_{z,z}
\end{array} \right] = \left[ \begin{array}{ccc}
var_{x} & cov_{y,x} & cov_{z,x} \\
cov_{x,y} & var_{y} & cov_{z,y} \\
cov_{x,z} & cov_{y,z} & var_{z}
\end{array} \right]$$`

---

##### Toujours en train de *récapituler*...

```r
x <- rnorm(5000, mean = 0, sd = 1)
y <- rnorm(5000, mean = 0, sd = 1)
z <- rnorm(5000, mean = 0, sd = 1)

xyz <- data.frame(x, y, z)

GGally::ggpairs(xyz)
```

<img src="workshop09-pres-fr_files/figure-html/xyz-norm-1.png" width="360" style="display: block; margin: auto;" />
]

```r
cov(xyz)
```

```
#          x        y       z
# x  0.98469 -0.00787 0.00880
# y -0.00787  1.00997 0.00286
# z  0.00880  0.00286 1.01408
```

Si les variables sont parfaitement indépendantes (ou non corrélées), la matrice de covariance `$C(x, y, z)$` est :

`$${C(x, y, z)} =  \left[ \begin{array}{ccc}
var_{x} & 0 & 0 \\
0 & var_{y} & 0 \\
0 & 0 & var_{z}
\end{array} \right]$$`

*i*.*e*. une covariance plus proche de `$1$` signifie que les variables sont *colinéaires*.

Et ici, `$var_{x} = var_{y} = var_{z} = 1$`.

]

---
## Transformations linéaires

Nous sommes souvent intéressés à l'observation des variables d'autres *formes*.

Pour cela, nous créons une nouvelle variable, disons `$x_{new}$`, en utilisant des constantes pour modifier la variable originale `$x$`. Par exemple :

.pull-left[
On peut transformer des distances mesurées en kilomètres `$d_{km}$` en miles :
`$$d_{mi} = 0.62 \times d_{km}$$`
]

`$$T_{C} = 0.5556\times T_{Fahr} - 17.778$$`
]

Ces exemples sont des **transformations linéaires** car les variables transformées sont liées linéairement aux variables d'origine et les formes de la distribution ne sont pas modifiées.

Deux types de transformations sont très importantes pour nous :

`$$x' = x_i - \bar{x}$$`
]

`$$x'' = \frac{x_i}{\sigma_x}$$`
]

???

Le centrage est essentiellement une technique où la moyenne des variables indépendantes est soustraite de toutes les valeurs. Cela signifie que toutes les variables indépendantes ont une moyenne nulle. La réduction est complementaire au centrage. Les variables prédictives sont divisées par leur écart type.

Le présentateur doit mentionner ici que le centrage amène la moyenne à zéro et la mise à l'échelle amène l'écart type à une unité. Ils doivent également mentionner que les variables deviennent comparables lors de la mise à l'échelle, car leur unité est perdue.

---

### Décomposition des matrices en éléments propres

**Les matrices carrées**, telles que la **matrice de covariance**, peuvent être décomposées en *valeurs propres* et en *vecteurs propres*.

Pour une matrice carrée, `$A_{n \times n}$`, un vecteur `$v$` est un *vecteur propre* de `$A$` s'il y a un *scalaire*, `$\lambda$`, pour lequel :

.center[
`$A_{n \times n} v_{n \times 1} = \lambda  v_{n \times 1}$`, or `$\left(\begin{matrix}a_{11}&\cdots&a_{1n}\\\vdots&\ddots&\vdots\\a_{1n}&\cdots&a_{nn}\\\end{matrix}\right)\left(\begin{matrix}v_1\\\vdots\\v_n\\\end{matrix}\right)=\lambda\left(\begin{matrix}v_1\\\vdots\\v_n\\\end{matrix}\right)$`
]

la valeur de `$\lambda$` étant la *valeur propre* correspondante.

c.-à-d. que la matrice `$A$` *étire* effectivement le vecteur propre `$v$` par la quantité spécifiée par la valeur propre (*scalaire*) `$\lambda$`.

Un *vecteur propre* est un vecteur dont la direction reste inchangée lorsqu'on lui applique une **transformation linéaire**.

---

### Décomposition des matrices en éléments propres

Représentons cela avec cet exemple simple.

Nous pouvons transformer un carré en un parallélogramme à l'aide d'une **transformation affine** à axe unique.

.pull-left[
Soit `$S$` le carré de vertices `$(0,0),\,(1,0),\,(1,1),\,(1,0)$` qui sera transformé par cisaillement en parallélogramme `$P$` de vertices `$(0,0),\,(1,0),\,(1,1.57),\,(1,0.57)$`.

Nous pouvons voir qu'après la transformation linéaire, la flèche violette n'a pas changé de direction, c.-à-d. qu'elle est un *vecteur propre* de `$S$`.

En revanche, la flèche rouge a changé de direction, et *n'est donc pas* un *vecteur propre* de `$S$`.
]

.pull-right[
<img src="workshop09-pres-fr_files/figure-html/unnamed-chunk-51-1.png" width="360" style="display: block; margin: auto;" />
]

---

##### Décomposition des matrices en éléments propres : implications

### Orthogonalité

Une propriété *fabuleuse* et *simple* des matrices *symétriques* que nous pouvons expliquer ici !

.pull-left[
Supposons que `$x$` soit un vecteur propre de `$A$` correspondant à la valeur propre `$λ_1$` et `$y$` un vecteur propre de `$A$` correspondant à la valeur propre `$λ_2$`, avec `$λ_1≠λ_2$`.

`$$Ax=\lambda_1x \\
Ay=\lambda_2y$$`

Multiplions chacun d'eux par l'autre *vecteur propre* transposé.
`$$y^{\intercal}Ax=\lambda_1y^{\intercal}x \\ x^{\intercal}A^{\intercal}y=\lambda_2x^{\intercal}y$$`

]

Maintenant, soustrayons la deuxième équation de la première et utilisons la commutativité du produit scalaire :

`$y^{\intercal}Ax-x^{\intercal}A^{\intercal}y=\lambda_1y^{\intercal}x - \lambda_2x^{\intercal}y \\ 0 = (\lambda_1 - \lambda_2)y^{\intercal}x$`

Parce que nous savons que `$\lambda_1-\lambda_2\neq 0$`, alors
`$y^{\intercal}x = 0$`, *c-*.*à*.*-d*., `$\mathbf{x}\perp\mathbf{y}$`, *c-*.*à*.*-d*. sont **orthogonaux** !

*Alors, que révèle la décomposition de la matrice de variance-covariance en elements propres ?*
]

???

L'explication de cette partie est très utile et assez simple, afin que chacun puisse comprendre ce qu'est l'orthogonalité. Il s'agit d'opérations simples d'équations et de soustractions.

---
##### Décomposition des matrices en éléments propres : implications

### Maximisation

.pull-left[
Si `$v_i' v_i = 1$`, alors `$Av_i=\lambda_iv_i$` peut être écrite comme :
`$$v_i' A v_i = \lambda_i$$`
En effet, `$v' A v$`
est la variance d'une combinaison linéaire avec des poids en `$v$`, *i*.*e*. `$\text{Var}(v_i'\,A)=v_i'\,\text{Var}(A)\,v_i$`.

*Donc, on peut relier les points !*

]

.pull-right[
Rappelez-vous que les *valeurs propres* de notre *matrice de variance-covariance* `$A$` sont directement liées à la variance !

Pour trouver un vecteur `$v$` qui maximise la variance, `$v' A v$`, il suffit de choisir le *vecteur propre* correspondant à la plus grande *valeur propre* `$\lambda_i$` !

De sorte que la variance maximale soit de `$\lambda_1$` !
]

La *variance expliquée* de chaque *vecteur propre* obéit à l'ordre : `$\lambda_1 > \lambda_2 > \dots > \lambda_k$`.

Cela nous permet de condenser un plus grand nombre de variables originales en un ensemble plus petit de vecteurs sélectionnés avec une perte minimale d'informations (c.-à-d. **une réduction de la dimensionnalité**).

---
# Méthodes d'ordination sans contrainte

Ceci est un bon point de départ pour nous mettre sur la voie des **méthodes d'ordination sans contrainte** que nous allons étudier aujourd'hui !

.pull-left[
Ils nous permettent de :
- Évaluer les relations *dans* un ensemble de variables (espèces ou variables environnementales);
- Trouver les composantes clés de la variation entre les échantillons, les sites, les espèces ;
- Réduire le nombre de dimensions dans les données multivariées tout en limitant la perte substantielle d'informations ;
- Créer de nouvelles variables à utiliser dans des analyses ultérieures.
]

1. **P**rincipal **C**omponent **A**nalysis;

2. **P**rincipal **Co**ordinate **A**nalysis;

3. **N**on-Metric **M**ulti**d**imensional **S**caling;
]
---
# Analyse en Composantes Principales

.small[
L'analyse en composantes principales (ACP) est une technique de réduction de la dimensionnalité *linéaire*, c.-à-d. qu'elle réduit les données fortement corrélées.

En bref, l'ACP transforme *linéairement* des variables originales vers des nouvelles variables contenant des **composantes principales**, qui expliquent la plupart de la variance dans l'ensemble de données, c.-à-d. qui maximisent la séparation entre les données.
]

--
.pull-left[

`$$Z_p = ∑_{j=1}^p ϕ_j * X_j$$`
]
.small[
où,
1. `$Z_p$` est la composante principale `$p$` ;
2. `$ϕ_j$` est le vecteur de charge comprenant les `$j$` charges pour la composante principale `$p$`, c.-à-d. les coefficients de la combinaison linéaire des variables originales à partir desquelles les composantes principales sont construites ;
3. `$X_j$` est le prédicteur normalisé, c.-à-d. avec une moyenne égale à zéro et un écart-type égal à un.
]

???

c.-à-d. que la reconstruction des données peut être donnée par une simple combinaison linéaire des composantes.

Les futures versions de ce document devraient inclure une figure (telle qu'une carte thermique) avec les données normalisées originales sur le côté gauche, et sur le côté droit, une carte thermique des chargements fois une carte thermique des composants (qui peut être considérée comme égale à la somme des matrices de rang un) étant égale à la carte thermique des données reconstruites.

---
# Analyse en Composantes Principales

L'ACP peut être calculée d'au moins *quatre* façons différentes.

Pour des raisons de simplicité, nous nous concentrerons ici sur la façon d'obtenir des composantes principales à partir d'une matrice de corrélation.

Nous apprendrons à le faire à partir de rien, puis à utiliser les paquets `R` pour calculer les composantes principales.

---

### Analyse en composantes principales : étape par étape

.pull-left3[
1. Point de départ : une matrice `$Y$` de `$n$` observations et `$p$` variables continues normalement distribuées ;
]

```r
data(varechem)

str(varechem)
# 'data.frame':	24 obs. of  14 variables:
#  $ N       : num  19.8 13.4 20.2 20.6 23.8 22.8 26.6 24.2 29.8 28.1 ...
#  $ P       : num  42.1 39.1 67.7 60.8 54.5 40.9 36.7 31 73.5 40.5 ...
#  $ K       : num  140 167 207 234 181 ...
#  $ Ca      : num  519 357 973 834 777 ...
#  $ Mg      : num  90 70.7 209.1 127.2 125.8 ...
#  $ S       : num  32.3 35.2 58.1 40.7 39.5 40.8 33.8 27.1 42.5 60.2 ...
#  $ Al      : num  39 88.1 138 15.4 24.2 ...
#  $ Fe      : num  40.9 39 35.4 4.4 3 ...
#  $ Mn      : num  58.1 52.4 32.1 132 50.1 ...
#  $ Zn      : num  4.5 5.4 16.8 10.7 6.6 9.1 7.4 5.2 9.3 9.1 ...
#  $ Mo      : num  0.3 0.3 0.8 0.2 0.3 0.4 0.3 0.3 0.3 0.5 ...
#  $ Baresoil: num  43.9 23.6 21.2 18.7 46 40.5 23 29.8 17.6 29.9 ...
#  $ Humdepth: num  2.2 2.2 2 2.9 3 3.8 2.8 2 3 2.2 ...
#  $ pH      : num  2.7 2.8 3 2.8 2.7 2.7 2.8 2.8 2.8 2.8 ...
```

]

---

#### Analyse en composantes principales : étape par étape

.pull-left3[
1. Point de départ : une matrice `$Y$` de `$n$` observations et `$p$` variables continues normalement distribuées ;
]

```r
data(varechem)

# Étape 1
Y <- varechem[, 1:2]

head(Y)
#       N    P
# 18 19.8 42.1
# 15 13.4 39.1
# 24 20.2 67.7
# 27 20.6 60.8
# 23 23.8 54.5
# 19 22.8 40.9
```

]

---
#### Analyse en composantes principales : étape par étape

.pull-left3[
1. Point de départ : une matrice `$Y$` de `$n$` observations et `$p$` variables continues normalement distribuées ;

2. Normalisation des observations, comme dans `$Y_{std} = \frac{y_i - \bar{y}}{\sigma_y}$` ; ce qui revient à centrer, comme dans `$y_c = [y_i - \bar{y}]$`, puis à mettre à l'échelle, comme dans `$y_s = \frac{y_i}{\sigma_y}$` ;
]

```r
data(varechem)

# Étape 1
Y <- varechem[, 1:2]

# Étape 2
Y_std <- as.matrix(scale(Y))

head(Y_std)
#              N          P
# 18 -0.46731082 -0.1993234
# 15 -1.62503567 -0.4000407
# 24 -0.39495301  1.5134640
# 27 -0.32259521  1.0518143
# 23  0.25626722  0.6303080
# 19  0.07537271 -0.2796103

round(apply(Y_std, 2, mean))
# N P 
# 0 0
round(apply(Y_std, 2, sd))
# N P 
# 1 1
```

]

---

#### Analyse en composantes principales : étape par étape

.pull-left3[
1. Point de départ : une matrice `$Y$` de `$n$` observations et `$p$` variables continues normalement distribuées ;

3. Calculer la matrice de variance-covariance `$R = cov(Y_{std})$` ;
]

```r
data(varechem)

# Étape 1
Y <- varechem[, 1:2]

# Étape 2
Y_std <- as.matrix(scale(Y))

# Étape 3
(Y_R <- cov(Y_std))
# N P
# N 1.0000000 -0.2511603
# P -0.2511603 1.0000000
```
]

---

#### Analyse en composantes principales : étape par étape

.pull-left3[
1. Point de départ : une matrice `$Y$` de `$n$` observations et `$p$` variables continues normalement distribuées ;

3. Calculer la matrice de variance-covariance `$R = cov(Y_{std})$` ;

4. Effectuer la décomposition de la matrice de covariance pour obtenir la matrice `$U$` des vecteurs propres, contenant les *composantes principales* ;
]

```r
data(varechem)

# Étape 1
Y <- varechem[, 1:2]

# Étape 2
Y_std <- as.matrix(scale(Y))

# Étape 3
Y_R <- cov(Y_std)

# Étape 4
(Eigenvalues <- eigen(Y_R)$values)
# [1] 1.2511603 0.7488397

(Eigenvectors <- eigen(Y_R)$vectors)
# [,1] [,2]
# [1,] -0.7071068 -0.7071068
# [2,] 0.7071068 -0.7071068
```
]

---

### Analyse en composantes principales : étape par étape

Les *vecteurs propres* sont ici les **composantes principales**, et comme nous l'avons vu, à chaque *vecteur propre* correspond une *valeur propre*.

Nous pouvons représenter les distances des observations au premier vecteur propre (`PC1`, en rouge).

La première composante principale est dessinée de façon à ce que la variation des valeurs le long de sa ligne soit maximale.

Les flèches sur les composantes principales sont obtenues en multipliant leurs *valeurs propres* par les *vecteurs propres*.
]

.pull-right[
<img src="workshop09-pres-fr_files/figure-html/unnamed-chunk-57-1.png" width="360" style="display: block; margin: auto;" />
]

???

Dans cette première représentation, nous pouvons observer qu'une première direction (ou la première composante linéaire) est tracée en essayant de maximiser la variance des données.

Les participants peuvent vous demander en quoi cela est différent d'une régression linéaire. L'une des principales différences réside dans la façon dont les carrés d'erreurs sont minimisés perpendiculairement à la ligne droite (90 degrés, ce qui la rend orthogonale), alors que dans la régression linéaire, les carrés d'erreurs sont minimisés dans la direction des ordonnées.

---

### Analyse en composantes principales : étape par étape

Les *vecteurs propres* sont ici les **composantes principales**, et comme nous l'avons vu, à chaque *vecteur propre* correspond une *valeur propre*.

.pull-left[
Nous pouvons alors représenter les distances des observations au deuxième vecteur propre (`PC2`, en orange).

La deuxième composante principale est également dessinée en maximisant la variance des données.

Notez comment les composantes principales sont orthogonales !
]

]

--
.pull-left[
*Nous avons représenté les vecteurs propres, c.-à-d. les composantes principales !*

*Mais, quelle est l'utilité des valeurs propres ?*
]

???

Ici, la deuxième direction (ou la deuxième composante linéaire) est tracée de manière à ce que la variance des données soit maximisée par rapport à cette deuxième composante.

---

#### Analyse en composantes principales : étape par étape

Nous avons vu que les *valeurs propres* représentent la magnitude (la variance) des composantes principales.

.pull-left[
En fait, la somme de toutes les *valeurs propres* est égale à la somme des variances, qui sont représentées sur la diagonale de la matrice de variance-covariance.
]

```r
sum(diag(cov(Y_std)))
# [1] 2
sum(eigen(cov(Y_std))$values)
# [1] 2
```
]

Intuitivement, on peut obtenir l'influence relative de chaque *vecteur propre* `$v_{k}$` (ou `$\text{PC}_{k}$`) en divisant leurs valeurs par la somme de toutes les *valeurs propres*.

`$$\text{Variance expliquée de}~v_{k} = \frac{\lambda_{v_k}}{\sum^p_{i=1}{\lambda_{v}}}$$`

En faisant cela, nous pouvons dire que le `$\text{PC}1$` explique 63% de la variance dans les données, tandis que `$\text{PC}2$` explique 37% de la variance.

Enfin, nous pouvons procéder à la dernière étape de notre analyse des composantes principales !

---

#### Analyse en composantes principales : étape par étape

.pull-left3[
1. Point de départ : une matrice `$Y$` de `$n$` observations et `$p$` variables continues normalement distribuées ;

3. Calculer la matrice de variance-covariance `$R = cov(Y_{std})$` ;

4. Effectuer la décomposition de la matrice de covariance pour obtenir la matrice `$U$` des vecteurs propres, contenant les *composantes principales* ;

5. Obtenir la *matrice des coordonnées* `$F$` en multipliant `$U$` par la matrice normalisée `$Y_{std}$`.
]

```r
# Étape 1
Y <- varechem[, 1:2]

# Étape 2
Y_std <- as.matrix(scale(Y))

# Étape 3
Y_R <- cov(Y_std)

# Étape 4
Eigenvalues <- eigen(Y_R)$values
Eigenvectors <- eigen(Y_R)$vectors

# Étape 5
F_PrComps <- Y_std %*% Eigenvectors
head(F_PrComps)
# [,1] [,2]
# 18 0.1894957 0.4713816
# 15 0.8662023 1.4319452
# 24 1.3494546 -0.7909067
# 27 0.9718543 -0.5156358
# 23 0.2644868 -0.6269034
# 19 -0.2510109 0.1444178
```
]

---

### Analyse en composantes principales : étape par étape

La *matrice des coordonnées*, `$F$`, (objet `F_PrComps`) permet de *rotationner* le nouvel espace de données, afin qu'il soit représenté par rapport aux composantes principales.

<img src="workshop09-pres-fr_files/figure-html/unnamed-chunk-61-1.png" width="252" style="display: block; margin: auto;" />
]

<img src="workshop09-pres-fr_files/figure-html/unnamed-chunk-62-1.png" width="252" style="display: block; margin: auto;" />
]

???
Les titres des axes ne se font pas imprimer. Je les ai inclus en haut de chaque graphique, mais ce problème n'est pas résolu.

Le présentateur ou la présentatrice devrait souligner la rotation et parler de ce que sont les points. Vous pouvez survoler les points et leur montrer quelle était la position des points dans les graphiques " sans rotation ", et maintenant dans les graphiques " avec rotation ", en soulignant que maintenant, les " nouveaux axes " sont PC1 et PC2.

Cette compréhension sera utile lorsque les participants utiliseront les fonctions PCA implémentées dans R.

---

### Analyse en composantes principales : étape par étape

L'ACP peut aussi être calculée en utilisant les fonctions `stats::prcomp()`, `stats::princomp()`, `vegan::rda()`, et `ade4::dudi.pca()`.

```r
data(varechem)

Y <- varechem[, 1:2]
Y_std <- as.matrix(scale(Y))
Y_R <- cov(Y_std)

Eigenvalues <- eigen(Y_R)$values
Eigenvectors <- eigen(Y_R)$vectors

F_PrComps <- Y_std %*% Eigenvectors

head(F_PrComps)
#          [,1]       [,2]
# 18  0.1894957  0.4713816
# 15  0.8662023  1.4319452
# 24  1.3494546 -0.7909067
# 27  0.9718543 -0.5156358
# 23  0.2644868 -0.6269034
# 19 -0.2510109  0.1444178
```
]

```r
PCA_prcomp <- prcomp(Y,
 center = TRUE,
 scale = TRUE)

# or PCA_prcomp <- prcomp(Y_std)

head(PCA_prcomp$x)
#           PC1        PC2
# 18 -0.1894957 -0.4713816
# 15 -0.8662023 -1.4319452
# 24 -1.3494546  0.7909067
# 27 -0.9718543  0.5156358
# 23 -0.2644868  0.6269034
# 19  0.2510109 -0.1444178
```
]

---

### Analyse en composantes principales : étape par étape

L'ACP peut aussi être calculée en utilisant les fonctions `stats::prcomp()`, `stats::princomp()`, `vegan::rda()`, et `ade4::dudi.pca()`.

```r
data(varechem)

Y <- varechem[, 1:2]
Y_std <- as.matrix(scale(Y))
Y_R <- cov(Y_std)

Eigenvalues <- eigen(Y_R)$values
Eigenvectors <- eigen(Y_R)$vectors

F_PrComps <- Y_std %*% Eigenvectors

```r
PCA_princomp <- princomp(Y_std)

head(PCA_princomp$scores)
#        Comp.1     Comp.2
# 18 -0.1894957 -0.4713816
# 15 -0.8662023 -1.4319452
# 24 -1.3494546  0.7909067
# 27 -0.9718543  0.5156358
# 23 -0.2644868  0.6269034
# 19  0.2510109 -0.1444178
```
]

---

### Analyse en composantes principales : étape par étape

L'ACP peut aussi être calculée en utilisant les fonctions `stats::prcomp()`, `stats::princomp()`, `vegan::rda()`, et `ade4::dudi.pca()`.

```r
data(varechem)

Y <- varechem[, 1:2]
Y_std <- as.matrix(scale(Y))
Y_R <- cov(Y_std)

Eigenvalues <- eigen(Y_R)$values
Eigenvectors <- eigen(Y_R)$vectors

F_PrComps <- Y_std %*% Eigenvectors

```r
PCA_vegan_rda <- rda(Y_std)

scores(PCA_vegan_rda,
       display = "sites",
       scaling = 1,
       choices = seq_len(PCA_vegan_rda$CA$rank),
       const = sqrt(PCA_vegan_rda$tot.chi * (nrow(PCA_vegan_rda$CA$u) - 1)))[1:5, ]
#           PC1        PC2
# 18 -0.1894957 -0.4713816
# 15 -0.8662023 -1.4319452
# 24 -1.3494546  0.7909067
# 27 -0.9718543  0.5156358
# 23 -0.2644868  0.6269034
```

`vegan::rda()` utilise des mises à l'échelle alternatives. Vous pouvez étudier la `vignette("decision-vegan")`.
]

???

Dites aux participants que le nom `rda` fait référence à un type différent de technique d'ordination sous contrainte, mais que si nous exécutons `rda()` avec une seule variable, il exécutera une ACP.

---

# Analyse en composantes principales

Nous avons implémenté l'ACP sur un ensemble de données à deux variables pour plus de simplicité.

Avançons et appliquons-la à notre jeu de données sur les espèces de poissons.

Pour cela, nous allons utiliser la fonction `vegan::rda()` sur les données de poissons *transformées par Hellinger* et résumer les résultats :

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

# summary(spe.h.pca)
```

---

# Analyse en composantes principales

.pull-left[
Les premières lignes de `summary.rda()` nous renseignent sur la *variance totale* et la *variance sans contrainte* de notre modèle.
]

```
# [1] "Partitioning of variance:"        "              Inertia Proportion"
# [3] "Total          0.5025          1" "Unconstrained  0.5025          1"
```
]

```
# [1] "Importance of components:"                                                     
# [2] "                         PC1     PC2     PC3     PC4     PC5     PC6     PC7"  
# [3] "Eigenvalue            0.2580 0.06424 0.04632 0.03850 0.02197 0.01675 0.01472"  
# [4] "Proportion Explained  0.5133 0.12784 0.09218 0.07662 0.04371 0.03334 0.02930"  
# [5] "Cumulative Proportion 0.5133 0.64118 0.73337 0.80999 0.85370 0.88704 0.91634"  
# [6] "                          PC14     PC15     PC16      PC17      PC18      PC19"
# [7] "Eigenvalue            0.001835 0.001455 0.001118 0.0008309 0.0005415 0.0004755"
# [8] "Proportion Explained  0.003651 0.002895 0.002225 0.0016535 0.0010776 0.0009463"
# [9] "Cumulative Proportion 0.988888 0.991783 0.994008 0.9956612 0.9967389 0.9976852"
```

]

Viennent ensuite les *valeurs propres*, et leur contribution à la variance.

En fait, si nous additionnons toutes nos *valeurs propres*, nous obtiendrons la quantité de variance sans contrainte expliquée par l'analyse !

```r
sum(spe.h.pca$CA$eig)
# [1] 0.5025103
```

]

???

Puisque nous n'avons pas contraint notre ordination, la variance proportionnelle non contrainte est égale à la variance totale.

Prenez un moment pour expliquer la proportion expliquée, et montrez que la proportion cumulée sera égale à 1 au 27ème PC.

---

# Analyse en composantes principales

Les informations suivantes sont liées à la *mise à l'échelle* (`scaling`), aux *scores d'espèces*, et aux *scores de sites*.

```
#  [1] "Eigenvalue            0.0003680 0.0002765 0.0002253 0.0001429 7.618e-05"
#  [2] "Proportion Explained  0.0007324 0.0005503 0.0004483 0.0002845 1.516e-04"
#  [3] "Cumulative Proportion 0.9984176 0.9989678 0.9994161 0.9997006 9.999e-01"
#  [4] "                          PC25      PC26      PC27"                     
#  [5] "Proportion Explained  9.93e-05 3.036e-05 1.814e-05"                     
#  [6] "Cumulative Proportion 1.00e+00 1.000e+00 1.000e+00"                     
#  [7] "Scaling 2 for species and site scores"                                  
#  [8] "* Species are scaled proportional to eigenvalues"                       
#  [9] "* Sites are unscaled: weighted dispersion equal on all dimensions"      
# [10] "* General scaling constant of scores:  1.93676 "                        
# [11] ""                                                                       
# [12] ""                                                                       
# [13] "Species scores"                                                         
# [14] "BCO -0.20174  0.08807 -0.067086 -0.0529106  0.0737228  0.037312"        
# [15] "PCH -0.14717  0.05829 -0.067311 -0.0458414  0.0501013  0.031605"        
# [16] "GAR -0.35245 -0.14076  0.168014  0.0185946  0.0213462 -0.129788"        
# [17] "BBO -0.24317  0.03679 -0.082731 -0.0384489  0.0939828  0.063369"        
# [18] "ABL -0.42536 -0.26155 -0.054190  0.1021959 -0.0078085  0.044540"        
# [19] "ANG -0.20631  0.11889 -0.062079 -0.0175733  0.0718743 -0.001956"        
# [20] ""                                                                       
# [21] ""                                                                       
# [22] "Site scores (weighted sums of species scores)"
```
]

---

# Analyse en composantes principales

.pull-left[
`Species` fait référence à vos descripteurs (c.-à-d. les colonnes de votre ensemble de données), qui sont ici les espèces de poissons.

`Scores` font référence à la position de chaque espèce le long des composantes principales.
]

```
# [1] "Cumulative Proportion 1.00e+00 1.000e+00 1.000e+00"               
# [2] "Scaling 2 for species and site scores"                            
# [3] "* Species are scaled proportional to eigenvalues"                 
# [4] "* Sites are unscaled: weighted dispersion equal on all dimensions"
# [5] "* General scaling constant of scores:  1.93676 "                  
# [6] ""                                                                 
# [7] ""                                                                 
# [8] "Species scores"
```
]

```
# [1] "PCH -0.14717  0.05829 -0.067311 -0.0458414  0.0501013  0.031605"
# [2] "GAR -0.35245 -0.14076  0.168014  0.0185946  0.0213462 -0.129788"
# [3] "BBO -0.24317  0.03679 -0.082731 -0.0384489  0.0939828  0.063369"
# [4] "ABL -0.42536 -0.26155 -0.054190  0.1021959 -0.0078085  0.044540"
# [5] "ANG -0.20631  0.11889 -0.062079 -0.0175733  0.0718743 -0.001956"
# [6] ""                                                               
# [7] ""                                                               
# [8] "Site scores (weighted sums of species scores)"
```
]

.pull-right2[
`Sites` représentent les lignes de votre jeu de données, qui sont ici les différents sites le long de la rivière *Doubs*.
]

```r
scores(spe.h.pca,
       display = "species" or "sites")
```

]

---
### Analyse en composantes principales : condensation des données

Nous avons 27 composantes principales. Cependant, nous pouvons appliquer des algorithmes pour sélectionner le plus petit nombre de composantes principales qui rendent encore compte d'une grande variance dans les données.

#### Critère de Kaiser-Guttman

.pull-left[
Nous pouvons sélectionner les composantes principales qui capturent plus de variance que l'explication moyenne de toutes les composantes principales. Nous le faisons en

1. Extraire les *valeurs propres* associées aux composantes principales ;

2. Sélectionner les *valeurs propres* au-dessus de la *valeur propre* moyenne :

```r
ev <- spe.h.pca$CA$eig
# ev[ev > mean(ev)]
```
]

```r
n <- length(ev)
barplot(ev, main = "Eigenvalues", col = "grey", las = 2)
abline(h = mean(ev), col = "red3", lwd = 2)
legend("topright", "Average eigenvalue",
 lwd = 2, col = "red3" , bty = "n")
```

<img src="workshop09-pres-fr_files/figure-html/unnamed-chunk-78-1.png" width="720" style="display: block; margin: auto;" />
]

---
### Analyse en composantes principales : condensation des données

#### Modèle *broken-stick*

.pull-left[
Le modèle *broken-stick* retient les composantes qui expliquent plus de variance que ce qui serait attendu en divisant aléatoirement la variance en `$p$` parties.

```r
head(bstick(spe.h.pca))
#        PC1        PC2        PC3        PC4        PC5        PC6 
# 0.07242581 0.05381432 0.04450858 0.03830475 0.03365187 0.02992957
```
]

```r
screeplot(spe.h.pca,
          bstick = TRUE, type = "lines")
```

<img src="workshop09-pres-fr_files/figure-html/unnamed-chunk-80-1.png" width="324" style="display: block; margin: auto;" />
]

---
## Analyse en composantes principales

Il ne reste plus qu'à discuter de la *mise à l'échelle* et à *visualiser* nos résultats.

Pratiquons et calculons une ACP sur les variables environnementales standardisées pour le même ensemble de données.

```r
env.pca <- rda(env.z)
# summary(env.pca, scaling = 2)
```

Déterminons notre ensemble de *valeurs propres* et leurs *vecteurs propres* correspondants :

```r
ev <- env.pca$CA$eig
```

```r
ev[ev>mean(ev)]
#      PC1      PC2      PC3 
# 6.097995 2.167126 1.037603
```
]

.pull-right[
<img src="workshop09-pres-fr_files/figure-html/unnamed-chunk-84-1.png" width="576" style="display: block; margin: auto;" />
]

---

# Analyse en composantes principales : `plot()`

L'information calculée par l'ACP peut être représentée par des *biplots*.

Nous pouvons produire un biplot rapide de l'ACP en utilisant la fonction `plot()` en base `R`.

```r
plot(spe.h.pca)
```

---

# Analyse en composantes principales : `plot()`

`biplot()` de `base` `R` permet une meilleure interprétation.

```r
biplot(spe.h.pca)
```

<img src="workshop09-pres-fr_files/figure-html/unnamed-chunk-86-1.png" width="468" style="display: block; margin: auto;" />
]

Les flèches sont tracées pour montrer la directionnalité et l'angle des descripteurs dans l'ordination.

Les descripteurs situés à 180 degrés les uns des autres sont corrélés négativement ;

Les descripteurs situés à 90 degrés l'un de l'autre ont une corrélation nulle ;

Les descripteurs situés à 0 degré l'un de l'autre sont positivement corrélés.

]

---
# Analyse en composantes principales : *mise à l'échelle*

.small[
*Mise à l'échelle de type 2* (`scaling = 2`, défaut) : les distances entre les objets ne sont pas des approximations des distances euclidiennes ; les angles entre les vecteurs des descripteurs (espèces) reflètent leurs corrélations.
]

.small[
*Mise à l'échelle de type 1* (`scaling = 1`) : tente de préserver la distance euclidienne (dans un espace multidimensionnel) entre les objets (sites) : les angles entre les vecteurs des descripteurs (espèces) ne sont pas pertinents.
]

```r
biplot(spe.h.pca, scaling = 1)
```

<img src="workshop09-pres-fr_files/figure-html/unnamed-chunk-87-1.png" width="324" style="display: block; margin: auto;" />
]

```r
biplot(spe.h.pca, scaling = 2)
```

<img src="workshop09-pres-fr_files/figure-html/unnamed-chunk-88-1.png" width="324" style="display: block; margin: auto;" />
]

???
2 : **Meilleur pour l'interprétation des relations entre les descripteurs (espèces)**.

1 : **Meilleur pour l'interprétation des relations entre les objets (sites)!**

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

En utilisant tout ce que vous avez appris, calculez une ACP sur les données d'abondance des espèces d'acariens.

```r
data(mite)
```

Soyez prêt à discuter et à répondre :
- Quelles sont les composantes principales *les plus pertinentes*, c.-à-d. les sous-ensembles ?
- Quels groupes de sites pouvez-vous identifier ?
- Quels groupes d'espèces sont liés à ces groupes de sites ?

---

# Défi #2: Solution

Calculer l'ACP sur les données d'espèces transformées par Hellinger.

```r
mite.spe.hel <- decostand(mite,
 method = "hellinger")

mite.spe.h.pca <- rda(mite.spe.hel)
```

```r
ev <- mite.spe.h.pca$CA$eig
ev[ev>mean(ev)]
n <- length(ev)
barplot(ev, main = "Eigenvalues",
 col = "grey", las = 2)
abline(h = mean(ev),
 col = "red3", lwd = 2)
legend("topright",
 "Average eigenvalue",
 lwd = 2,
 col = "red3", bty = "n")
```
]

.pull-right[
<img src="workshop09-pres-fr_files/figure-html/unnamed-chunk-92-1.png" width="360" style="display: block; margin: auto;" />
]

---
# Défi #2: Solution

```r
biplot(mite.spe.h.pca,
       col = c("red3", "grey15"))
```

---
# Analyse des coordonnées principales

La **ACoP** (*PCoA*, en anglais) est similaire dans son esprit à l'ACP, mais elle prend des *dissimilarités* comme données d'entrée !

Elle vise à représenter fidèlement les distances avec l'espace dimensionnel le plus bas possible.

Elle commence par le (i) calcul d'une matrice de distance pour les `$p$` éléments, puis (ii) le centrage de la matrice par lignes et colonnes, et enfin, la (iii) *decomposition* de la matrice de distance centrée en éléments propres.

Pour calculer une ACoP, nous pouvons utiliser les fonctions `cmdscale()` ou `pcoa()` des paquets `stats` et `ape` :

```r
library(ape)
spe.h.pcoa <- pcoa(dist(spe.hel))
summary(spe.h.pcoa)
# Length Class Mode 
# correction 2 -none- character
# note 1 -none- character
# values 5 data.frame list 
# vectors 783 -none- numeric 
# trace 1 -none- numeric
```

---

# Analyse des coordonnées principales

```r
head(spe.h.pcoa$values)
#   Eigenvalues Relative_eig Broken_stick Cumul_eig Cumul_br_stick
# 1   7.2228939   0.51334374   0.14412803 0.5133437      0.1441280
# 2   1.7987449   0.12783995   0.10709099 0.6411837      0.2512190
# 3   1.2970423   0.09218307   0.08857247 0.7333668      0.3397915
# 4   1.0780684   0.07662021   0.07622679 0.8099870      0.4160183
# 5   0.6150273   0.04371107   0.06696753 0.8536980      0.4829858
# 6   0.4691296   0.03334186   0.05956013 0.8870399      0.5425459
```

---

# Analyse des coordonnées principales

Nous pouvons également voir les *vecteurs propres* associés à chaque *valeur propre* contenant les coordonnées dans l'espace euclidien pour chaque site.

```r
head(spe.h.pcoa$vectors)[, 1:5]
#         Axis.1      Axis.2      Axis.3     Axis.4       Axis.5
# 1 -0.509824403 -0.27654372  0.64011383 -0.3393734  0.207330880
# 2 -0.698794880 -0.03935586  0.11324989 -0.2328859 -0.157730682
# 3 -0.640690642  0.01566707  0.03835044 -0.2669706 -0.125293094
# 4 -0.413985947  0.10477084 -0.15728486 -0.2851828 -0.001250382
# 5  0.003083242  0.05284310 -0.32206098 -0.2730693  0.315944703
# 6 -0.295314224  0.05778805 -0.32395301 -0.2262902  0.056493824
```

---

# Analyse des coordonnées principales : `biplot.pcoa()`

Nous pouvons afficher les distances entre les sites en utilisant la fonction `biplot.pcoa()`, ainsi que représenter les espèces associées à chaque site.

```r
biplot.pcoa(spe.h.pcoa, spe.hel)
```

---

##### Analyse des coordonnées principales : distances non métriques

La ACoP peut également être utilisée pour capturer les informations contenues dans les distances non métriques, telles que la populaire distance de Bray-Curtis. Essayons :

```r
spe.bray.pcoa <- pcoa(spe.D.Bray)
```

```r
spe.bray.pcoa$values$Eigenvalues
#  [1]  3.695331e+00  1.098472e+00  7.104740e-01  4.149729e-01  3.045604e-01
#  [6]  1.917884e-01  1.569703e-01  1.319099e-01  1.294251e-01  8.667896e-02
# [11]  4.615780e-02  3.864487e-02  2.745800e-02  1.306508e-02  7.087896e-03
# [16]  4.039469e-03  1.300594e-03  0.000000e+00 -3.534426e-05 -3.940676e-03
# [21] -8.956051e-03 -1.461149e-02 -1.598905e-02 -2.145686e-02 -3.017013e-02
# [26] -3.432671e-02 -3.760052e-02 -6.037087e-02 -6.880061e-02
```
]

Cela est dû au fait que les distances non métriques ne peuvent pas être représentées dans l'espace euclidien sans corrections (*voir* Legendre & Legendre 2012 pour plus de détails à ce sujet) :

```r
spe.bray.pcoa <- pcoa(spe.D.Bray,
 correction = "cailliez")
```
]

---
##### Analyse des coordonnées principales : distances non métriques

Les valeurs propres corrigées sont maintenant sur une nouvelle colonne !

```r
spe.bray.pcoa$values$Corr_eig
#  [1] 5.20461437 1.60465006 1.09152082 0.68985417 0.52129425 0.38710929
#  [7] 0.36447367 0.29671255 0.27559544 0.21600584 0.15419396 0.15378333
# [13] 0.11812808 0.08848541 0.07304055 0.06999353 0.05712927 0.05587583
# [19] 0.05432215 0.04912221 0.04100207 0.03777775 0.03451234 0.02959507
# [25] 0.02436729 0.01902747 0.00284371 0.00000000 0.00000000
```

]

```r
biplot.pcoa(spe.bray.pcoa)
```

]

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

Calculez une ACoP sur les données d'abondance des espèces d'acariens transformées par Hellinger.

Soyez prêt à répondre :

- Quels sont les *vecteurs propres* et les *valeurs propres* significatifs ?
- Quels groupes de sites pouvez-vous identifier ?
- Quels groupes d'espèces sont liés à ces groupes de sites ?
- Comment les résultats de la ACoP se comparent-ils à ceux de l'ACP ?

---
# Défi #3: Solution

- Transformation de Hellinger des données sur les espèces

```r
mite.spe <- mite
mite.spe.hel <- decostand(mite.spe, method = "hellinger")
```

- Calcul de la ACoP

```r
mite.spe.h.pcoa <- pcoa(dist(mite.spe.hel))
```

---
# Défi #3: Solution

- Construisez un biplot pour visualiser les données :

```r
biplot.pcoa(mite.spe.h.pcoa, mite.spe.hel)
```

---
# Positionnement multidimensionnel non-métrique

- Dans l'ACP et la ACoP, les objets sont ordonnés selon un petit nombre de dimensions (généralement > 2) ;

- Les biplots 2D peuvent ne pas représenter toute la variation au sein de l'ensemble de données ;

- Parfois, nous cherchons à représenter les données dans un plus petit nombre de dimensions spécifié ;

- Comment pouvons-nous tracer l'espace d'ordination pour représenter la plus grande variation possible dans les données ?

Nous pouvons essayer d'utiliser le *positionnement multidimensionnel non-métrique*!

* Le PMnM (en anglais, *nMDS* de *non-metric multidimensional scaling*) est la contrepartie non-métrique de ACoP ;
* Il utilise un algorithme d'optimisation itératif pour trouver la meilleure représentation des distances dans un espace réduit ;

---
# Positionnement multidimensionnel non-métrique

- Le PMnM (ou nMDS) applique une procédure itérative qui tente de positionner les objets dans le nombre de dimensions demandé de manière à minimiser une fonction de contrainte (échelonnée de `$0$` à `$1$`), qui mesure la qualité de l'ajustement de la distance dans la configuration de l'espace réduit.

- Par conséquent, plus la valeur de stress est faible, meilleure est la représentation des objets dans l'espace d'ordination.

- nMDS est implémenté dans `vegan` comme `metaMDS()` où :
  - `distance` spécifie la métrique de distance à utiliser ;
  - `k` spécifie le nombre de dimensions.

```r
spe.nmds <- metaMDS(spe, distance = 'bray', k = 2)
```

---
# Positionnement multidimensionnel non-métrique: *qualité de l'ajustement*

Le diagramme de *Shepard* et les valeurs de contraintes peuvent être obtenus avec `stressplot()` :

```r
spe.nmds$stress # [1] 0.07376229
stressplot(spe.nmds, main = "Shepard plot")
```

.pull-left[
<img src="workshop09-pres-fr_files/figure-html/unnamed-chunk-109-1.png" width="360" style="display: block; margin: auto;" />
]
.pull-right[
 
Le graphique de Shepard identifie une forte corrélation entre la dissimilarité observée et la distance d'ordination `$(R^2 > 0.95)$`, ce qui met en évidence la qualité élevée de l'ajustement du PMnM.

]

---
# Positionnement multidimensionnel non-métrique : `biplot()`

Construire un biplot

```r
plot(spe.nmds, type = "none",
     main = paste("NMDS/Bray - Stress =",
                  round(spe.nmds$stress, 3)),
     xlab = c("NMDS1"), ylab = "NMDS2")

points(scores(spe.nmds, display = "sites",
              choiches = c(1,2),
              pch = 21,
              col = "black",
              g = "steelblue",
              cex = 1.2))
text(scores(spe.nmds, display = "species", choices = c(1)),
            scores(spe.nmds, display = "species", choices = c(2)),
            labels = rownames(scores(spe.nmds, display = "species")),
            col = "red", cex = 0.8)
```

---
##### Positionnement multidimensionnel non-métrique : `biplot()`

Le biplot du PMnM montre un groupe de sites fermés caractérisés par les espèces BLA, TRU, VAI, LOC, CHA et OMB, tandis que les autres espèces forment un groupe de sites dans la partie supérieure droite du graphique.

Quatre sites dans la partie inférieure du graphique sont fortement différents des autres.
]

.pull-right[
<img src="workshop09-pres-fr_files/figure-html/unnamed-chunk-111-1.png" width="360" style="display: block; margin: auto;" />
]

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

Exécutez le PMnM de l'abondance des espèces d'acariens `mite.spe` en 2 dimensions basé sur une distance de Bray-Curtis.

Évaluez la qualité de l'ajustement de l'ordination et interprétez le biplot.

Rappelez-vous de ces fonctions utiles:

```r
?metaMDS
?stressplot
```

---
# Défi #4: Solution

Exécutez le PMnM de l'abondance des espèces d'acariens `mite.spe` en 2 dimensions basé sur une distance de Bray-Curtis.

```r
mite.nmds <- metaMDS(mite.spe, distance = 'bray', k = 2)
```

---
# Défi #4: Solution

```r
plot(mite.nmds, type = "none",
     main = paste("NMDS/Bray - Stress =",
                  round(mite.nmds$stress, 3)),
     xlab = c("NMDS1"), ylab = "NMDS2")
points(scores(mite.nmds, display = "sites",
              choices = c(1,2)),
              pch = 21,
              col = "black",
              bg = "steelblue",
              cex = 1.2)
text(scores(mite.nmds, display = "species", choices = c(1)),
     scores(mite.nmds, display = "species", choices = c(2)),
     labels = rownames(scores(mite.nmds, display = "species")),
     col = "red", cex = 0.8)
```

---
# Défi #4: Solution

Il n'y a pas vraiment de regroupement évident de sites dans le biplot NMDS. Cela nous indique que les espèces sont présentes dans la plupart des sites. Seuls quelques sites contiennent des communautés spécifiques.

---
# Défi #4: Solution

```r
stressplot(mite.nmds)
```

La dissimilarité observée est fortement corrélée avec la distance d'ordination, et la valeur de stress NMDS est relativement faible, ce qui nous indique que l'ordination NMDS est relativement précise.

---
# Conclusion

.alert[De nombreuses techniques d'ordination existent, mais leur spécificité doit guider vos choix sur les méthodes à utiliser.]

| | Distance préservée | Variables | Nombre maximal d'axes |
|---|---------|--------------|------|
|PCA| Euclidienne | Données quantitatives, relations linéaires | p|
|CA| Chi2 | Non négatif, données quantitatives homogènes, données binaires | p-1 |
|PCoA| Définie par l'utilisateur | Données quantitatives, semi-quantitatives, mixtes| p-1|
|NMDS| Définie par l'utilisateur | Données quantitatives, semi-quantitatives et mixtes | Définie par l'utilisateur |

---
# C'est l'heure du quiz !

Analyse en composantes principales

.alert[Laquelle de ces méthodes est la meilleure pour visualiser les *distances* entre composition des communautés de différents sites?]

Analyse en coordonnées principales (ACoP)

La proportion de variance capturée par une composante principale

---
# C'est l'heure du quiz!

Trouvez l'erreur dans cette figure!

![:scale 90%](images/Chall7.png)

--
.alert[
- Les données ne sont pas centrées. Beurk!
]

---
class: inverse, center, bottom

# Merci pour votre participation à cet atelier!

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