Série d'ateliers en R du CSBQ

Cette série de 10 ateliers guide les participants à travers les étapes requises afin de maîtriser le logiciel R pour une grande variété d’analyses statistiques pertinentes en recherche en biologie et en écologie. Ces ateliers en libre accès ont été créés par des membres du CSBQ à la fois pour les membres du CSBQ et pour la grande communauté d’utilisateurs de R.

Le contenu de cet atelier a été révisé par plusieurs membres du CSBQ. Si vous souhaitez y apporter des modifications, veuillez SVP contacter les coordonnateurs actuels de la série, listés sur la page d'accueil

AVIS IMPORTANT: MISES À JOUR MAJEURES

Mise à jour de mars 2021: Ce wiki n'est plus activement développé ou mis à jour. Le matériel mis à jour pour la série d'ateliers R du CSBQ est maintenant hébergé sur la page GitHub des ateliers R du CSBQ.

Le matériel disponible inclut;

1. La présentation Rmarkdown pour cet atelier;
2. Un script R qui suit la présentation - *en construction*.
3. Le matériel écrit qui accompagne la présentation en format bookdown - *en construction*.

Atelier 10: Analyses multivariées avancées

Développé par: Monica Granados, Emmanuelle Chrétien, Bérenger Bourgeois, Amanda Winegardner and Xavier Giroux-Bougard. (Le matériel des scripts R est adapté de: Borcard, Gillet & Legendre (2011). Numerical Ecology with R. Springer New York.)

Résumé: Durant cet atelier, vous apprendrez à réaliser des analyses multivariées avancées sur des données de communauté. Cet atelier se concentre sur les méthodes sous contraintes, telles que l'analyse canonique de redondances (RDA), l'arbre de régression multivarié (MRT) et l'analyse discriminante linéaire (LDA) afin d'explorer comment les variables environnementales peuvent expliquer les patrons de composition en espèces à travers différents sites.

Lien vers la nouvelle présentation Rmarkdown

Lien vers la présentation Prezi associée : Prezi

Téléchargez le code R, les librairies et données requises pour cet atelier:

• Code R
• données DoubsEnv
• données DoubsSpe
• données test pour l'analyse linéaire discriminante
• librairie mvpart
• librairie MVPARTwrap
• librairie rdaTest

Assurez-vous d'importer les librairies suivantes dans R Studio (procédure fournie dans le code R):

