Atelier 10: Analyses multivariées avancées

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

---

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

Mathieu Vaillancourt

Katherine Hébert

Pedro Henrique P. Braga

Gabriel Muñoz

Kevin Cazelles

Marie Hélène Brice

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

Marie Hélène-Brice

Pedro Henrique P. Braga

Katherine Hébert

]
]

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

Monica Granados

Emmanuelle Chrétien

Bérenger Bourgeois

Amanda Winegardner

Xavier Giroux-Bougard

Vincent Fugère

]
]

.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/).

Vous devez également avoir installé ces librairies R:
* [Hmisc](https://cran.r-project.org/package=Hmisc)
* [labdsv](https://cran.r-project.org/package=labdsv)
* [MASS](https://cran.r-project.org/package=MASS)
* [vegan](https://cran.r-project.org/package=vegan)
* [ggplot2](https://cran.r-project.org/package=ggplot2)

```R
install.packages(c('Hmisc', 'labdsv', 'MASS', 'vegan', 'ggplot2'))
```

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

???

Paquets non-utilisés dans la version de 3 heures
[`mvpart`](https://cran.r-project.org/package=mvpart)
remotes::install_github("cran/mvpart")

---
# Objectifs d'apprentissage

__Objectif__: Réaliser des analyses multivariées avancées sur des données de communauté.

En utilisant des méthodes d'ordination sous contrainte, nous apprendrons à:

- Explorer comment les variables environnementales __expliquent__ des patrons de composition en espéces à travers différents sites.

- __Décrire et prédire__ ces tendances entre la composition des communautés et des variables environnementales.

---
class: inverse, center, middle

# Introduction

---

# Introduction

Cet atelier est une continuation de l'atelier précédent (#9), qui a offert un aperçu des analyses multivariées non-contraintes:

* Mesures de distances et transformations de données
* Groupement hiérarchique
* Ordinations non-contraintes (PCA, PCoA, CA, nmDS)

Celles-ci permettent de relever des ***patrons*** dans la structure des communautés d'espèces ou des descripteurs.

L'atelier #10 montre des analyses permettant d'explorer comment les variables environnementales ***expliquent*** ces patrons.

---
# Introduction

Cet atelier se concentrera sur les analyses .alert[**sous contraintes**]:

* Analyse canonique de redondances (ACR ou RDA)
* Analyse canonique de redondances partielle (ACR partielle ou RDA partielle)
* Partitionnement de la variation
* Analyse linéaire discriminante (ALD ou LDA)

Ces analyses nous permettrons de **décrire** et de **prédire** les relations entre la structure des communautés et les variables environnementales.

On pourra alors ***tester des hypothèses***!

???

Pas présenté dans la version de 3 heures:
Arbre de régression multivarié (ARM ou MRT)

---
# Code et données

Lien vers le matériel de l'atelier:[ github.com/QCBSRworkshops/workshop10](https://github.com/QCBSRworkshops/workshop10).

Téléchargez le code R et les données requises pour cet atelier:
* [Code R](https://qcbsrworkshops.github.io/workshop10/book-fr/workshop10-script-fr.R)
* Données:
  * [DoubsEnv](https://raw.githubusercontent.com/QCBSRworkshops/workshop10/main/pres-fr/data/doubsenv.csv)
  * [DoubsSpe](https://raw.githubusercontent.com/QCBSRworkshops/workshop10/main/pres-fr/data/doubsspe.csv)
  * [DoubsSpa](https://raw.githubusercontent.com/QCBSRworkshops/workshop10/main/pres-fr/data/doubsspa.csv)
  * [Données test pour l'analyse linéaire discriminante](https://raw.githubusercontent.com/QCBSRworkshops/workshop10/main/pres-fr/data/classifyme.csv)

.small[.comment[Pour télécharger les données, faire `clic droit + Sauvegarder` sur la page qui ouvre.]]

---
# Librairies R

Assurez-vous de télécharger, d'installer et de charger les librairies suivants:
* [Hmisc](https://cran.r-project.org/package=Hmisc)
* [labdsv](https://cran.r-project.org/package=labdsv)
* [MASS](https://cran.r-project.org/package=MASS)
* [vegan](https://cran.r-project.org/package=vegan)

???

Unused packages in 3-hr version:
.comment[mvpart] package .alert[(voir code R)]

---
# Suivre l'atelier

Quelques conseils:
* Créez votre propre code (ou commentez le code R fourni);
* Évitez de copier-coller, ou d'exécuter le code directement du script fourni;
* N'oubliez pas de bien définir votre répertoire de travail vers le dossier contenant les fichiers requises pour l'atelier.

---
class: inverse, center, middle

# Exploration et préparation des données

---
# Introduction aux données

#### Rivière Doubs (Verneaux 1973)

.pull-left[Données d'abondances d'espèces de communautés de poissons de la rivière Doubs.
 * 27 espèces
 * 30 sites
 * 11 variables environnementales

![:scale 100%](images/Doubs_Laissey.jpg)
 ]

![:scale 100%](images/Doubs_carte.jpg)]

---
# Charger les données

Chargez la matrice d'abondances d'espèces (`doubsspe.csv`).

```r
# Assurez vous que les fichiers se trouvent dans votre répertoire de travail!
spe <- read.csv("data/doubsspe.csv", row.names = 1)
spe <- spe[-8,] # Supprimer site 8 (pas d'espèces).
```

Chargez la matrice de données environnementales (`doubsenv.csv`).

```r
env <- read.csv("data/doubsenv.csv", row.names = 1)
env <- env[-8,] # Supprimer site 8 (pas d'espèces).
```

---
# Données d'abondances d'espèces

Explorons la matrice des abondances de poissons.

```r
names(spe) # noms d'objets (espèces)
#  [1] "CHA" "TRU" "VAI" "LOC" "OMB" "BLA" "HOT" "TOX" "VAN" "CHE" "BAR" "SPI"
# [13] "GOU" "BRO" "PER" "BOU" "PSO" "ROT" "CAR" "TAN" "BCO" "PCH" "GRE" "GAR"
# [25] "BBO" "ABL" "ANG"
dim(spe) # dimensions de la matrice
# [1] 29 27
```

Et, si on veut plus de détails sur les espèces:

```r
head(spe) # 5 premières lignes
str(spe) # structure d'objets de la matrice
summary(spe) # statistiques descriptives des objets (min, moyenne, max, etc.)
```

---
# Distribution des abondances d'espèces

Explorons la structure de la communauté.

```r
# Compter la fréquence d'espèces dans chaque classe d'abondance
ab <- table(unlist(spe))
# Visualiser cette distribution
barplot(ab, las = 1,
 xlab = "Abundance class", ylab = "Frequency",
 col = grey(5:0/5))
```

<img src="workshop10-pres-fr_files/figure-html/unnamed-chunk-5-1.png" width="432" style="display: block; margin: auto;" />
.alert[Note]: Il y a beaucoup de zéros.

---
# Distribution des abondances d'espèces

Combien y a-t-il de zéros?

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

Quelle proportion de l'ensemble des données cela représente-t-il ?

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

---
# Transformation des données d'abondances

.alert[**Plus de 50%**] des données d'abondances sont en fait des absences. C'est élevé, mais pas inhabituel pour ce type de données.

Par contre, on ne veut pas que les **double zéros** amplifient artificiellement la similarité entre sites.
* Pour éviter ceci, on peut .alert[transformer] les données d'abondances.

Nous appliquerons alors une transformation .alert[Hellinger]:

```r
# la fonction decostand() dans la libraire vegan nous facilite la tâche:
library(vegan)
spe.hel <- decostand(spe, method = "hellinger")
```

---
# Données environnementales

Explorons les données environnementales.

```r
names(env) # noms des objets (variables environnementales)
#  [1] "das" "alt" "pen" "deb" "pH"  "dur" "pho" "nit" "amm" "oxy" "dbo"
dim(env) # dimensions de la matrice
# [1] 29 11
head(env) # 5 premières lignes
#    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
```

Et, si on veut plus de détails sur les variables environnementales:

```r
str(env) # structure des objets
summary(env) # statistiques descriptives (min, moyenne, max, etc.)
```

---
# Colinéarité

```r
# On peut également détecter (visuellement) les colinéarités entres variables:
heatmap(abs(cor(env)), # corrélation de Pearson (note: ce sont des valeurs absolues!)
        col = rev(heat.colors(6)),
        Colv = NA, Rowv = NA)
legend("topright",
       title = "R de Pearson",
       legend =  round(seq(0,1, length.out = 6),1),
       y.intersp = 0.7, bty = "n",
       fill = rev(heat.colors(6)))
```

<img src="workshop10-pres-fr_files/figure-html/unnamed-chunk-11-1.png" width="576" style="display: block; margin: auto;" />
]
.alert[Note:] Colinéarité entre quelques variables... .small[(das vs. alt, das vs. deb, das vs. dur, das vs. nit, oxy vs. dbo, etc.)]

---
# Standardisation des données

Il est impossible de comparer les effets de variables  qui ont des unités différentes.

Avant d'effectuer les analyses qui suivent, les données doivent donc être .alert[standardisées].

```r
# standardiser les données
env.z <- decostand(env, method = "standardize")

# centrer les données (moyenne ~ 0)
round(apply(env.z, 2, mean), 1)
# das alt pen deb  pH dur pho nit amm oxy dbo 
#   0   0   0   0   0   0   0   0   0   0   0

# réduire les données (écart type = 1)
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
```

---

# Analyses canoniques

---
# Analyses canoniques

Les analyses canoniques nous permettent:
- d'identifier les **relations** entre un ensemble de variables réponses et un ensemble de variables explicatives
- de tester des .alert[hypothèses écologiques] à propos de ces relations
- de faire des prédictions

---

# Analyses canoniques
## Analyse canonique de redondances (RDA)

---
# Analyse canonique de redondances (RDA)

L'analyse canonique de redondances est une ordination sous contraintes.
- La RDA est une extension directe de la régression multiple.
- Elle modélise l’effet d’une matrice explicative sur une matrice réponse .comment[(au lieu d'une seule variable réponse)].

Les variables peuvent être quantitatives, qualitatives, ou binaires (0/1).
- .alert[transformez] et .alert[standardisez] les variables avant d'effectuer une RDA.

---
# Effectuer une RDA dans R

**Étape 1**: Transformer et/ou standardiser les données.

```r
# On utilisera nos données explicatives standardisées
# Enlever la variable "distance from the source" (colinéarité avec autres variables)
env.z <- subset(env.z, select = -das)
```

--
**Étape 2**: Effectuer une RDA.

```r
# Modèlise l'effect de tous les variables environnementales sur la composition en espèces des communautés
spe.rda <- rda(spe.hel ~ ., data = env.z)
```

**Étape 3**: Extraire les résultats de la RDA.

```r
summary(spe.rda)
...
# Partitioning of variance:
#               Inertia Proportion
# Total          0.5025     1.0000
# Constrained    0.3689     0.7341
# Unconstrained  0.1336     0.2659
...
```

---
# RDA output in R

Pour interpréter les résultats d'une RDA, on peut d'abord se concentrer sur cette partie clé de la sortie:

```
...
# Partitioning of variance:
#               Inertia Proportion
# Total          0.5025     1.0000
# Constrained    0.3689     0.7341
# Unconstrained  0.1336     0.2659
...
```

* **Constrained Proportion**: variance de `$Y$` expliquée par `$X$` .alert[(73.41%)]
* **Unconstained Proportion**: variance in `$Y$` non expliquée par .alert[(26.59%)]

Comment présenteriez-vous ces résultats?

* .comment[Les variables environnementales mesurées expliquent .alert[73.41%] de la variation dans la composition en espèces des communautés de poissons dans la rivière Doubs.]

---
# Sélection de variables

Une .alert[sélection progressive] peut être effectuée afin de sélectionner les variables explicatives significatives.

**Quelles variables contribuent de façon significative au pouvoir explicatif du modèle?**
--

```r
# Sélection progressive de variables:
fwd.sel <- ordiR2step(rda(spe.hel ~ 1, data = env.z), # modèle le plus simple
 scope = formula(spe.rda), # modèle "complet"
 direction = "forward",
 R2scope = TRUE, # limité par le R2 du modèle "complet"
 pstep = 1000,
 trace = FALSE) # mettre TRUE pour voir le processus du sélection!
```

.small[Essentiellement, on ajoute une variable à la fois au modèle, et on retient la variable si elle augmente significativement le `$R^2$` ajusté du modèle.]

???

Notez que la sélection progressive est un outil d'aide à la décision, mais ne remplace pas votre jugement éclairé! La sélection finale des variables doit s'appuyer sur un raisonnement écologique. Ce n'est pas parce qu'une variable d'intérêt écologique n'est pas sélectionnée par sélection progressive qu'il faut nécessairement l'exclure.

---
# Sélection de variables

- Quelles variables ont été sélectionnées?

```r
fwd.sel$call
# rda(formula = spe.hel ~ alt + oxy + dbo, data = env.z)
```

- Quel est le R2 ajusté d'une RDA incluant seulement les variables significatives?

```r
# Écrire notre nouveau modèle
spe.rda.signif <- rda(spe.hel ~ alt + oxy + dbo, data = env.z)
# vérifier son R2 ajusté
RsquareAdj(spe.rda.signif)
# $r.squared
# [1] 0.5894243
# 
# $adj.r.squared
# [1] 0.5401552
```

---
# Tester la significativité d'une RDA

Utilisez `anova.cca()` pour tester la significativité globale de notre RDA.

```r
anova.cca(spe.rda.signif, permutations = 1000)
...
#          Df Variance      F   Pr(>F)    
# Model     3  0.29619 11.963 0.000999 ***
# Residual 25  0.20632                    
# ---
...
```

On peut aussi tester la significativité de chaque variable avec `by = 'term'`!

```r
anova.cca(spe.rda.signif, permutations = 1000, by = "term")
...
# Model: rda(formula = spe.hel ~ alt + oxy + dbo, data = env.z)
#          Df Variance       F   Pr(>F)    
# alt       1 0.164856 19.9759 0.000999 ***
# oxy       1 0.082426  9.9877 0.000999 ***
# dbo       1 0.048909  5.9264 0.001998 ** 
# Residual 25 0.206319                     
...
```

---
# Représentation graphique des RDAs

Une RDA permet la **visualisation simultanée** des variables réponses et explicatives .comment[(i.e. espèces et variables environnementales)].
 
--
 
Comme pour la PCA, on doit choisir entre deux options de cadrage:

| Type 1                                           | Type 2                                            |
| -------------------------------------------------|---------------------------------------------------|
| distances entre objets ≈ distances euclidiennes  |    angles entre variables ≈ leur corrélation      |

---
# Triplot RDA: Cadrage de type 1

```r
ordiplot(spe.rda.signif,
         scaling = 1,
         type = "text")
```

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

.pull-right[
######Le cadrage 1 permet d'interpréter les distances entre les objets dans la **matrice réponse**.

.small[
* Les communautés dans les sites (chiffres) .alert[plus rapprochés] ont des compositions *plus similaires*.
* Les espèces .alert[plus rapprochées] occupent souvent les *mêmes sites*.
]
]

---
# Triplot RDA: Cadrage de type 2

```r
ordiplot(spe.rda.signif,
         scaling = 2,
         type = "text")
```

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

.pull-right[
######Le cadrage 2 montre les effets des **variables explicatives**.
.small[
* .alert[Longues] flèches = cette variable explique fortement la variation dans la matrice d'abondances.
* Flèches pointant des .alert[directions opposées] = montrent une relation *négative*.
* Flèches pointant la .alert[même direction] = montrent une relation *positive*.
]
]

---
# Configuration des triplots RDA

Les fonctions `plot()` et `ordiplot()` produisent des triplots rapidement et facilement, mais on peut aussi configurer les graphiques avec l'extraction de scores avec `scores()` et leur visualisation avec `points()`, `text()`, et `arrows()`.

Voir le script et le [livre](https://qcbsrworkshops.github.io/workshop10/book-fr/index.html) de l'atelier 10 pour coder la figure!

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

Effectuez une RDA pour modéliser les effets des variables environnementales sur l’abondance des espèces d'acariens.

```r
# Charger les données d'abondance des espèces d'acariens
data("mite")

# Charger les données environnementales
data("mite.env")
```

Rappel de fonctions utiles:

```r
decostand()
rda()
ordiR2step()
anova.cca()
ordiplot()
```
]

---
# Défi 1: Solution

**Étape 1:** Transformer et standardiser les données.

```r
# Transformer les données d'abondances
mite.spe.hel <- decostand(mite, method = "hellinger")

# Standardiser les données environmentales quantiatives
mite.env$SubsDens <- decostand(mite.env$SubsDens, method = "standardize")
mite.env$WatrCont <- decostand(mite.env$WatrCont, method = "standardize")
```

---
# Défi 1: Solution

**Étape 2:** Sélectionner les variables environnementales.

```r
# RDA avec tous les variables environnementales
mite.spe.rda <- rda(mite.spe.hel ~ ., data = mite.env)

# Sélection progressive des variables environnementales significatives
fwd.sel <- ordiR2step(rda(mite.spe.hel ~ 1, data = mite.env),
 scope = formula(mite.spe.rda),
 direction = "forward",
 R2scope = TRUE, pstep = 1000, trace = FALSE)
fwd.sel$call
# rda(formula = mite.spe.hel ~ WatrCont + Shrub + Substrate + Topo, 
# data = mite.env)
```

---
# Défi 1: Solution

**Étape 3:** Effectuer la RDA et extraire le R2 ajusté.

```r
# Refaire la RDA avec seulement les variables significatives
mite.spe.rda.signif <- rda(mite.spe.hel ~ WatrCont + Shrub +
 Substrate + Topo + SubsDens,
 data = mite.env)

# Calculer le R2 ajusté
RsquareAdj(mite.spe.rda.signif)$adj.r.squared
# [1] 0.4367038
```

---
# Défi 1: Solution

**Étape 4:** Tester la significativité globale du modèle.

```r
anova.cca(mite.spe.rda.signif, step = 1000)
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = mite.spe.hel ~ WatrCont + Shrub + Substrate + Topo + SubsDens, data = mite.env)
#          Df Variance     F Pr(>F)    
# Model    11  0.20759 5.863  0.001 ***
# Residual 58  0.18669                 
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

Les variables environnementales sélectionnées expliquent .alert[43.7% (p = 0.001)] de la variation dans la composition de communautés des acariens.

---
# Défi 1: Solution

**Étape 5:** Visualiser le modèle à l'aide d'un triplot!
.pull-left[

```r
ordiplot(mite.spe.rda.signif,
         scaling = 1,
         frame = F,
         main = "Cadrage 1")
```

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

```r
ordiplot(mite.spe.rda.signif,
         scaling = 2,
         frame = F,
         main = "Cadrage 2")
```

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

---

# Analyses canoniques
## RDA partielle

---
# RDA partielle

La RDA partielle est un cas particulier de la RDA qui permet de tenir compte de **covariables**.

En d'autres mots, on peut modéliser les effets linéaires de la matrice `$X$` sur la matrice `$Y$`, tout en contrôlant pour l'effet de la matrice de covariables `$W$`.

---
# Applications de la RDA partielle

Puisque la RDA partielle contrôle pour l'effet de covariables, on peut:

* Évaluer l’effet de variables environnementales sur la composition des communautés en prenant en compte des covariables de **moindre intérêt**.
* **Isoler** les effets d'un ou plusieurs groupes de variables explicatives.

--
.center[.comment[Par exemple, la RDA partielle est souvent utilisée pour séparer les effets de variables environnementales des effets de variables spatiales.]]

---
# Exemple: RDA partielle sur les données Doubs

Dans `R`, on peut faire une RDA partielle avec la fonction `rda()`.

Par exemple, évaluons l'effet de la chimie de l'eau sur l’abondance des poissons (`spe.hel`) en tenant compte de covariables topographiques.

```r
# Divisez le tableau de données environnementales en deux:
# variables topographiques et chimiques
env.topo <- subset(env.z, select = c(alt, pen, deb))
env.chem <- subset(env.z, select = c(pH, dur, pho, nit, amm, oxy, dbo))

# Faire la RDA partielle
spe.partial.rda <- rda(spe.hel, env.chem, env.topo)
```

.small[
**Note:** On peut aussi utiliser une syntaxe de formule comme `Y ~ X + Condition(W)`, où `Condition()` permet de tenir compte de covariables.

```r
spe.partial.rda <- rda(spe.hel ~ pH + dur + pho + nit + amm + oxy + dbo +
 Condition(alt + pen + deb), # covariables ici
 data = env.z)
```
]

---
# Interprétation de la sortie d'une RDA partielle

```r
summary(spe.partial.rda)
```

```
...
# Partitioning of variance:
# Inertia Proportion
# Total 0.5025 1.0000
# Conditioned 0.2087 0.4153
# Constrained 0.1602 0.3189
# Unconstrained 0.1336 0.2659
...
```

* **Conditioned Proportion**: variance de `$Y$` expliquée par `$W$` .alert[(41.53%)]
* **Constrained Proportion**: variance de `$Y$` expliquée par `$X$` .alert[(31.89%)]
* **Unconstained Proportion**: variance de `$Y$` non expliquée .alert[(26.59%)]
 
 
.center[Comment présenteriez-vous ces résultats?]

.center[.comment[La chimie de l’eau explique .alert[31.89%] de l’abondance des espèces de poissons, tandis que la topographie explique .alert[41.53%] de la variation en abondances des poissons.]]

---
# Interprétation d'une RDA partielle

Il nous manque encore quelques détails importants pour notre interprétation du modèle!
 
 
.pull-left2[
1\. Quel est le **pouvoir explicatif** du modèle?

```r
# Extraire le R2 ajusté du modèle
RsquareAdj(spe.partial.rda)$adj.r.squared
# [1] 0.2413464
```
]

.pull-right2[
 
.small[.comment[Notre modèle explique .alert[24.1%] de la variation en abondance de poissons entre sites.]]
]

```r
# Évaluer la significativité statistique du modèle
anova.cca(spe.partial.rda, step = 1000)
...
# Permutation test for rda under reduced model
# Number of permutations: 999
# 
# Model: rda(X = spe.hel, Y = env.chem, Z = env.topo)
#          Df Variance      F Pr(>F)    
# Model     7  0.16024 3.0842  0.001 ***
# Residual 18  0.13360                  
...
```
]

---
# Représentation graphique

On peut visualiser les effets des variables environnementales sur la communauté de poissons avec la fonction `ordiplot()`.

```r
ordiplot(spe.partial.rda,
         scaling = 2,
         main = "Rivière Doubs - Cadrage 2")
```

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

???

Le cadrage 2 montre les effets des variables explicatives, donc de la matrice X sur la matrice Y. L'effet des covariables est 'pris en compte', et n'est pas réellement d'intérêt dans le modèle.

Rappel: En cadrage 2, la longueur des flèches montre l'amplitude de l'effet, leur direction montre la direction de la relation (direction opposée = négative, même direction = positive).

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

.small[
Effectuez une RDA partielle de l’abondance des espèces de mites (`mite.spe.hel`) en fonction des variables environnementales tenant compte de l’effet du substrat (`SubsDens`, `WaterCont` and `Substrate`).
* Quel pourcentage de variance est expliqué par les variables environnementales?
* Le modèle est-il significatif?
* Quels sont les axes significatifs?

Rappel des données et fonctions utiles:

```r
rda()
summary()
RsquareAdj()
anova.cca() # voir l'argument 'by' dans ?anova.cca
```
]

---
# Défi 2: Solution

**Étape 1:** Transformer et standardiser les données.

**Étape 2:** Faire la RDA partielle:

```r
mite.spe.subs <- rda(mite.spe.hel ~ Shrub + Topo
 + Condition(SubsDens + WatrCont + Substrate),
 data = mite.env)

# Extraire les résultats
summary(mite.spe.subs)
...
# Partitioning of variance:
#               Inertia Proportion
# Total         0.39428    1.00000
# Conditioned   0.16891    0.42839
# Constrained   0.03868    0.09811
# Unconstrained 0.18669    0.47350
...
```
--
Shrub et Topo expliquent .alert[9.8%] de la variation de l’abondance de mites, tandis que le substrat explique .alert[42.8%] de cette variation.

---
# Défi 2: Solution

**Étape 3:** Interpréter les résultats!

* Quel pourcentage de la variance est expliqué par les variables environnementales?

```r
RsquareAdj(mite.spe.subs)$adj.r.squared
# [1] 0.08327533
```
--
 
* Le modèle est-il significatif?

```r
anova.cca(mite.spe.subs, step = 1000)
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = mite.spe.hel ~ Shrub + Topo + Condition(SubsDens + WatrCont + Substrate), data = mite.env)
#          Df Variance     F Pr(>F)    
# Model     3 0.038683 4.006  0.001 ***
# Residual 58 0.186688                 
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---
# Défi 2: Solution

* Quels axes sont significatifs?

```r
anova.cca(mite.spe.subs, step = 1000, by = "axis")
# Permutation test for rda under reduced model
# Forward tests for axes
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = mite.spe.hel ~ Shrub + Topo + Condition(SubsDens + WatrCont + Substrate), data = mite.env)
#          Df Variance      F Pr(>F)    
# RDA1      1 0.027236 8.4618  0.001 ***
# RDA2      1 0.008254 2.5643  0.022 *  
# RDA3      1 0.003193 0.9919  0.389    
# Residual 58 0.186688                  
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---

# Analyses canoniques
## Partitionnement de la variation

---
# Partitionnement de la variation

Divise la variation d’une matrice de variables réponses en 2, 3, ou 4 matrices de variables explicatives.
* e.g. variables locales vs. à large échelle
* e.g. abiotique vs. biotique

---
# Partitionnement de la variation

---
# Partitionnement de la variation dans R

Pour démontrer le partitionnement de la variation dans `R`, nous allons partitionner la variation de la composition des espèces de poissons entre les variables chimiques et topographiques.

```r
spe.part.all <- varpart(spe.hel, env.chem, env.topo)
spe.part.all$part # extraire résultats
# No. of explanatory tables: 2 
# Total variation (SS): 14.07 
# Variance: 0.50251 
# No. of observations: 29 
# 
# Partition table:
# Df R.squared Adj.R.squared Testable
# [a+c] = X1 7 0.60579 0.47439 TRUE
# [b+c] = X2 3 0.41526 0.34509 TRUE
# [a+b+c] = X1+X2 10 0.73414 0.58644 TRUE
# Individual fractions 
# [a] = X1|X2 7 0.24135 TRUE
# [b] = X2|X1 3 0.11205 TRUE
# [c] 0 0.23304 FALSE
# [d] = Residuals 0.41356 FALSE
# ---
# Use function 'rda' to test significance of fractions of interest
```
]

---
# Diagramme de Venn

```r
plot(spe.part.all,
     Xnames = c("Chem", "Topo"), # noms des matrices explicatives
     bg = c("seagreen3", "mediumpurple"), alpha = 80,
     digits = 2,
     cex = 1.5)
```

---
# Tester la significativité

* La significativité de la fraction partagée [b] ne peut .alert[pas] être testée.
* Mais, on peut tester la significativité des autres fractions!

---
# Significativité: X1 [a+b]

```r
anova.cca(rda(spe.hel, env.chem))
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(X = spe.hel, Y = env.chem)
#          Df Variance      F Pr(>F)    
# Model     7  0.30442 4.6102  0.001 ***
# Residual 21  0.19809                  
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```
]

---
# Significativité: X2 [b+c]

```r
anova.cca(rda(spe.hel, env.topo))
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(X = spe.hel, Y = env.topo)
#          Df Variance     F Pr(>F)    
# Model     3  0.20867 5.918  0.001 ***
# Residual 25  0.29384                 
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```
]

---
# Significativité: Fractions individuelles

.center[![:scale 30%](images/VarPart_small.png)]
.small[
[a] Chimie (ajusté pour tenir compte de topographie)

```r
anova.cca(rda(spe.hel, env.chem, env.topo))
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(X = spe.hel, Y = env.chem, Z = env.topo)
#          Df Variance      F Pr(>F)    
# Model     7  0.16024 3.0842  0.001 ***
# Residual 18  0.13360                  
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---
# Significativité: Fractions individuelles

```r
anova.cca(rda(spe.hel, env.topo, env.chem))
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(X = spe.hel, Y = env.topo, Z = env.chem)
#          Df Variance      F Pr(>F)    
# Model     3 0.064495 2.8965  0.001 ***
# Residual 18 0.133599                  
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```
]

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

.small[
**Partitionnez la variation de l’abondance des espèces de mites entres des variables de substrat (`SubsDens`, `WatrCont`) et des variables spatiales significatives.**
* Quelle est la proportion de variance expliquée par le substrat? par l'espace?
* Quelles sont les fractions significatives?
* Diagramme de Venn des résultats!

Chargez les variables spatiales:

```r
data("mite.pcnm")
```

Rappel de fonctions utiles:

```r
ordiR2step()
varpart()
anova.cca(rda())
plot()
```
]

---
# Défi 3: Solution

**Étape 1:** Sélection de variables spatiales significatives

```r
# Modèle RDA avec tous les variables spatiales
full.spat <- rda(mite.spe.hel ~ ., data = mite.pcnm)

# Sélection progressive des variables spatiales
spat.sel <- ordiR2step(rda(mite.spe.hel ~ 1, data = mite.pcnm),
 scope = formula(full.spat),
 R2scope = RsquareAdj(full.spat)$adj.r.squared,
 direction = "forward",
 trace = FALSE)
spat.sel$call
# rda(formula = mite.spe.hel ~ V2 + V3 + V8 + V1 + V6 + V4 + V9 + 
# V16 + V7 + V20, data = mite.pcnm)
```

---
# Défi 3: Solution

**Étape 2:** Créer des sous-groupes de variables explicatives.

```r
# Variables de substrat
mite.subs <- subset(mite.env, select = c(SubsDens, WatrCont))

# Variables spatiales significatives
mite.spat <- subset(mite.pcnm,
 select = names(spat.sel$terminfo$ordered))
 # pour rapidement accèder aux variables sélectionnées
```

---
# Défi 3: Solution

**Étape 3:** Partitionnement de la variation de la matrice d'abondances.

```r
mite.part <- varpart(mite.spe.hel, mite.subs, mite.spat)
mite.part$part$indfract # extraire résultats
# Df R.squared Adj.R.squared Testable
# [a] = X1|X2 2 NA 0.05901929 TRUE
# [b] = X2|X1 10 NA 0.19415929 TRUE
# [c] 0 NA 0.24765221 FALSE
# [d] = Residuals NA NA 0.49916921 FALSE
```

* Quelle est la proportion de variance expliquée par le substrat?
 * .alert[5.9%]

* Quelle est la proportion de variance expliquée par l'espace?
  * .alert[19.4%]

---
# Défi 3: Solution

**Étape 4:** Quelles sont les fractions significatives?

**[a]: Substrat seulement**

```r
anova.cca(rda(mite.spe.hel, mite.subs, mite.spat))
...
# Model: rda(X = mite.spe.hel, Y = mite.subs, Z = mite.spat)
#          Df Variance      F Pr(>F)    
# Model     2 0.025602 4.4879  0.001 ***
# Residual 57 0.162583                  
...
```

**[c]: Espace seulement**

```r
anova.cca(rda(mite.spe.hel, mite.spat, mite.subs))
...
# Model: rda(X = mite.spe.hel, Y = mite.spat, Z = mite.subs)
#          Df Variance      F Pr(>F)    
# Model    10  0.10286 3.6061  0.001 ***
# Residual 57  0.16258                  
...
```

---
# Défi 3: Solution

**Étape 5:** Visualiser les résultats avec un diagramme de Venn.

```r
plot(mite.part,
     digits = 2,
     Xnames = c("Subs", "Space"), # titre des fractions
     cex = 1.5,
     bg = c("seagreen3", "mediumpurple"), # ajoutez des couleurs!
     alpha = 80)
```

---
# Défi 3: Solution

**Alors, quels sont les effets du substrat et de l'espace sur les abondances d'espèces de mites?**

En groupes, résumez et interprétez ces résultats comme si vous écriviez un article. **Indice:** Pourquoi trouve-t-on un effet si important de l'espace?

.comment[
L'espace explique la plupart de la variation dans la communauté: elle explique **19.4%** de la variation à elle seule, alors que le substrat à lui seul n'en explique que **~6%**. **24.8%** de la variation est expliqué *conjointement* par l'espace et le substrat.

L'effet de l'espace survient peut-être parce que:
* un processus spatial influence la communauté, comme la dispersion
* ou encore il manque une variable environnementale importante dans notre modèle, qui varie elle-même dans l'espace
]

???

Salles de répartition: 5-10 minutes.

Réponse: Cet effet important de l'espace pourrait être un signe qu'un processus écologique spatial est important ici (comme la dispersion, par exemple). Cependant, il pourrait également nous indiquer qu'il manque une variable environnementale importante dans notre modèle, qui varie elle-même dans l'espace !

Notez également que la moitié de la variation n'est pas expliquée par les variables que nous avons incluses dans le modèle (regardez les résidus !), le modèle pourrait donc théoriquement être amélioré.

---
exclude: true

---
exclude: true

# Arbre de régression multivarié (MRT)

L'arbre de régression multivarié (ARM ou MRT) est une méthode de groupement hiérarchique sous contrainte.
* Partitionne une matrice réponse quantitative (Y) en sous-groupes sous la contrainte d'une matrice de variables explicatives (X).

---
exclude: true

# Arbre de régression multivarié (MRT)

L'arbre de régression multivarié est formé de:
* .alert[Branches]: chacune des divisions d'un noeud
* .alert[Noeuds]: points où les données se divisent en 2 groupes (caractérisés par une valeur limite d'une variable explicative)
* .alert[Feuilles]: groupes terminaux de sites

---
exclude: true

# Arbre de régression multivarié (MRT)

Plusieurs avantages:
* Ne présume pas l'existence d'une relation linéaire entre les matrices Y et X
* Facile à interpréter et à visualiser
* Robuste en présence de valeurs manquantes ou de colinéarité(s) entre les descripteurs
* Les valeurs brutes peuvent être utilisées (sans transformation)

---
exclude: true

# MRT: La méthode

La méthode implique deux volets s'effectuant en parallèle:
1. .comment[Partitionnement] des données sous contrainte
2. .comment[Validation croisée] pour identifier l'arbre ayant le meilleur pouvoir prédictif.

Choisissez l'arbre selon les objectifs de votre étude. Généralement, on veut un arbre:
* .comment[parcimonieux]
* mais avec un nombre .comment[informatif] de groupes
* Essentiellement: quel arbre répond à votre question?

---
exclude: true

# MRT dans R

Dans ce qui suit, nous allons utiliser .alert[`mvpart` qui est archivé sur le CRAN]. Nous l'installons depuis GitHub avec le package remotes:

```r
remotes::install_github("cran/mvpart")

library(mvpart)
```

---
exclude: true

# MRT dans R

```r
# Enlever la variable “distance from source”
# Remove this chunk when LDA is back in. It is needed for the LDA (just switch next chunk to true).
env <- subset(env, select = -das)
```

```r
# Enlever la variable “distance from source”
env <- subset(env, select = -das)

# Créer l'arbre de regression multivarié
# library(mvpart)
doubs.mrt <- mvpart(as.matrix(spe.hel) ~ ., data = env,
 xv = "pick", # selection graphique intéractive
 xval = nrow(spe.hel), # nombre de validations
 xvmult = 100, # nombre de validations multiples
 which = 4, # identifier les noeuds
 legend = FALSE, margin = 0.01, cp = 0)
```

---
exclude: true

# MRT dans R: Processus de sélection

.small[
* Points verts: erreur relative
* Points bleus: erreur relative de validation croisée (CVRE)
* Point rouge: arbre avec la valeur minimale de CVRE
* Point orange: l'arbre le plus petit ayant un CVRE à 1 écart type du CVRE minimal
* Barres vertes: nombre de fois que chaque taille d'arbre a été choisie
]

---
exclude: true

# MRT dans R: Processus de sélection

.small[
* Cliquez sur le point bleu correspondant à la taille de l'arbre choisie!
* Puisqu'on ne sait pas *a priori* comment partitionner ces données, on choisira .comment[l'arbre le plus petit ayant un CVRE à 1 écart type du CVRE minimal] (i.e. le point orange).
]

---
exclude: true

# MRT dans R: Visualisation

.small[
* La matrice d'abondances est partitionnée selon un seuil d'.alert[altitude (361.5)].
  * "Barplots": abondances d'espèces incluses dans chaque groupe
* Erreur résiduelle = 0.563, alors le R2 du modèle est .alert[43.7%]
]

---
exclude: true

# MRT dans R: Comparaison d'arbres

Pour choisir un arbre, on peut aussi comparer plusieurs solutions possibles.

Par exemple, considérons une solution à 10 groupes!

.small[
* L'interprétation est .alert[plus difficile].
* Plus grand pouvoir explicatif, MAIS le pouvoir prédictif (CV Error = 0.671) ressemble à la solution précédente (CV Error = 0.673).
]

---
exclude: true

# MRT dans R: Comparaison d'arbres

Considérons une solution avec moins de groupes (4)!
.tiny[

]
.small[
* Plus facile à interpréter!
* Plus grand pouvoir explicatif .alert[(Error)] que notre solution originale
* .alert[Plus grand pouvoir prédictif] que les 2 solutions précédentes (CV Error)
]

---
exclude: true

# MRT dans R: Paramètre de complexité

Le .comment[paramètre de complexité (CP)] représente la variance expliquée par chaque noeud.

```r
doubs.mrt$cptable
```

* CP @ nsplit 0 = R2 de l'arbre au complet
* CP @ autres neouds = R2 de chaque noeud (voir sommaire complet pour la valeur de seuil de chaque noeud)

---
exclude: true

# MRT dans R: Sommaire des résultats

Pour accéder au sommaire des résultats:
.tiny[

```r
summary(doubs.mrt)
```
]

---
exclude: true

# MRT dans R: Espèces discriminantes

Pour déterminer les espèces discriminantes .alert[significatives] pour chaque groupe de sites:

```r
library(labdsv)

# Calcul d'une valeur indval pour chaque espèce
doubs.mrt.indval <- indval(spe.hel, doubs.mrt$where)

# Extraire les espèces indicatrices à chaque noeud
doubs.mrt.indval$maxcls[which(doubs.mrt.indval$pval <= 0.05)]

# Extraire leur valeur indval
doubs.mrt.indval$indcls[which(doubs.mrt.indval$pval <= 0.05)]
```
]

* Les principales espèces discriminantes au premier noeud sont TRU, VAI et ABL.
* TRU et VAI contribuent beaucoup à la branche de gauche, alors que ABL est davantage indicatrice des sites à basse altitude (<361.5m).

---
exclude: true

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

.small[
Créez un arbre de régression multivarié pour les données .comment[mite].
* Choisir l'arbre le plus petit à 1 écart type du CVRE minimal.
* Quelle est la variance totale expliquée par cet arbre?
* Combien y a-t-il de feuilles?
* Quelles sont les 3 principales espèces discriminantes?
]
 
.small[
Rappel: chargez les données!

```r
data("mite")
data("mite.env")
```

Rappel de fonctions utiles:

```r
?mvpart() # argument 'xv'!
summary()
```
]
---
exclude: true

# Défi 4: Solution

**Étape 1:** Créer un arbre de régression multivarié.

```r
mite.mrt <- mvpart(as.matrix(mite.spe.hel) ~ ., data = mite.env,
 xv = "1se",
 xval = nrow(mite.spe.hel),
 xvmult = 100,
 which = 4, legend = FALSE, margin = 0.01, cp = 0,
 prn = FALSE)
```

---
exclude: true

# Défi 4: Solution

.small[
* Quelle est la variance totale expliquée (R2) par cet arbre?
  * 1 - Error = 0.252, alors l'arbre explique .alert[25.2%] de la variation dans la matrice d'abondances.

* Combien y a-t-il de feuilles?
  * 2 feuilles
]

---
exclude: true

# Défi 4: Solution

**Étape 2**: Identifier les espèces indicatrices.

Quelles sont les .alert[espèces indicatrices] pour chaque groupe de sites?

```r
# Calcul d'une valeur indicatrice pour chaque espèce
mite.mrt.indval <- indval(mite.spe.hel, mite.mrt$where)

# Extraire les espèces indicatrices à chaque noeud
mite.mrt.indval$maxcls[which(mite.mrt.indval$pval <= 0.05)]

# Extraire leur valeur indval
mite.mrt.indval$indcls[which(mite.mrt.indval$pval <= 0.05)]
```

---

---
# Analyse discriminante linéaire (LDA)

Une LDA divise une variable réponse en groupes selon un facteur en trouvant une combinaison de variables qui donne la **meilleure séparation possible entre les groupes**. Ceci est utile parce que:
* La LDA détermine si une matrice de variables indépendantes explique bien un groupement établi *a priori*;
* Et permet donc de prédire la classification de nouvelles données!

.comment[.small[e.g. prédire l’appartenance d’une espèce de poisson à un groupe (marin vs. eau douce) selon sa morphologie]]

???
Le regroupement est effectué en maximisant la dispersion entre les groupes par rapport à la dispersion à l'intérieur des groupes.

---
# LDA dans R: Rivière Doubs

Généralement, les variables environnementales changent avec la latitude. On peut donc poser la question suivante:

--

.center[
**Si on classifie les sites de la rivière Doubs en fonction de leur latitude, à quel point les variables environnementales expliquent-elles ces groupements?**

]

---
# LDA dans R: Rivière Doubs

Commençons par charger les données spatiales des sites:

```r
# charger les données spatiales des sites Doubs:
spa <- read.csv("data/doubsspa.csv", row.names = 1)
spa$site <- 1:nrow(spa) # assigner un chiffre par site
spa <- spa[-8,] # enlever le site #8
```

Ensuite, on peut classifier les sites dans 3 groupes de latitudes:

```r
spa$group <- NA # créer colonne "group"
spa$group[which(spa$y < 82)] <- 1
spa$group[which(spa$y > 82 & spa$y < 156)] <- 2
spa$group[which(spa$y > 156)] <- 3
```

---
# LDA dans R: Rivière Doubs

Visualisons ces regroupements par latitude:

---
# LDA dans R

.alert[Note]: Normalement, nous devons vérifier que les matrices de covariance des variables explicatives sont homogènes (voir Borcard et al. 2011). Au besoin, on peut utiliser la fonction `betadisper()` du paquet `vegan`.

Pour les besoins de l'atelier, on passera directement à la LDA:

```r
# charger la libraire requise
library(MASS)

# faire la LDA
LDA <- lda(env, spa$group)
```

---
# LDA dans R

Nos sites sont maintenant réorganisés en groupes qui sont les plus distincts possibles, à partir des variables environnementales.

```r
# prédire les groupes à partir de la LDA
lda.plotdf <- data.frame(group = spa$group, lda = predict(LDA)$x)
```
 
<img src="workshop10-pres-fr_files/figure-html/unnamed-chunk-78-1.png" width="468" style="display: block; margin: auto;" />

---
# LDA dans R: Exactitude du groupement

On peut ensuite voir comment les sites sont classifiés en groupes, et si cette classification est exacte.

```r
# classification des objets en fonction de la LDA
spe.class <- predict(LDA)$class

# probabilités que les objets appartiennent à chaque groupe a posteriori
spe.post <- predict(LDA)$posterior

# tableau des classifications a priori et prédites
(spe.table <- table(spa$group, spe.class))
# spe.class
# 1 2 3
# 1 7 0 0
# 2 0 6 0
# 3 0 0 16

# proportion de classification correcte
diag(prop.table(spe.table, 1))
# 1 2 3 
# 1 1 1
```
]

Tous les sites ont été correctement classifiés dans leur groupe de latitude en fonction des variables environnementales.

---
# LDA dans R: Prédictions

On peut maintenant utiliser la LDA pour classifier de **nouveaux** sites dans les groupes de latitude.

Essayons de **prédire la classification** de cinq nouveaux sites à l'aide de la LDA:
.small[

```r
# charger les nouvelles données
classify.me <- read.csv("data/classifyme.csv", header = TRUE)
# enlever das
classify.me <- subset(classify.me, select = -das)

# prédire le groupement des nouvelles données
predict.group <- predict(LDA, newdata = classify.me)

# prédire la classification pour chaque site
predict.group$class
# [1] 1 1 1 3 3
# Levels: 1 2 3
```
]

???

Nos nouveaux sites, dans l'ordre, ont été classés dans les groupes 1, 1, 1, 3 et 3.

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

.small[Créez quatre groupes de latitude avec des étendues égales à partir des données `mite.xy`. Ensuite, faites une LDA sur les données environnementales `mite.env` des acariens (`SubsDens` et `WatrCont`).

**Quelle proportion de sites ont été classifiés correctement dans le groupe 1? Le groupe 2?**

Pour commencer, chargez les données `mite.xy`:

```r
data(mite.xy)
```

Rappel de fonctions utiles:

```r
lda()
predict()
table()
diag()
```
]

???

Si vous manquez de temps, passez directement à la prochaine diapo pour montrer les groupes de latitudes.

---
# Défi 4: Solution

**Étape 1:** Créer quatre groupes de latitude avec des étendues égales.

```r
# numéroter les sites
mite.xy$site <- 1:nrow(mite.xy)

# trouver une étendue égale de latitudes par groupe
(max(mite.xy[,2])-min(mite.xy[,2]))/4
# [1] 2.4

# classifier les sites dans 4 groupes de latitude
mite.xy$group <- NA # nouvelle colonne "group"
mite.xy$group[which(mite.xy$y < 2.5)] <- 1
mite.xy$group[which(mite.xy$y >= 2.5 & mite.xy$y < 4.9)] <- 2
mite.xy$group[which(mite.xy$y >= 4.9 & mite.xy$y < 7.3)] <- 3
mite.xy$group[which(mite.xy$y >= 7.3)] <- 4
```

---
# Défi 4: Solution

**Étape 2:** Faire la LDA.
.small[

```r
LDA.mite <- lda(mite.env[,1:2], mite.xy$group)
```
]

**Étape 3:** Vérifier l'exactitude du groupement.

```r
# classification des objects en fonction de la LDA
mite.class <- predict(LDA.mite)$class
# tableeau de classifications (prior versus predicted)
(mite.table <- table(mite.xy$group, mite.class))
# mite.class
# 1 2 3 4
# 1 9 4 2 0
# 2 2 11 4 0
# 3 1 2 14 2
# 4 0 0 3 16
# proportion de classifications exactes
diag(prop.table(mite.table, 1))
# 1 2 3 4 
# 0.6000000 0.6470588 0.7368421 0.8421053
```
]

---
# Défi 4: Solution

On peut maintenant répondre à la question du défi avec cette partie du code:

```r
# proportion of correct classification
diag(prop.table(mite.table, 1))
#         1         2         3         4 
# 0.6000000 0.6470588 0.7368421 0.8421053
```

**Alors, quelle proportion de sites ont été classifiés correctement dans le groupe 1? Le groupe 2?**

* .alert[60%] des sites ont été classifiés correctement dans le groupe 1, et .alert[64.7%] dans le groupe 2.

---

# Merci d'avoir participé à cet atelier!

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