• vegan (pour analyses multivariées)
• labdsv (pour identification d'espèces indicatrices pour l'arbre de régression multivarié)
• plyr (classification pour l'analyse linéaire discriminante)
• MASS (pour l'analyse linéaire discriminante)
• mvpart
• MVPARTwrap
• rdatest

| Installer et importer les librairies

install.packages("vegan")
install.packages("mvpart")
install.packages("labdsv")
install.packages("plyr")
install.packages("MASS")
 
# Pour les deux librairies suivantes, importez le fichier d'archivage donné sur la page wiki.
# Pour ce faire, cliquez sur l'onglet "Packages" de la section en bas à droite sur R Studio
# Cliquez sur "Install Packages"
# Choisissez d'installer depuis "Package Archive file" et importez les deux fichiers d'archivage
install.packages("MVPARTwrap")
install.packages("rdaTest")
 
library(vegan)
library(mvpart)
library(MVPARTwrap)
library(rdaTest)
library(labdsv)
library(plyr)
library(MASS)

Introduction aux analyses multivariées avancées

L'atelier précédent donnait un aperçu des analyses multivariées de base:

• comment choisir les mesures de distances et les transformations appropriées selon le type de données
• groupement hiérarchique
• ordinations sans contraintes (analyses en composantes principales, analyse en coordonnées principales, analyse de correspondances, positionnement multidimensionnel non-métrique)

En se basant sur ces acquis, le présent atelier se concentrera sur les analyses sous contraintes. Toutes les méthodes vues lors du précédent atelier ont permis de relever des tendances dans la structure de communautés d'espèces ou des descripteurs par rapport à des sites, mais pas d'explorer comment les variables environnementales pouvaient expliquer ces tendances. Avec des analyses sous contraintes telles l'analyse canonique de redondances (ACR), l'analyse linéaire discriminante (ALD) et l'arbre de régression multivarié (ARM), il sera possible de décrire et de prédire les relations entre la structure des communautés et les variables environnementales.

Introduction aux données

Encore une fois, nous utiliserons les données de la rivière Doubs. “DoubsSpe.csv” est une matrice de données d'abondance d'espèces de communautés de poissons dans laquelle la première colonne contient les noms des sites de 1 à 30 et les colonnes subséquentes correspondent aux différentes espèces de poissons. “DoubsEnv.csv” est une matrice de données environnementales pour les mêmes sites. La première colonne contient donc les noms des sites de 1 à 30 et les colonnes suivantes les mesures de 11 variables abiotiques. Notez que les données utilisées pour les analyses d'ordination sont généralement en format long (en anglais).

| Chargez les données

#Matrice d'abondances d'espèces: “DoubsSpe.csv”
spe<- read.csv(file.choose(), row.names=1)
spe<- spe[-8,] #Pas d'espèces dans le site 8, supprimer site 8. Exécuter cette ligne une seule fois 
 
#Matrice de données environnementales: “DoubsEnv.csv”
env<- read.csv(file.choose(), row.names=1)
env<- env[-8,] #Supprimer le site 8 puisqu'on l'a supprimé de la matrice d'abondance. N'exécuter qu'une seule fois. 

1. Exploration et préparation des données

1.1 Données d'abondances d'espèces

Nous pouvons utiliser les fonctions de résumé R pour explorer les données “Spe” (données d'abondances de poissons) et découvrir les caractéristiques telles que les dimensions de la matrice, les noms des colonnes et les statistiques descriptives de ces colonnes (révision de l'atelier 2).

| Explorez DoubsSpe

names(spe) #voir les noms des colonnes spe
dim(spe) #dimensions de spe; nombre de lignes et de colonnes
str(spe) #la structure interne de la matrice
head(spe) #les premières lignes
summary(spe) #statistiques descriptives  

Regardez la distribution des espèces.

| La distribution des espèces (DoubsSpe)

(ab <- table(unlist(spe))) #en utilisant les parenthèses de cette façon, la sortie s'affichera immédiatement dans la console
 
barplot(ab, las=1, xlab="Abundance class", ylab="Frequency", col=grey(5:0/5))

Il y a une grande fréquence de zéros dans les données d'abondance.

Calculez le nombre d'absences.

| Absences dans les données de poissons

sum(spe==0) 

Regardez la proportion de zéros dans les données d'abondances de poissons.

| Proportion de zéros

sum(spe==0)/(nrow(spe)*ncol(spe))

La proportion de zéros dans la matrice est de ~0.5. C'est élevé mais pas inhabituel pour des données d'abondances d'espèces. Cependant, pour éviter que les double zéros soient considérés comme une similarité entre sites, nous appliquerons une transformation aux données d'abondances d'espèces.Legendre et Gallagher (2001) ont proposé cinq transformations appropriées aux données d'abondances d'espèces pour les analyses multivariées et quatre d'entre elles sont disponibles dans la librairie vegan sous la fonction decostand().

La transformation d'Hellinger sera appliquée aux données d'abondances de poissons. Elle exprime l'abondance en tant que la racine carrée de l'abondance pondérée sur l'abondance totale à chaque site (Borcard et al. 2011).

| transformation de Hellinger

spe.hel <- decostand(spe, method="hellinger") # vous pouvez également utiliser method="hell" 

1.2 Données environnementales

Explorez les données environnementales et détectez les colinéarités.

| Exploration de la matrice Doubs Env

names(env)
dim(env)
str(env)
head(env)
summary(env)
pairs(env, main="Bivariate Plots of the Environmental Data" ) 

Dans ce cas, les données environnementales sont toutes dans des unités différentes et doivent donc être standardisées avant de calculer les mesures de distance utilisées pour effectuer la plupart des analyses d'ordination. La standardisation des données (11 variables) peut être effectuée en utilisant la fonction decostand () de vegan.

| Decostand

env.z <- decostand(env, method="standardize")
apply(env.z, 2, mean) # les données sont centrées (moyenne~0)
apply(env.z, 2, sd)   # les données sont réduites (écart type =1)

2. Analyses canoniques

Décrites à l’origine par Rao (1964), les analyses canoniques rassemblent de nombreuses méthodes statistiques combinant les concepts d’ordination et de régression et partageant un but commun, à savoir identifier les relations entre un ensemble de variables réponse (matrice Y décrivant généralement la composition en espèces de différentes communautés) et un ensemble de variables explicatives (matrice X contenant le plus souvent des variables environnementales). Les analyses canoniques permettent de tester des hypothèses écologiques relatives aux déterminants environnementaux de la composition en espèces des communautés. Parmi l’ensemble de méthodes d’analyses canoniques existantes, nous insisterons ici principalement sur l’analyse canonique de redondance (RDA).

2.1 Analyse canonique de redondances (ACR ou RDA)

L’analyse canonique de redondance correspond à une extension directe de la régression multiple puisqu’elle modélise l’effet d’une matrice X de variables explicatives (n x p) sur une matrice Y de variables réponses (n x m) en effectuant une ordination de Y afin d’obtenir des axes d’ordination qui correspondent à des combinaisons linéaires des variables de la matrice X. Dans une RDA, les axes d’ordination sont calculés à partir d’une PCA de la matrice Yfit, calculée en ajustant les variables Y sur les variables X par régression linéaires multivariée. Notons que les variables explicatives X peuvent être quantitatives, qualitatives ou binaires. Avant d’effectuer une RDA, les variables explicatives de X doivent être centrées, standardisées (si les variables explicatives sont d’unités différentes), transformées (afin de limiter l’asymétrie des variables explicatives) ou normalisées (afin de linéariser leurs relations) suivant les mêmes principes que dans une PCA. La colinéarité entre variables explicatives doit également être minimale pour effectuer une RDA. Afin d’obtenir le meilleur modèle de RDA, les variables explicatives peuvent êtres sélectionnées par sélections progressives ou régressives afin d’éliminer les variables explicatives non significatives.

Le calcul d’une RDA s’articule en deux étapes Dans une première étape, la matrice de valeurs ajustées Yfit est calculée via l’équation linéaire :

Yfit = X[X'X]-1 [X'Y]

Dans une seconde étape, une PCA de la matrice Yfit est implémentée (voir équations sur la figure ci-dessous) afin de calculer ses valeurs propres et vecteurs propres ainsi que la matrice Z contenant les axes d’ordination. Ces axes correspondent aux combinaisons linéaires des variables explicatives de X. La linéarité des combinaisons de variables X est une propriété fondamentale de la RDA. Lors de l’analyse de composition en espèces des communautés, ces axes d’ordination sont souvent interprétés comme des gradients environnementaux complexes.

Une RDA est un processus en deux étapes (d’après Legendre et Legendre, 2012)

Plusieurs valeurs statistiques peuvent être extraites de la RDA, en particulier :

- Le R² mesure la force de la relation canonique entre Y et X en calculant la proportion de la variation de Y expliquée par les variables X, - Le R² ajusté mesure également la force de la relation entre Y et X, mais applique une correction du R² afin de tenir compte du nombre de variables explicatives, - La statistique F correspond à un test global de la significativité de la RDA en comparant le modèle étudié à un modèle nul.

Ce test est basé sur l’hypothèse nulle selon laquelle la force de la relation calculé par le R² n’est pas supérieure à la valeur qui serait obtenue pour des matrices X et Y de même taille sans aucune relation statistique. Notons que la statistique de F peut également être utilisée pour tester la significativité de chaque axe d’ordination de manière séquentielle.

Avec R, la RDA peut être calculée en utilisant la fonction rda du package vegan, comme suit:

| effectuer une RDA avec rda() de vegan

#Préparer les données
env.z <- subset(env.z, select = -das) # Enlever la variable "distance from the source" 
 
#Faire la RDA
?rda
spe.rda <- rda(spe.hel~., data=env.z)
 
### Extraire les résultats
summary(spe.rda, display=NULL) 
 
#Les résultats peuvent être vus en utilisant summary():
summary(spe.rda, display=NULL) #display = NULL optional  

Le sommaire ressemble à ceci:

Ces résultats incluent la proportion de la variance de Y expliquée par les variables X (proportion contrainte ou “constrained proportion” 72.71% ici), la variance de Y non expliqué par X (proportion non contrainte “unconstrained proportion”, 27.29% ici) et résument les valeurs propres, les proportions expliquées et les proportions cumulées de variance pour chaque axe d’ordination.

Afin de sélectionner les variables explicatives significatives, une sélection progressive peut être effectuée via la fonction ordiR2step de vegan (ou la function forward.sel du package packfor):

| ordiR2step() pour la sélection progressive

?ordiR2step
ordiR2step(rda(spe.hel~1, data=env.z), scope= formula(spe.rda), direction= "forward", R2scope=TRUE, pstep=1000)
env.signif <- subset(env.z, select = c("alt", "oxy", "dbo"))

Dans ce cas, trois variables sont retenues par la sélection progressive, soit les variables alt, oxy et dbo. Ces trois variables peuvent alors être placées dans un nouveau tableau de données env.signif afin d’effectuer une nouvelle RDA contenant uniquement les variables X significatives :

| RDA avec seulement les variables significatives

spe.rda.signif <- rda(spe.hel~., data=env.signif)
summary(spe.rda.signif, display=NULL)

Les variables explicatives expliquent désormais 59% de la variance de Y.

Le R² ajusté de cette RDA est calculé à partir de la fonction RsquareAdj :

| R2 ajusté

(R2adj <- RsquareAdj(spe.rda.signif)$adj.r.squared)

La significativité globale de cette RDA et celle de chaque axe d’ordination peuvent être testées en utilisant la fonction anova (ceci est différent de la sélection des variables significatives. Ici, nous testons si les axes de la RDA sont significatifs):

| anova.cca pour tester la significativité des axes

?anova.cca
anova.cca(spe.rda.signif, step=1000)
anova.cca(spe.rda.signif, step=1000, by="axis")
#Le modèle global est significatif ainsi que les trois axes canoniques.

Afin de visualiser les résultats de la RDA, il est possible de tracer des triplots en utilisant la fonction plot(). Comme pour la PCA, en choisissant le cadrage (scaling) 1, les distances entre les objets sont des approximations de leurs distances euclidiennes, alors qu'en cadrage 2 les angles entre les variables X et Y reflètent leur corrélation. Les triplots en cadrage 1 peuvent donc être utilisés pour interpréter les distances entre objets alors que les triplots en cadrage 2 permettent d’interpréter les relations entre variables X et Y.

Pour obtenir un triplot scaling 1 d’une RDA, le code suivant peut être employé:

| représentations graphiques des RDAs

#Pour faire des graphiques rapidement et facilement
windows()
plot(spe.rda.signif, scaling=1, main="Triplot RDA (scaling 1)")
windows()
plot(spe.rda.signif, scaling=2, main="Triplot RDA (scaling 2)")
 
#Code avancé pour graphiques plus attrayants
#cadrage 1
windows()
plot(spe.rda.signif, scaling=1, main="Triplot RDA - scaling 1", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
points(scores(spe.rda.signif, display="sites", choices=c(1,2), scaling=1),
       pch=21, col="black", bg="steelblue", cex=1.2)
arrows(0,0,
       scores(spe.rda.signif, display="species", choices=c(1), scaling=1),
       scores(spe.rda.signif, display="species", choices=c(2), scaling=1),
       col="black",length=0)
text(scores(spe.rda.signif, display="species", choices=c(1), scaling=1),
     scores(spe.rda.signif, display="species", choices=c(2), scaling=1),
     labels=rownames(scores(spe.rda.signif, display="species", scaling=1)),
     col="black", cex=0.8)    
arrows(0,0,
       scores(spe.rda.signif, display="bp", choices=c(1), scaling=1),
       scores(spe.rda.signif, display="bp", choices=c(2), scaling=1),
       col="red")
text(scores(spe.rda.signif, display="bp", choices=c(1), scaling=1)+0.05,
     scores(spe.rda.signif, display="bp", choices=c(2), scaling=1)+0.05,
     labels=rownames(scores(spe.rda.signif, display="bp", choices=c(2), scaling=1)),
     col="red", cex=1) 
 
#cadrage 2
windows()
plot(spe.rda.signif, scaling=2, main="Triplot RDA - scaling 2", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
points(scores(spe.rda.signif, display="sites", choices=c(1,2), scaling=2),
       pch=21, col="black", bg="steelblue", cex=1.2)
arrows(0,0,
       scores(spe.rda.signif, display="species", choices=c(1), scaling=2)*2,
       scores(spe.rda.signif, display="species", choices=c(2), scaling=2)*2,
       col="black",length=0)
text(scores(spe.rda.signif, display="species", choices=c(1), scaling=2)*2.1,
     scores(spe.rda.signif, display="species", choices=c(2), scaling=2)*2.1,
     labels=rownames(scores(spe.rda.signif, display="species", scaling=2)),
     col="black", cex=0.8)    
arrows(0,0,
       scores(spe.rda.signif, display="bp", choices=c(1), scaling=2),
       scores(spe.rda.signif, display="bp", choices=c(2), scaling=2),
       col="red")
text(scores(spe.rda.signif, display="bp", choices=c(1), scaling=2)+0.05,
     scores(spe.rda.signif, display="bp", choices=c(2), scaling=2)+0.05,
     labels=rownames(scores(spe.rda.signif, display="bp", choices=c(2), scaling=2)),
     col="red", cex=1)  

Les triplots finaux:

Défi 1: Effectuer la RDA de l’abondance des espèces du tableau d'acariens en fonction des variables environnementales mite.env

| Chargez les données d'acariens

#Ces données sont disponibles dans la librairie vegan
data(mite) 
mite.spe<-mite
mite.spe.hel <- decostand(mite.spe, method="hellinger")
 
data(mite.env)

Quelles sont les variables explicatives significatives ? Quel pourcentage de variance est expliqué par les variables significatives ? Quels sont les axes significatifs ? Quels groupes de sites pouvez-vous identifier ? Quelles espèces sont liées à chaque groupe de sites ?

Défi 1: Solution

2.2 RDA partielle

La RDA partielle est un cas particulier de la RDA dans lequel une matrice Y de variables réponses est ajustée à une matrice de variables explicatives X, en présence d’une matrice W de variables explicatives additionnelles, appelées co-variables. Comme pour une régression linéaire partielle, l’effet de X sur Y est ajusté pour tenir compte de l’effet des co-variables W. Pour cela, une RDA de l’effet des co-variables W sur Y est d’abord effectuée. Les résidus de cette RDA sont ensuite extraits, c’est-à-dire une matrice Yres|W contenant les variables réponses Y dans laquelle l’effet des co-variables W a été supprimé. La RDA partielle correspond alors à la RDA de Yres|W en fonction de X. Toutes les valeurs statistiques présentées auparavant pour la RDA (R2, R2 ajusté et statistique F) s’appliquent également à la RDA partielle.

La RDA partielle est donc un outil statistique particulièrement performant pour évaluer l’effet de variables environnementales sur la composition en espèces des communautés tout en prenant en compte la variation d’abondance des espèces due à d’autres variables environnementales de moindre intérêt. Cette méthode peut également être utilisée pour contrôler des effets linéaires connus, isoler l’effet d’une variables explicative ou pour analyser des échantillons appariés. Dans l’exemple ci-dessous, nous évaluerons l’effet de la physico-chimie de l’eau sur l’abondance des espèces de poisson en supprimant l’effet de la physiographie.

Avec R, la RDA partielle est effectuée en utilisant le fonction rda et en ajoutant un argument correspondant aux co-variables :

| RDA partielle avec rda()

#Divisez le tableau de données environnementales en deux: 
envtopo <- env[, c(1:3)] # Physiographie : tableau 1
names(envtopo)
envchem <- env[, c(4:10)] # Physico-chimie de l'eau : tableau 2
names(envchem)
 
#FAire la RDA partielle
spechem.physio=rda(spe.hel, envchem, envtopo)
summary(spechem.physio, display=NULL)
#ou
spechem.physio2=rda(spe.hel ~ pH + dur + pho + nit + amm + oxy + dbo 
                    + Condition(alt + pen + deb), data=env)
 
#Extraire les résultats
summary(spechem.physio, display=NULL)
 
#Calcul du R2 ajusté de la RDA partielle
(R2adj <- RsquareAdj(spechem.physio)$adj.r.squared)
 
#Test de la significativité des axes
anova.cca(spechem.physio, step=1000)
anova.cca(spechem.physio2, step=1000, by="axis")
 
#Faire les triplots
#Cadrage 1
windows(title="Partial RDA scaling 1")
plot(spechem.physio, scaling=1, main="Triplot partial RDA - scaling 1", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
points(scores(spechem.physio, display="sites", choices=c(1,2), scaling=1),
       pch=21, col="black", bg="steelblue", cex=1.2)
arrows(0,0,
      scores(spechem.physio, display="species", choices=c(1), scaling=1),
      scores(spechem.physio, display="species", choices=c(2), scaling=1),
      col="black")
text(scores(spechem.physio, display="species", choices=c(1), scaling=1),
     scores(spechem.physio, display="species", choices=c(2), scaling=1),
     labels=rownames(scores(spechem.physio, display="species", scaling=1)),
     col="black", cex=0.8)  
arrows(0,0,
      scores(spechem.physio, display="bp", choices=c(1), scaling=1),
      scores(spechem.physio, display="bp", choices=c(2), scaling=1),
      col="red")
text(scores(spechem.physio, display="bp", choices=c(1), scaling=1)+0.05,
     scores(spechem.physio, display="bp", choices=c(2), scaling=1)+0.05,
     labels=rownames(scores(spechem.physio, display="bp", choices=c(2), scaling=1)),
     col="red", cex=1) 
 
#Cadrage 2
windows(title="Partial RDA scaling 2")
plot(spechem.physio, scaling=2, main="Triplot partial RDA - scaling 2", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
points(scores(spechem.physio, display="sites", choices=c(1,2), scaling=2),
       pch=21, col="black", bg="steelblue", cex=1.2)
arrows(0,0,
      scores(spechem.physio, display="species", choices=c(1), scaling=2),
      scores(spechem.physio, display="species", choices=c(2), scaling=2),
      col="black")
text(scores(spechem.physio, display="species", choices=c(1), scaling=2),
     scores(spechem.physio, display="species", choices=c(2), scaling=2),
     labels=rownames(scores(spechem.physio, display="species", scaling=2)),
     col="black", cex=0.8)  
arrows(0,0,
       scores(spechem.physio, display="bp", choices=c(1), scaling=2),
       scores(spechem.physio, display="bp", choices=c(2), scaling=2),
       col="red")
text(scores(spechem.physio, display="bp", choices=c(1), scaling=2)+0.05,
     scores(spechem.physio, display="bp", choices=c(2), scaling=2)+0.05,
     labels=rownames(scores(spechem.physio, display="bp", choices=c(2), scaling=2)),
     col="red", cex=1)  

Cette RDA partielle est significative (p<0.001) de même que les deux premiers axes d’ordination. La physico-chimie de l’eau explique 31.89% de l’abondance des espèces de poissons tandis que les co-variables physiographiques expliquent 41.53% de la variation en abondances des poissons. La variance non expliquée s’élève à 26.59%. Le R2 ajusté de cette RDA partielle est de 24.13%.

Défi 2

Effectuer la RDA partielle de l’abondance des espèces du tableau mite en fonction des variables environnementales en supprimant l’effet dû aux paramètres du substrat (SubsDens, WaterCont and Substrate) Le modèle est-il significatif ? Quels sont les axes significatifs ? Interpréter le triplot obtenu.

Défi 2: Solution

2.3 Partitionnement de la variation par RDA partielles

Le partitionnement de la variation est un type d’analyse qui combine à la fois la RDA et la RDA partielle pour diviser la variation d’une matrice de variable réponse en deux, trois ou quatre jeux de données explicatives. Le résultat d’un partitionnement de la variation est généralement représenté par un diagramme de Venn sur lequel sont annotés les pourcentages de variance expliquée par chacun des jeux de données explicatives (ou par leurs interactions). Dans le cas de deux jeux de données explicatives X et W :

- La fraction a + b + c correspond à la variance expliquée par les deux jeux de données calculée par une RDA de Y en fonction de X + W, - La fraction d est la variance non expliquée par les deux jeux de données calculée par la même RDA que précédemment, - La fraction a est la variance expliquée par les variables X uniquement, calculée par une RDA partielle Y en fonction de X en utilisant les variables W comme co-variables, - La fraction c est la variance expliquée par les variables W uniquement, calculée par une RDA partielle Y en fonction de W en utilisant les variables X comme co-variables, - La fraction b correspond à l’interaction des deux jeux de données calculée par soustraction, i.e. b = (a + b + c) – a – b.

Le diagramme de Venn représente le partitionnement de la variation de Y en fonction de deux jeux de données explicatives X et W (d’après Legendre et Legendre, 2012)

Le partitionnement de la variation est donc une analyse toute indiquée pour expliquer l’abondance des espèces dans une communauté en fonction de différents types de variables explicatives, par exemple en présence de variables abiotiques vs biotiques, de variables locales vs à large échelle, etc. Dans l’exemple suivant, nous partitionnerons la variation d’abondance des espèces de poissons en fonction des variables physico-chimiques et des variables physiographiques.

Dans R, le partitionnement de la variation est calculé en utilisant la fonction varpart. Les diagrammes de Venn peuvent être tracés à l’aide de la fonction plot.

| varpart()

?varpart
vegandocs("partitioning.pdf")
 
#Partitionnement de la variation avec toutes les variables explicatives
spe.part.all <- varpart(spe.hel, envchem, envtopo)
spe.part.all
windows(title="Variation partitioning - all variables")
plot(spe.part.all, digits=2)  

Le sommaire ressemble à ceci:

Dans ce cas, les variables physico-chimiques expliquent 24.10% de la composition en espèces de poissons, les variables physiographiques 11.20% de la variation et l’interaction de ces deux types de variables explique 23.30% de la variation. Notons que la fonction varpart permet aussi d’identifier les fractions dont la significativité peut être testée à l’aide de la fonction anova.cca.

Il est également possible d’effectuer un partitionnement de la variation à partir de jeux de données explicatives contenant uniquement des variables significatives :

| varpart avec variables explicatives

#RDA des variables physico-chimiques
spe.chem <- rda(spe.hel~., data=envchem)
 
#Sélection des variables significatives
R2a.all.chem <- RsquareAdj(spe.chem)$adj.r.squared
ordiR2step(rda(spe.hel~1, data=envchem), 
           scope= formula(spe.chem), direction= "forward", R2scope=TRUE, pstep=1000)
names(envchem)
(envchem.pars <- envchem[, c( 4, 6, 7 )])
 
#RDA avec les autres variables significatives
spe.topo <- rda(spe.hel~., data=envtopo)
R2a.all.topo <- RsquareAdj(spe.topo)$adj.r.squared
ordiR2step(rda(spe.hel~1, data=envtopo), 
           scope= formula(spe.topo), direction= "forward", R2scope=TRUE, pstep=1000)
names(envtopo)
envtopo.pars <- envtopo[, c(1,2)]
 
#Varpart
spe.part <- varpart(spe.hel, envchem.pars, envtopo.pars)
windows(title="Variation partitioning - parsimonious subsets")
plot(spe.part, digits=2)
 
#Tests de significativité des fractions
anova.cca(rda(spe.hel, envchem.pars), step=1000) # Test of fractions [a+b]
anova.cca(rda(spe.hel, envtopo.pars), step=1000) # Test of fractions [b+c]
env.pars <- cbind(envchem.pars, envtopo.pars) 
anova.cca(rda(spe.hel, env.pars), step=1000) # Test of fractions [a+b+c]
anova.cca(rda(spe.hel, envchem.pars, envtopo.pars), step=1000) # Test of fraction [a]
anova.cca(rda(spe.hel, envtopo.pars, envchem.pars), step=1000) # Test of fraction [c]

Les variables physico-chimiques expliquent maintenant 25.30% de la variation d’abondances des poissons, les variables physiographiques 14.20% de la variation et l’interaction de ces deux types de variables 19.60% de la variation. Toutes ces fractions sont significatives (p<0.001).

Défi 3 Effectuer le partitionnement de la variation de l’abondance des espèces du tableau mite en utilisant un premier jeu de variables explicatives relatives au substrat (SubsDens, WaterCont and Substrate) et un second jeu de variables explicatives contenant les autres variables significatives (Shrud and Topo) Quelle est la proportion de variance expliquée par chaque jeu de données ? Quelles sont les fractions significatives ?

Défi 3: Solution

3. Arbre de régression multivarié

L'arbre de régression multivarié (ARM ou MRT) est une méthode de groupement hiérarchique sous contrainte (De’ath 2002). Le MRT fait le partitionnement d'une matrice réponse quantitative (ex. données d'abondances d'espèce) sous la contrainte d'une matrice de variables explicatives. Celle-ci détermine les points de séparation des données d'abondance en différents groupes, de manière à minimiser la somme des carrés des écarts intra-groupes (comme pour le groupement hiérarchique de Ward). La RDA et le MRT sont toutes deux des méthodes de régression, la première expliquant la structure globale des relations par un modèle linéaire, la dernière mettant davantage en lumière les structures locales et les interactions entre variables en produisant un arbre.

Avantages du MRT par rapport à la RDA:

• n'assume pas de relation linéaire entre les abondances d'espèces et les variables environnementales (quantitatives ou qualitatives),
• la méthode est robuste en présence de valeurs manquantes,
• la méthode est robuste en présence de colinéarité entre les descripteurs,
• les MRTs ne nécessitent pas la transformation des variables explicatives (les valeurs brutes peuvent être utilisées).
• l'arbre obtenu est facile à interpréter pour un public scientifique ou néophyte.

Le MRT divise les données en groupes ayant des compositions en espèce semblables et caractérisés par des variables environnementales. La méthode implique deux volets s'effectuant en parallèle: 1) la construction de l'arbre et 2) la sélection de la partition finale optimale par validation croisée. La fonction mvpart() de la librairie mvpart {} permet de créer les arbres de régression.

Vocabulaire lié aux MRTs:

Feuille: Groupe terminal de sites

Noeud: Point où les données se divisent en deux groupes. Chaque noeud est caractérisé par une valeur d'une variable explicative.

Branche: Chaque lignée formée par un noeud.

1- Partitionnement des données sous contrainte (construction de l'arbre)

Premièrement, la méthode calcule toutes les partitions des sites en deux groupes. Pour chaque variable environnementale quantitative, les sites seront classés en ordre croissant des valeurs; pour chaque variable qualitative (ou catégorique), les sites seront classés par niveaux. La méthode divise les données après le premier objet, après le second, et ainsi de suite et calcule à chaque fois la somme des carrés des écarts intra-groupes de la matrice réponse. La méthode choisira la partition qui minimisera la somme des carrés des écarts intra-groupes et le point de division défini par une valeur seuil d'une variable environnementale. Ces étapes seront répétées dans les deux groupes formés précédemment, jusqu'à ce que tous les objets forment leur propre groupe. En d'autre mots, jusqu'à ce que chaque feuille de l'arbre de contienne qu'un seul objet.

2- Validation croisée et élagage de l'arbre

La fonction effectue également une validation croisée et identifie l'arbre ayant le meilleur pouvoir prédictif. La validation croisée s'effectue en utilisant une partie des données pour construire l'arbre et le reste des données est classé dans les groupes créés. Dans un arbre ayant un bon pouvoir prédictif, les objets sont assignés aux groupes appropriés. L'erreur relative de validation croisée (ERVC ou CVRE) mesure l'erreur de prédiction. Sans validation croisée, le nombre de partitions retenu serait celui minimisant la variance non expliquée par l'arbre (i.e. l'erreur relative: la somme des carrés des écarts intra-groupes de toutes les feuilles divisée par la somme de carrée des écarts de toutes les données). Cette solution maximise le R2 et on obtiendrait donc un arbre explicatif plutôt que prédictif.

Créons un arbre de régression multivarié sur les données de la rivière Doubs.

| mvpart()

?mvpart
 
#Preparer les données: enlever la variable “distance from source”
env <- subset(env, select = -das)
 
# Créer l'arbre
doubs.mrt<-mvpart(data.matrix(spe.hel)~.,env,
legend=FALSE,margin=0.01,cp=0,xv="pick",xval=nrow(spe.hel),xvmult=100,which=4)

Vous devrez maintenant sélectionner la taille de l'arbre voulue sur le graphique qui s'affichera puisque nous avons utilisé l'argument xv=“pick”. En d'autres mots, il faut élaguer l'arbre. Un arbre complètement résolu n'est pas la solution idéale. Il est plus intéressant d'obtenir un arbre contenant seulement des groupes qu'il est possible d'interpréter. Il est possible d'avoir une idée à priori du nombre de groupes voulu, d'opter pour la taille minimisant la CVRE ou d'opter pour la taille minimale dont la CVRE est à un écart-type de la CVRE minimale (Breiman et al. 1984).

Le graphique montre l'erreur relative (RE, en vert) et l'erreur relative de validation croisée (en bleu) d'arbres de tailles croissantes. Le point rouge indique la solution avec la valeur minimale de CVRE et le point orange montre l'arbre le plus petit dont la valeur de CVRE est à 1 écart type de de la valeur CVRE minimale. Breiman et al. (1984) suggèrent de choisir cette dernière option car cet arbre a à la fois une erreur relative de validation croisée près de la plus faible et il contient un nombre restreint de groupe, ce qui en fait un choix parcimonieux. Les barres vertes en haut du graphique indiquent le nombre de fois que chaque taille d'arbre a été choisi durant le processus de validation croisée.

Ce graphique est interactif. Il faudra donc cliquer sur le point bleu correspondant à la taille de l'arbre choisie. Après avoir cliqué, l'arbre de régression multivarié s'affichera. En cliquant sur le point orange, la figure suivante apparaîtra.

Des statistiques sont affichées au bas de l'arbre: l'erreur résiduelle (dont la réciproque est le R2 du modèle, dans le cas présent 43.7%), l'erreur de validation croisée et l'erreur type. Cet arbre n'est constitué que de deux branches séparées par un noeud. Ce noeud divise les données en deux groupes à la valeur seuil d'altitude de 361.5m.

Sous chaque feuille, on retrouve un petit diagramme à bandes montrant les abondances des espèces des sites retrouvés dans la branche, le nombre de sites et l'erreur relative.

Nous pouvons comparer cet arbre avec la solution minimisant le CVRE (contenant 10 feuilles), ou avec une solution intermédiaire (ex. avec 4 feuilles).

| Comparaison des arbres

# En utilisant le critère CVRE
doubs.mrt.cvre<-mvpart(data.matrix(spe.hel)~.,env,
legend=FALSE,margin=0.01,cp=0,xv="pick",xval=nrow(spe.hel),xvmult=100,which=4)
 
# En choisissant un arbre avec 4 feuilles
doubs.mrt.4<-mvpart(data.matrix(spe.hel)~.,env,
legend=FALSE,margin=0.01,cp=0,xv="pick",xval=nrow(spe.hel),xvmult=100,which=4)

L'arbre de 10 feuilles a un très grand pouvoir explicatif, mais son pouvoir prédictif n'est que marginalement plus élevé que l'arbre à 2 feuilles. L'arbre à 4 feuilles semble être un compromis entre les deux.

Le sommaire de la fonction contient beaucoup d'informations.

| sommaire du MRT

summary(doubs.mrt)

CP est une abréviation pour “complexity parameter” (paramètre de complexité) qui est en fait l'équivalent de la variance expliquée à chaque noeud. Le CP au nombre de divisions 0 est en quelque sorte le R2 de l'arbre. Le sommaire indique ensuite, pour chaque noeud, la meilleure valeur seuil à utiliser pour diviser les données. Bien que très informatif, ce sommaire est dense. La fonction MRT() de la librairie MVPARTwrap donne accès à un sommaire plus détaillé, permet d'extraire davantage d'informations et identifie les espèces discriminantes pour chaque branche.

| Trouver des espèces discriminantes et indicatrices

# Trouver les espèces discriminantes de l'arbre de régression
doubs.mrt.wrap<-MRT(doubs.mrt,percent=10,species=colnames(spe.hel))
summary(doubs.mrt.wrap)
 
# Extraire les p-values des valeurs indval
doubs.mrt.indval<-indval(spe.hel,doubs.mrt$where)
doubs.mrt.indval$pval
 
# Extraire les espèces indicatrices à chaque noeud et leur valeur indval
doubs.mrt.indval$maxcls[which(doubs.mrt.indval$pval<=0.05)]
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). Le sommaire indique également quels sites composent chaque branche.

La deuxième partie du code permet de tester la significativité des valeurs de l'indice IndVal pour chaque espèce à l'aide d'un test par permutation. Pour chaque espèce indicatrice significative, la fonction a extrait la branche pour laquelle l'espèce est indicatrice ainsi que la valeur de l'indice IndVal. Dans le cas présent, TRU, VAI et LOC sont toutes des espèces indicatrices de la branche de gauche, mais TRU a la valeur Indval la plus élevée (0.867).

Défi 4: Faire un arbre de régression multivarié pour les données sur les acariens. Sélectionnez l'arbre de taille minimale à 1 écart-type de la valeur minimale CVRE. Quelle est la proportion de variance expliquée par l'arbre? Combien de feuilles contient cet arbre? Quelles en sont les espèces discriminantes?

Défi 4 - Solution

4. Analyse discriminante linéaire

L’analyse discriminante linéaire (ADL ou LDA) est une méthode sous contrainte permettant de déterminer si une matrice de variables indépendantes explique bien un groupement préétabli. Ce groupement peut avoir été obtenu par une méthode de groupement effectuée au préalable (voir Atelier 8) ou par une hypothèse (ex. regroupement de sites/objets selon la latitude ou selon différents traitements expérimentaux). Une LDA peut aussi être utilisée pour classifier de nouvelles données dans ces groupes prédéfinis. On pourrait par exemple être intéressé à prédire l’appartenance d’une espèce de poisson à un groupe selon sa morphologie. On pourrait aussi déterminer si un nouvel article concerne un écosystème terrestre, marin ou d’eau douce selon une classification existante d’articles dans ces biomes effectuée à partir de mots clés de résumés.

La LDA compile des fonctions discriminantes à partir de descripteurs centrés-réduits. Les coefficients obtenus quantifient la contribution relative des variables explicatives sur la discrimination des objets. Les fonctions d’identification peuvent être générées à partir des descripteurs originaux pour classifier de nouvelles données dans les groupes pré-définis.

Nous effectuerons une LDA sur les données de la rivière Doubs. Au préalable, nous devons nous assurer que les matrices de covariances des variables explicatives sont homogènes.

Premièrement, nous voulons effectuer une classification à priori qui est indépendante des variables environnementales. Généralement, les variables environnementales changent avec la latitude (Budyko 1969). Nous classifierons ici les données d’abondance de poissons de la rivière Doubs en fonction de la latitude pour déterminer à quel point les variables environnementales expliquent ces groupements par latitude. Les groupes sont déterminés en divisant l’étendue des latitudes également entre trois groupes/classes et en assignant chaque site à sa classe de latitude.

| Charger les données spatiales et classifier les abondances de poissons par site en fonction de la latitude

#Charger les données spatiales pour déterminer les groupes
spa <- read.csv ('http://www.davidzeleny.net/anadat-r/data-download/DoubsSpa.csv', row.names = 1)
spa <- spa[,-8]
 
#Visualiser la matrice de données 
View (spa)
 
#Ajouter les numéros de site
numbers<-(1:30)
numbers<-numbers[!numbers%in%8] 
spa$site<-numbers
 
#Faire des groupes en fonction de la latitudey<82=group1, 82<y<156=group2, y>156=group3
spa.group<-ddply(.data=spa, .variables=.(x, y, site), .fun= summarise, group = if(y <= 82) 1 else if (y <= 156) 2 else 3)
 
#Ordonner par site
spa.group<-spa.group[with(spa.group, order(site)), ]

Normalement, nous devons vérifier que les matrices de covariance des variables explicatives sont homogènes. Pour les besoins de l'atelier nous n'effectuerons pas cette étape mais vous pourrez consulter Borcard et al. (2011) pour la suite.

Après avoir fait la LDA nous pouvons utiliser les résultats pour déterminer 1- où les sites sont classifiés par rapport aux variables environnementales et 2- quelles sont les probabilités a posteriori que les sites appartiennent aux groupes et 3- le pourcentage de classification correctes basées sur les classes de latitude.

| Faire l'analyse discriminante linéaire (LDA)

#faire la LDA
LDA<-lda(env,spa.group[,4])
 
#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[,4], spe.class)
 
#proportion de classification correcte
diag(prop.table(spe.table, 1))

Les résultats suggèrent que les facteurs environnementaux expliquent le premier et le troisième groupe parfaitemnet, mais le classement des sites dans le deuxième groupe était prédit à 83%. Que peut-on retenir de notre classification? Il y a possiblement des contrastes plus évidents entre les latitudes les plus basses et les plus élevés et que le groupe du milieu est un mélange des deux?

Tentons maintenant de classifier de nouveaux sites en se basant sur les liens que nous avons établi entre nos classes de latitude et les variables environnementales en utilisant la LDA. En utilisant la fonction predict(), nous pouvons charger une nouvelle matrice de sites et les classifier en utilisant les objets de la LDA.

Chargez le fichier classifyme.csv. Celui-ci contient des données fictives de 5 nouveaux sites.

| Charger les données fictives et prédire le groupement des nouvelles données

#prédire la classification des nouvelles données 
#charger les nouvelles données 
classify.me<-read.csv("classifyme.csv", header = T)
 
#prédire le groupement des nouvelles données
predict.group<-predict(LDA, newdata=classify.me)
 
#donner la classification pour chaque site
group.new<-predict.group$class

Nos nouveaux sites, en ordre, ont été classifiés dans les groupes 1, 1, 1, 3 et 3 respectivement.

Défi 5: Faire une LDA sur les données environnementales des acariens (deux premières variables) en se basant sur 4 groupes de latitude à créer à partir des données mite.xy. Quelle proportion de sites ont été classifiés correctement au groupe 1? Au groupe 2?

Défi 5: Solution

5. Autres méthodes d'ordination utiles

| autres méthodes

?cca # Analyse canonique des correspondances (Constrained Correspondence Analysis, CCA)
 
# Méthode d’analyse canonique similaire à la RDA préservant les distance de Chi-carré entre objets (au lieu des 
# distances euclidiennes dans le cas d’une RDA). Cette méthode est mieux adaptée à l’étude de longs gradients
# que la RDA.
 
 
 
?CCorA # Analyse canonique de corrélations (Canonical Correlation Analysis, CCorA)
 
# Cette analyse canonique diffère d’une RDA puisque les deux matrices étudiées sont considérées symétriques
# tandis que dans une RDA la matrice Y dépend toujours de X. Cette méthode est principalement utilisée pour
# tester la significativité des corrélations entre deux jeux de données multivariés et explorer la structure des
# données présentes. 
 
help(coinertia, package=ade4) # Coinertia Analysis 
 
help(coinertia, package=ade4) # Analyse de co-inertie (Coinertia Analysis, CoIA)
 
# Méthode d’analyse canonique symétrique permet de comparer des jeux de données jouant des rôles
# équivalents lors de l’analyse. Cette méthode calcule un espace d’ordination commun dans lequel les objets et les
# variables des deux jeux de données sont projetées et comparées. Par rapport à une CCorA, la CoIA n’impose pas
# de contraintes vis-à-vis du nombre de variables des deux jeux de données, et permet donc de comparer des
# communautés, même si elles sont riches en espèces. En revanche, la CoIA n’est pas adaptée à l’analyse de jeux de 
# données appariés.
 
 
help(mfa, package=ade4) # Analyse factorielle multiple (Multiple Factorial Analysis, MFA)
 
# Méthode permettant de comparer plusieurs jeux de données décrivant les mêmes objets en projetant les objets 
# et variables sur une PCA globale calculée à partir de tous les jeux de données. 
 
 
# Les packages AEM et PCNM permettent d’effectuer diverses formes d’analyses spatiales :
https://r-forge.r-project.org/R/?group_id=195

Références

Alday & Marrs (2014). A simple test for alternative states in ecological restoration: the use of principal response curves. Journal of Vegetation Science, 17, 302-311.

Borcard, Gillet & Legendre (2011). Numerical Ecology with R. Springer New York.

Breiman, L., J. H. Friedman, et al. (1984). Classification and Regression Trees. Belmont, California, USA, Wadsworth International Group.

Budyko, M.I. (1969) The effect of solar radiation variations on the climate of the Earth. Tellus, 21(5), 611-619.

Clarke & Warwick (2001). Change in Marine Communities: An Approach to Statistical Analysis and Interpretation 2nd edition. Primer-E Ltd.

De'ath, G. (2002). Multivariate regression trees : a new technique for modeling species-environment relationships. Ecology, 83(4), 1105–1117.

Gotelli & Ellison (2004). A Primer of Ecological Statistics. Sinaeuer Associates Inc., Sunderland MA.

Legendre & Legendre (2012). Numerical Ecology 3rd edition. Elsevier Science BV, Amsterdam.

Poulin, Andersen & Rochefort (2013) A new approach for tracking vegetation change after restoration: a case study with peatlands. Restoration Ecology, 21, 363-371.