Differences
This shows you the differences between two versions of the page.
| Both sides previous revision Previous revision Next revision | Previous revision | ||
|
r_atelier9 [2016/04/05 18:52] emmanuelle.chretien [4. Analyse discriminante linéaire] |
r_atelier9 [2021/10/14 03:52] (current) lsherin |
||
|---|---|---|---|
| Line 1: | Line 1: | ||
| - | ======= Série d'ateliers en R du CSBQ ======= | + | <WRAP group> |
| + | <WRAP centeralign> | ||
| + | <WRAP important> | ||
| + | <wrap em> __AVIS IMPORTANT__ </wrap> | ||
| - | [[http://qcbs.ca/fr/|{{:logo_text.png?nolink&500|}}]] | + | <wrap em> Depuis l'automne 2021, ce wiki a été discontinué et n'est plus activement développé. </wrap> |
| - | Cette série de [[r|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. | + | <wrap em> Tout le matériel mis à jour et les annonces pour la série d'ateliers R du CSBQ se trouvent maintenant sur le [[https://r.qcbs.ca/fr/workshops/r-workshop-09/|site web de la série d'ateliers R du CSBQ]]. Veuillez mettre à jour vos signets en conséquence afin d'éviter les documents périmés et/ou les liens brisés. </wrap> |
| - | ====== Atelier 9: Analyses multivariées avancées ====== | + | <wrap em> Merci de votre compréhension, </wrap> |
| - | Développé par: Monica Granados, Emmanuelle Chrétien, Bérenger Bourgeois, Amanda Winegardner and Xavier Giroux-Bougard. (Matériel des codes R adaptés de: Borcard, Gillet & Legendre (2011). //Numerical Ecology with R//. Springer New York.) | + | <wrap em> Vos coordonnateurs de la série d’ateliers R du CSBQ. </wrap> |
| + | </WRAP> | ||
| + | </WRAP> | ||
| + | <WRAP clear></WRAP> | ||
| - | **Résumé:** Durant cet atelier, vous apprendrez à faire des analyses multivariées avancées sur vos données. Cet atelier se concentre sur les méthodes sous contraintes telles que l'analyse canonique de redondances, l'arbre de régression multivarié et l'analyse discriminante linéaire afin d'explorer comment les variables environnementales peuvent expliquer la composition en espèces entre différentes sites. | + | ======= Ateliers R du CSBQ ======= |
| + | [[http://qcbs.ca/fr/|{{:logo_text.png?nolink&500|}}]] | ||
| - | Lien vers le Prezi associé: [[https://prezi.com/iuaffmpxhie7/csbq-atelier-r-9/|Prezi]] | + | Cette série de [[r|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. |
| - | Téléchargez le code R, les librairies et données requises pour cet atelier: | + | //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// |
| - | * [[http://qcbs.ca/wiki/_media/code_atelier9.r| Code R]] | + | |
| - | * [[http://qcbs.ca/wiki/_media/DoubsEnv.csv|données DoubsEnv]] | + | |
| - | * [[http://qcbs.ca/wiki/_media/DoubsSpe.csv|données DoubsSpe]] | + | |
| - | * [[http://qcbs.ca/wiki/_media/classifyme.csv|données test pour l'analyse linéaire discriminante]] | + | |
| - | * [[http://qcbs.ca/wiki/_media/mvpart_1.6-2.tar.gz|librairie mvpart]] | + | |
| - | * [[http://qcbs.ca/wiki/_media/MVPARTwrap_0.1-9.tar.gz|librairie MVPARTwrap]] | + | |
| - | * [[http://qcbs.ca/wiki/_media/rdaTest_1.10.tar.gz|librairie rdaTest]] | + | |
| + | <wrap em>AVIS IMPORTANT: MISES À JOUR MAJEURES</wrap> | ||
| - | Assurez-vous d'importer les librairies suivantes dans R Studio (procédure fournie dans le code R): | + | **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 [[https://github.com/QCBSRworkshops/workshop09|la page GitHub]] des ateliers R du CSBQ. |
| - | * [[http://cran.r-project.org/web/packages/vegan/index.html|vegan (pour analyses multivariées)]] | + | |
| - | * [[http://cran.r-project.org/web/packages/vegan/index.html|labdsv (pour identification d'espèces indicatrices pour l'arbre de régression multivarié)]] | + | |
| - | * [[http://cran.r-project.org/web/packages/vegan/index.html|plyr (classification pour l'analyse linéaire discriminante)]] | + | |
| - | * [[http://cran.r-project.org/web/packages/vegan/index.html|MASS (pour l'analyse linéaire discriminante)]] | + | |
| - | * mvpart | + | |
| - | * MVPARTwrap | + | |
| - | * rdatest | + | |
| + | Le matériel disponible inclut; | ||
| + | - La [[https://qcbsrworkshops.github.io/workshop09/pres-fr/workshop09-pres-fr.html|présentation Rmarkdown]] pour cet atelier; | ||
| + | - Un [[https://qcbsrworkshops.github.io/workshop09/book-fr/workshop09-script-fr.R|script R]] qui suit la présentation - //*en construction*//. | ||
| + | - [[https://qcbsrworkshops.github.io/workshop09/book-fr/index.html|Le matériel écrit]] qui accompagne la présentation en format bookdown - //*en construction*//. | ||
| + | ====== Atelier 9: Analyses multivariées ====== | ||
| - | <code rsplus | Installer et importer les librairies> | + | Développé par : Bérenger Bourgeois, Xavier Giroux-Bougard, Amanda Winegardner, Emmanuelle Chrétien et Monica Granados. |
| - | install.packages("vegan") | + | (Le script R est en parti issu de: Borcard, Gillet & Legendre (2011). //Numerical Ecology with R//. Springer New York). |
| - | 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. | + | **Résumé:** Dans cet atelier, vous apprendrez les bases des analyses multivariées qui vous |
| - | # Pour ce faire, cliquez sur l'onglet "Packages" de la section en bas à droite sur R Studio | + | permettront de révéler les patrons de diversité dans vos données de communautés. Vous apprendrez d'abord comment choisir les mesures de distance et les transformations appropriées pour ensuite réaliser plusieurs types d'analyses multivariées: des groupements, des Analyses en Composantes Principales (PCA), des Analyses de Correspondance (CA), des Analyses en Coordonnées Principales (PCoA) et des Positionnements Multidimensionnels Non-Métriques (NMDS). |
| - | # 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) | + | **Lien vers la nouvelle [[https://qcbsrworkshops.github.io/workshop09/workshop09-fr/workshop09-fr.html|présentation Rmarkdown]]** |
| - | library(mvpart) | + | |
| - | library(MVPARTwrap) | + | |
| - | library(rdaTest) | + | |
| - | library(labdsv) | + | |
| - | library(plyr) | + | |
| - | library(MASS) | + | |
| - | </code> | + | |
| + | Lien vers la présentation Prezi associée : [[https://prezi.com/puivxp0qd-tz/?utm_campaign=share&utm_medium=copy | ||
| + | | Prezi]] | ||
| + | **Téléchargez le script R et les données pour cet atelier:** | ||
| + | * [[http://qcbs.ca/wiki/_media/multivar1_f.r| script R]] | ||
| + | * [[http://qcbs.ca/wiki/_media/DoubsEnv.csv|DoubsEnv]] | ||
| + | * [[http://qcbs.ca/wiki/_media/DoubsSpe.csv|DoubsSpe]] | ||
| + | * [[http://qcbs.ca/wiki/_media/coldiss.R|Coldiss (fonction R)]] | ||
| + | **Télechargez les paquets R pour cet atelier:** | ||
| + | * [[http://cran.r-project.org/web/packages/vegan/index.html|vegan]] | ||
| + | * [[http://cran.r-project.org/web/packages/gclus/index.html|gclus]] | ||
| + | * [[http://cran.r-project.org/web/packages/ape/index.html|ape]] | ||
| - | ======Introduction aux analyses multivariées avancées====== | + | <code rsplus | Chargez les paquets et les fonctions nécessaires> |
| + | install.packages("vegan") | ||
| + | install.packages("gclus") | ||
| + | install.packages("ape") | ||
| + | library(vegan) | ||
| + | library(gclus) | ||
| + | library(ape) | ||
| + | source(file.choose()) #coldiss.R | ||
| + | </code> | ||
| + | ======Qu'est-ce que l'ordination?====== | ||
| - | L'atelier précédent donnait un aperçu des analyses multivariées de base: | + | L'ordination est un ensemble de méthodes pour décrire des échantillons dans de multiples dimensions (Clarke et Warwick 2011). Les méthodes d'ordination sont donc très utiles pour simplifier et interpréter des données multivariées. Les écologistes parlent souvent de "faire une PCA" face à des données multidimensionnelles complexes et désordonnées. Programmer des méthodes d'ordination à partir de R est relativement simple. L'interprétation des analyses d'ordination peut par contre être plus difficile, surtout si vous n'êtes pas sûr des questions biologiques que vous souhaitez explorer avec la méthode d'ordination que vous utilisez. Un examen attentif des objectifs de ces méthodes et de leur cadre d'application est nécessaire pour obtenir de bons résultats ! |
| - | * 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. | + | Lorsque vous utilisez une méthode d'ordination, un ensemble de variables est utilisé pour ordonner des échantillons (objets, sites, etc.) le long d'axes principaux représentant des combinaison de variables (Gotelli et Ellison 2004). L'ordination permet donc de réduire ou de simplifier les données en créant de nouveaux axes intégrant la majeure partie de la variation présentes dans les données. A titre d'exemple, un ensemble de données avec 24 variables peut être réduit à cinq composantes principales qui représentent les principaux gradients de variation entre les échantillons. Les méthodes d'ordination sans contrainte ne sont pas adaptées au test d'hypothèses biologiques, mais permettent l'analyse exploratoire des données. Voir ([[http://ordination.okstate.edu/|the Ordination Website]]) pour un aperçu des différents méthodes (en Anglais). |
| + | ======Introduction aux données====== | ||
| - | ======Introduction aux données====== | + | Nous allons utiliser deux principaux ensembles de données dans la première partie de cet atelier. "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 [[http://en.wikipedia.org/wiki/Wide_and_narrow_data|format long (en anglais)]]. |
| - | 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 [[http://en.wikipedia.org/wiki/Wide_and_narrow_data|format long (en anglais)]]. | + | <code rsplus | Chargez DoubsSpe et DoubsEnv> |
| - | + | ||
| - | <code rsplus | Chargez les données> | + | |
| #Matrice d'abondances d'espèces: “DoubsSpe.csv” | #Matrice d'abondances d'espèces: “DoubsSpe.csv” | ||
| spe<- read.csv(file.choose(), row.names=1) | 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 | + | spe<- spe[-8,] #Pas d'espèces dans le site 8, supprimer site 8. |
| + | #Exécutez cette ligne une seule fois. | ||
| - | #Matrice de données environnementales: “DoubsEnv.csv” | + | #L'environnement: “DoubsEnv.csv” |
| env<- read.csv(file.choose(), row.names=1) | 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. | + | env<- env[-8,] #Supprimer site 8 puisqu'on l'a retiré de la matrice d'abondances. |
| + | #Exécutez cette ligne une seule fois. | ||
| </code> | </code> | ||
| + | ======1. Exploration des données====== | ||
| + | =====1.1 Données sur les espèces===== | ||
| - | ======1. Exploration et préparation des données====== | + | 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). |
| - | + | ||
| - | =====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). | + | |
| <code rsplus | Explorez DoubsSpe> | <code rsplus | Explorez DoubsSpe> | ||
| - | names(spe) #voir les noms des colonnes spe | + | names(spe) # Les noms des colonnes |
| - | dim(spe) #dimensions de spe; nombre de lignes et de colonnes | + | dim(spe) # Le nombre de lignes et de colonnes. |
| - | str(spe) #la structure interne de la matrice | + | str(spe) # La structure interne de la matrice. |
| - | head(spe) #les premières lignes | + | head(spe) # Les premières lignes. |
| - | summary(spe) #statistiques descriptives | + | summary(spe) # Les statistiques descriptives. |
| </code> | </code> | ||
| Line 105: | Line 105: | ||
| <code rsplus | La distribution des espèces (DoubsSpe)> | <code rsplus | 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 | + | (ab<-table(unlist(spe))) #Les parenthèses signifient que la sortie s'affiche immédiatement |
| - | barplot(ab, las=1, xlab="Abundance class", ylab="Frequency", col=grey(5:0/5)) | + | barplot(ab, las=1, xlab=”Abundance class”, ylab=”Frequency”, col=grey(5:0/5)) |
| </code> | </code> | ||
| {{:spe_barplot.png?300|}} | {{:spe_barplot.png?300|}} | ||
| - | Il y a une grande fréquence de zéros dans les données d'abondance. | + | Pouvez-vous voir qu'il y a une grande fréquence de zéros dans les données d'abondance ? |
| Calculez le nombre d'absences. | Calculez le nombre d'absences. | ||
| - | <code rsplus | Absences dans les données de poissons> | + | <code rsplus | Absences> |
| sum(spe==0) | sum(spe==0) | ||
| </code> | </code> | ||
| - | Regardez la proportion de zéros dans les données d'abondances de poissons. | + | Regardez la proportion de zéros dans les données de la communauté de poissons. |
| <code rsplus | Proportion de zéros> | <code rsplus | Proportion de zéros> | ||
| sum(spe==0)/(nrow(spe)*ncol(spe)) | sum(spe==0)/(nrow(spe)*ncol(spe)) | ||
| </code> | </code> | ||
| + | La proportion de zéros dans la matrice est de ~0.5 | ||
| - | 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(). | + | Calculez le nombre de sites où chaque espèce est présente. |
| + | <code rsplus | Nombre de sites où chaque espèce est présente> | ||
| + | spe.pres<- colSums(spe>0) # Somme des sites où chaque espèce est présente. | ||
| + | hist(spe.pres, main=”Cooccurrence des espèces”, las=1, xlab=”Fréquence”, breaks=seq(0,30, by=5), col=”grey”) | ||
| + | </code> | ||
| + | Le plus grand nombre d'espèces se retrouvent dans un nombre intermédiaire de sites. | ||
| - | 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). | + | Calculez le nombre d'espèces présentes à chaque site. Ici, nous utilisons une façon simpliste de calculer la richesse en espèces. Dans certains cas, le nombre d'espèces présentes dans un site varie entre les sites car il peut y avoir une relation entre le nombre d'individus comptés (l'abondance) dans ce site et la richesse en espèces. [[http://cc.oulu.fi/~jarioksa/softhelp/vegan/html/diversity.html|richesse d'espèces raréfiée (en anglais)]] est souvent une mesure plus appropriée que la richesse totale. La fonction rarefy () de vegan peut être utilisée pour calculer la richesse raréfiée des espèces. |
| - | + | ||
| - | <code rsplus | transformation de Hellinger> | + | <code rsplus | Richesse en espèces> |
| - | spe.hel <- decostand(spe, method="hellinger") # vous pouvez également utiliser method="hell" | + | site.pres<- rowSums(spe>0) # Nombre d'espèces présentes dans chaque site |
| + | hist(site.pres, main=”Richesse en espèces”, las=1, xlab=”Fréquence des sites”, ylab=”Nombre d'espèces”, breaks=seq(0,30, by=5), col=”grey”) | ||
| </code> | </code> | ||
| + | =====1.2 Données sur l'environnement===== | ||
| - | =====1.2 Données environnementales===== | + | Explorez les données environnementales pour détecter les colinéarités : |
| - | Explorez les données environnemntales et détectez les colinéarités. | + | |
| - | <code rsplus | Exploration de la matrice Doubs Env> | + | <code rsplus | Exploration de DoubsEnv> |
| names(env) | names(env) | ||
| dim(env) | dim(env) | ||
| Line 142: | Line 148: | ||
| head(env) | head(env) | ||
| summary(env) | summary(env) | ||
| - | pairs(env, main="Bivariate Plots of the Environmental Data" ) | + | pairs(env, main="Données environnementales" ) |
| </code> | </code> | ||
| - | 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. | + | 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. |
| - | <code rsplus | Decostand> | + | <code rsplus | Standardisation des données> |
| - | env.z <- decostand(env, method="standardize") | + | env.z<-decostand(env, method="standardize") |
| - | apply(env.z, 2, mean) # les données sont centrées (moyenne~0) | + | apply(env.z, 2, mean) # Les données sont maintenant centrées (moyennes~0)... |
| - | apply(env.z, 2, sd) # les données sont réduites (écart type =1) | + | apply(env.z, 2, sd) # et réduites (écart-type=1) |
| </code> | </code> | ||
| + | ======2. Mesures d'association====== | ||
| + | L'algèbre matricielle est à la base des méthodes d'ordination. Une matrice est constituée de données (ex. valeurs mesurées) réparties en lignes et colonnes. Les analyses d'ordination sont effectuées sur des matrices d'association calculées à partir des matrices de données écologiques (telles que DoubsEnv ou DoubsSpe). La création d'une matrice d'association permet le calcul de similarité et de distance entre les objets ou les descripteurs (Legendre et Legendre 2012). Avant de se lancer dans des analyses d'ordination, il est important de passer du temps sur vos matrices de données. Explorer les mesures possibles d'association qui peuvent être générées à partir de vos données avant de faire une ordination peut vous aider à mieux comprendre quelles mesures de distance sont appropriées pour les méthodes d'ordination. Il peut être difficile de voir l'objectif de chaque indice de dissimilarité, mais cette connaissance sera nécessaire pour mieux comprendre les méthodes d'ordination canoniques présentées par la suite. | ||
| - | ======2. Analyses canoniques====== | + | EN RÉSUMÉ: Pour l'ordination d'objets, il faut calculer les distances entre eux. Ces distances peuvent être calculées de plusieurs façons, en prenant en compte l'abondance ou les données de présence / absence. Plus important encore, plusieurs propriétés sont de grande importance pour les mesures de distances, et elles seront explorées dans les exemples ci-dessous. Pour plus d'informations sur les propriétés des mesures de distance et certains termes clés, cliquez sur la section cachée. |
| - | 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). | + | <hidden> |
| - | =====2.1 Analyse canonique de redondances (ACR ou RDA)===== | + | **Termes clés:** |
| - | 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. | + | **Association -** «Terme général pour décrire toute mesure ou coefficient servant à quantifier la ressemblance ou différence entre les objets ou les descripteurs. Dans une analyse entre descripteurs, zéro signifie l'absence d'association.» (Legendre et Legendre 2012). |
| - | 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 : | + | **Similarité -** une mesure dont «le maximum (S = 1) est atteint lorsque deux objets sont identiques et le minimum lorsque deux objets sont complètement différents.» (Legendre et Legendre 2012). |
| - | + | ||
| - | Yfit = X[X'X]-1 [X'Y] | + | |
| + | **Distance (aussi «dissimilarité») -** une mesure dont «le maximum (D=1) est atteint lorsque deux objects sont complètement différent.» (Legendre et Legendre 2012). Distance ou dissimilarité (D) = 1-S | ||
| - | 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. | + | Le choix d'une mesure d'association dépend de vos données, mais aussi de ce que vous en savez d'un point de vue écologique. Par exemple, la distance euclidienne est une mesure de distance très commune, facile à utiliser et utile pour comprendre comment les différences entre deux échantillons sont basées sur la cooccurrence des espèces. Le calcul de la distance euclidienne prend en compte les zéros dans les données, ce qui signifie que deux échantillons ou sites ne contenant aucune espèce en commun (double-absence) peuvent sembler plus similaires que deux sites partageant quelques espèces. Dans ce cas, la distance euclidienne peut être trompeuse et il est souvent préférable de choisir une mesure de distance différente si beaucoup d'espèces ont une abondance nulle dans votre matrice. Cette propriété est communément appelée le problème des «double zéros» en ordination. |
| - | {{ :constrained_ord_diagram.png |}} | + | Quelques mesures de distance (d'après Gotelli et Ellison 2004): |
| + | ^Mesure^Propriété^Description^ | ||
| + | ^Euclidienne^Métrique^Distance entre deux points dans un espace en 2D.^ | ||
| + | ^Manhattan^Métrique^Distance entre deux points - la distance est la somme des différences entre coordonnées cartésiennes.^ | ||
| + | ^Corde^Métrique^Généralement utilisée pour déterminer les différences dues à la dérive génétique.^ | ||
| + | ^Mahalanobis^Métrique^Distance entre un point et une distribution, où la distance est le nombre d'écart-types du point correspondant à la moyenne de la distribution.^ | ||
| + | ^Chi-carré^Métrique^Similaire à la distance euclidienne.^ | ||
| + | ^Bray-Curtis^Semi-métrique^Dissimilarité entre deux échantillons (ou sites) où la somme des valeurs minimales des espèces présentes dans les deux sites sont divisées par la somme des espèces répertoriées dans chaque site.^ | ||
| + | ^Jaccard^Métrique^[[http://en.wikipedia.org/wiki/Jaccard_index|Description]] ^ | ||
| + | ^Sorensen's^Semi-métrique^Bray-Curtis correspond à 1 - Sorensen^ | ||
| + | </hidden> | ||
| - | Une RDA est un processus en deux étapes (d’après Legendre et Legendre, 2012) | + | =====2.1 Mesures de distance===== |
| - | Plusieurs valeurs statistiques peuvent être extraites de la RDA, en particulier : | + | **Les données quantitatives des espèces** |
| + | Nous pouvons utiliser la fonction vegdist () pour calculer des indices de dissimilarité sur des données de composition de la communauté. Ceux-ci peuvent ensuite être visualisés sous forme de matrice si désiré. | ||
| - | - 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, | + | <code rsplus | vegdist> |
| - | - 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, | + | spe.db<-vegdist(spe, method="bray") # distance de Bray (avec des données de présence-absence, correspond à Sorensen) |
| - | - La statistique F correspond à un test global de la significativité de la RDA en comparant le modèle étudié à un modèle nul. | + | spe.dj<-vegdist(spe, method="jac") # distance de Jaccard |
| + | spe.dg<-vegdist(spe, method="gower") # distance de Gower | ||
| + | spe.db<-as.matrix(spe.db) # réarranger en format matrice (pour visualisation, ou pour exporter en .csv) | ||
| + | </code> | ||
| - | 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. | + | Une version condensée de la matrice spe.db représentant la distance entre les trois premières espèces de DoubsSpe ressemblerait à ceci: |
| + | ^ ^Espèce 1^Espèce 2^Espèce 3^ | ||
| + | ^Espèce 1^0.0^0.6^0.68^ | ||
| + | ^Espèce 2^0.6^0.0^0.14^ | ||
| + | ^Espèce 3^0.68^0.14^0.0^ | ||
| + | Vous pouvez voir que lorsque l'on compare une espèce à elle-même (par exemple, Espèce 1 à Espèce 1), la distance est de 0 puisque les espèces sont identiques. | ||
| - | Avec R, la RDA peut être calculée en utilisant la fonction rda du package vegan, comme suit: | + | Ces mêmes mesures de distance peuvent être calculées à partir des données de présence-absence par l'utilisation de l'argument binary=TRUE dans la fonction vegdist(). Cela donnera des mesures de distance légèrement différentes. |
| + | Vous pouvez également créer des représentations graphiques de ces matrices d'association en utilisant la fonction coldiss. | ||
| + | <hidden> | ||
| + | Cette fonction peut être sourcée à partir du script coldiss.R: | ||
| + | <code rsplus | coldiss> | ||
| + | windows() | ||
| + | coldiss(spe.db, byrank=FALSE, diag=TRUE) # Carte des points chauds Bray-Curtis | ||
| + | windows() | ||
| + | coldiss(spe.dj, byrank=FALSE, diag=TRUE) # Carte des points chauds Jaccard | ||
| + | windows() | ||
| + | coldiss(spe.dg, byrank=FALSE, diag=TRUE) # Carte des points chauds Gower | ||
| + | </code> | ||
| + | {{:coldiss_Bray.png?800|}} | ||
| - | <code rsplus | effectuer une RDA avec rda() de vegan> | + | La figure montre une matrice de dissimilarité dont les couleurs reflètent la mesure de distance. La couleur violet est associée aux zones de fortes dissimilarités. |
| - | #Préparer les données | + | </hidden> |
| - | 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 | + | **Données environnementales quantitatives** |
| - | summary(spe.rda, display=NULL) | + | Regardons //les associations// entre les variables environnementales (aussi appelée mode Q): |
| - | + | <code rsplus | distances mesurées à partir des données environnementales> | |
| - | #Les résultats peuvent être vus en utilisant summary(): | + | env.de<-dist(env.z, method = "euclidean") # matrice de distances euclidiennes des données env. standardisées |
| - | summary(spe.rda, display=NULL) #display = NULL optional | + | windows() # crée une nouvelle fenêtre graphique |
| + | coldiss(env.de, diag=TRUE) | ||
| </code> | </code> | ||
| + | Nous pouvons ensuite regarder //la dépendance// entre les variables environnementales (aussi appelée mode R): | ||
| + | <code rsplus | corrélations entre les données environnementales> | ||
| + | (env.pearson<-cor(env)) # coefficient r de corrélation de Pearson | ||
| + | round(env.pearson, 2) # arrondit les coefficients à deux décimales | ||
| + | (env.ken<-cor(env, method="kendall")) # coefficient tau de corrélation de rang de Kendall | ||
| + | round(env.ken, 2) | ||
| + | </code> | ||
| - | Le sommaire ressemble à ceci: | + | La corrélation de Pearson mesure la corrélation linéaire entre deux variables. La corrélation de Kendall est une corrélation de rang qui quantifie la relation entre deux descripteurs ou deux variables lorsque les données sont ordonnées au sein de chaque variable. |
| + | Dans certains cas, il peut y avoir des types mixtes de variables environnementales. Le mode Q peut alors être utilisé pour trouver des associations entre variables environnementales. C'est ce que nous allons faire avec l'exemple fictif suivant: | ||
| + | <code rsplus | Exemple> | ||
| + | var.g1<-rnorm(30, 0, 1) | ||
| + | var.g2<-runif(30, 0, 5) | ||
| + | var.g3<-gl(3, 10) | ||
| + | var.g4<-gl(2, 5, 30) | ||
| + | (dat2<-data.frame(var.g1, var.g2, var.g3, var.g4)) | ||
| + | str(dat2) | ||
| + | summary(dat2) | ||
| + | </code> | ||
| - | {{ :rda_output_1.png |}} | + | Une matrice de dissimilarité peut être générée pour ces variables mixtes en utilisant la distance de Gower: |
| + | <code rsplus | daisy> | ||
| + | ?daisy #Cette fonction peut gérer la présence de NA dans les données | ||
| + | (dat2.dg<-daisy(dat2, metric="gower")) | ||
| + | coldiss(dat2.dg) | ||
| + | </code> | ||
| - | 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. | + | **Défi 1 - Niveau intermédiaire** |
| + | Discutez avec votre voisin: Comment pouvons-nous dire si des objets sont similaires avec un jeu de données multivariées? Faites une liste de toutes vos suggestions. | ||
| - | 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): | + | **Défi 1 - Solution** |
| + | <hidden> | ||
| + | Discussion avec le groupe. | ||
| + | </hidden> | ||
| + | **Défi 1 - Niveau avancé** | ||
| + | Calculer à la mitaine // sans utiliser la fonction decostand () // les distances de Bray-Curtis et de Gower pour l'abondance des espèces CHA, TRU et VAI dans les sites 1, 2 et 3. | ||
| - | <code rsplus | ordiR2step() pour la sélection progressive> | + | **Défi 1 - Solution** |
| - | ?ordiR2step | + | <hidden> |
| - | 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")) | + | |
| - | </code> | + | |
| + | Formule pour calculer la distance de Bray-Curtis: d[jk] = (sum abs(x[ij]-x[ik]))/(sum (x[ij]+x[ik])) | ||
| - | 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 : | + | Réduisez le jeu de données aux espèces CHA, TRU et VAI et aux sites 1, 2 et 3 |
| + | <code rsplus | Subset> | ||
| + | spe.challenge<-spe[1:3,1:3] # les 3 premières lignes et 3 premières espèces (colonnes) | ||
| + | </code> | ||
| - | <code rsplus | RDA avec seulement les variables significatives> | + | Déterminer l'abondance totale des espèces pour chaque site d'intérêt (somme des trois lignes) qui correspondra au dénominateur de la distance de Bray-Curtis. |
| - | spe.rda.signif <- rda(spe.hel~., data=env.signif) | + | <code rsplus | Abondance par site> |
| - | summary(spe.rda.signif, display=NULL) | + | >(Abund.s1<-sum(spe.challenge[1,])) |
| + | (Abund.s2<-sum(spe.challenge[2,])) | ||
| + | (Abund.s3<-sum(spe.challenge[3,])) | ||
| </code> | </code> | ||
| + | Maintenant, calculez la différence de l'abondance des espèces pour chaque paire de sites. Par exemple, quelle est la différence entre l'abondance de CHA et TRU dans le site 1? Vous devez calculer les différences suivantes: | ||
| + | CHA et TRU site 1 | ||
| + | CHA et VAI site 1 | ||
| + | TRU et VAI site 1 | ||
| + | CHA et TRU site 2 | ||
| + | CHA et VAI site 2 | ||
| + | TRU et VAI site 2 | ||
| + | CHA et TRU site 3 | ||
| + | CHA et VAI site 3 | ||
| + | TRU et VAI site 3 | ||
| - | Les variables explicatives expliquent désormais 59% de la variance de Y. | + | <code rsplus | Différences d'abondance> |
| + | Spec.s1s2<-0 | ||
| + | Spec.s1s3<-0 | ||
| + | Spec.s2s3<-0 | ||
| + | for (i in 1:3) { | ||
| + | Spec.s1s2<-Spec.s1s2+abs(sum(spe.challenge[1,i]-spe.challenge[2,i])) | ||
| + | Spec.s1s3<-Spec.s1s3+abs(sum(spe.challenge[1,i]-spe.challenge[3,i])) | ||
| + | Spec.s2s3<-Spec.s2s3+abs(sum(spe.challenge[2,i]-spe.challenge[3,i])) } | ||
| + | </code> | ||
| - | Le R² ajusté de cette RDA est calculé à partir de la fonction RsquareAdj : | + | Maintenant, utilisez les différences calculées comme numérateur et l'abondance totale de l'espèce comme dénominateur pour retrouver l'équation de la distance de Bray-Curtis. |
| + | <code rsplus | Distance de Bray-Curtis> | ||
| + | (db.s1s2<-Spec.s1s2/(Abund.s1+Abund.s2)) #1 comparé à 2 | ||
| + | (db.s1s3<-Spec.s1s3/(Abund.s1+Abund.s3)) #1 comparé à 3 | ||
| + | (db.s2s3<-Spec.s2s3/(Abund.s2+Abund.s3)) #2 comparé à 3 | ||
| + | </code> | ||
| - | + | Vérifiez vos résultats en utilisant la fonction vegdist () : | |
| - | <code rsplus | R2 ajusté> | + | <code rsplus | Vérification des résultats> |
| - | (R2adj <- RsquareAdj(spe.rda.signif)$adj.r.squared) | + | (spe.db.challenge<-vegdist(spe.challenge, method="bray")) |
| </code> | </code> | ||
| + | Une matrice comme celle-ci est calculée et devrait être correspondre à vos calculs manuels: | ||
| + | ^ ^Site 1^Site 2^ | ||
| + | ^Site 2^0.5^--^ | ||
| + | ^Site 3^0.538^0.0526^ | ||
| - | 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): | + | Pour la distance de Gower, procédez de la même façon, mais utiliser l'équation appropriée: |
| + | Distance de Gower: d[jk] = (1/M) sum(abs(x[ij]-x[ik])/(max(x[i])-min(x[i]))) | ||
| + | <code rsplus | Distance de Gower> | ||
| + | # Calculer le nombre de colonnes | ||
| + | M<-ncol(spe.challenge) | ||
| + | # Calculer les différences d'abondance de chaque espèce entre paires de sites | ||
| + | Spe1.s1s2<-abs(spe.challenge[1,1]-spe.challenge[2,1]) | ||
| + | Spe2.s1s2<-abs(spe.challenge[1,2]-spe.challenge[2,2]) | ||
| + | Spe3.s1s2<-abs(spe.challenge[1,3]-spe.challenge[2,3]) | ||
| + | Spe1.s1s3<-abs(spe.challenge[1,1]-spe.challenge[3,1]) | ||
| + | Spe2.s1s3<-abs(spe.challenge[1,2]-spe.challenge[3,2]) | ||
| + | Spe3.s1s3<-abs(spe.challenge[1,3]-spe.challenge[3,3]) | ||
| + | Spe1.s2s3<-abs(spe.challenge[2,1]-spe.challenge[3,1]) | ||
| + | Spe2.s2s3<-abs(spe.challenge[2,2]-spe.challenge[3,2]) | ||
| + | Spe3.s2s3<-abs(spe.challenge[2,3]-spe.challenge[3,3]) | ||
| - | <code rsplus | anova.cca pour tester la significativité des axes> | + | # Calculer l'étendue d'abondance de chaque espèces parmi les sites |
| - | ?anova.cca | + | Range.spe1<-max(spe.challenge[,1]) - min (spe.challenge[,1]) |
| - | anova.cca(spe.rda.signif, step=1000) | + | Range.spe2<-max(spe.challenge[,2]) - min (spe.challenge[,2]) |
| - | anova.cca(spe.rda.signif, step=1000, by="axis") | + | Range.spe3<-max(spe.challenge[,3]) - min (spe.challenge[,3]) |
| - | #In this case, the RDA model is highly significant (p=0.001) as well as the two significant axes. | + | |
| + | # Calculer la distance de Gower | ||
| + | (dg.s1s2<-(1/M)*((Spe2.s1s2/Range.spe2)+(Spe3.s1s2/Range.spe3))) | ||
| + | (dg.s1s3<-(1/M)*((Spe2.s1s3/Range.spe2)+(Spe3.s1s3/Range.spe3))) | ||
| + | (dg.s2s3<-(1/M)*((Spe2.s2s3/Range.spe2)+(Spe3.s2s3/Range.spe3))) | ||
| + | |||
| + | # Vérifier vos résultats | ||
| + | (spe.db.challenge<-vegdist(spe.challenge, method="gower")) | ||
| </code> | </code> | ||
| + | </hidden> | ||
| - | 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é: | + | =====2.2 Transformations des données de composition des communautés===== |
| + | Les données de composition des communautés peuvent également être standardisées ou transformées. La fonction decostand () de vegan fournit des options de standardisation et de transformation de ce type de données | ||
| - | <code rsplus | représentations graphiques des RDAs> | + | Transformer les abondances en données de présence-absence: |
| - | #Pour faire des graphiques rapidement et facilement | + | <code rsplus | decostand, présence-absence> |
| - | windows() | + | spe.pa<-decostand(spe, method="pa") |
| - | 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) | + | |
| </code> | </code> | ||
| + | D'autres transformations peuvent être utilisées pour corriger l'influence d'espèces rares, par exemple, la transformation de Hellinger: | ||
| + | <code rsplus | Hellinger et Chi-carré> | ||
| + | #La transformation Hellinger | ||
| + | spe.hel<-decostand(spe, method="hellinger") # vous pouvez aussi simplement écrire "hel" | ||
| - | Les triplots finaux: | + | #Transformation de chi-carré |
| + | spe.chi<-decostand(spe, method="chi.square") | ||
| + | </code> | ||
| - | {{ :doubs_rda1.png |}} | + | **Défi option 2 - Niveau avancé** |
| - | {{ :doubs_rda2.png |}} | + | Calculez les distances de Hellinger et de Chi-carré sur les données "spe" sans utiliser decostand (). |
| - | **Défi 1**: Effectuer la RDA de l’abondance des espèces du tableau d'acariens en fonction des variables environnementales mite.env | + | **Défi option 2 - Solution** |
| + | <hidden> | ||
| + | La transformation de Hellinger est une transformation qui diminue l'importance accordée aux espèces rares. | ||
| + | <code rsplus | Solution> | ||
| + | # Hellinger | ||
| + | # Calculer l'abondance des espèces par site | ||
| + | (site.totals=apply(spe, 1, sum)) | ||
| - | <code rsplus | Chargez les données d'acariens> | + | # Réduire les abondances d'espèces en les divisant par les totaux par sites |
| - | #Ces données sont disponibles dans la librairie vegan | + | (scale.spe<-spe/site.totals) |
| - | data(mite) | + | |
| - | mite.spe<-mite | + | |
| - | mite.spe.hel <- decostand(mite.spe, method="hellinger") | + | |
| - | data(mite.env) | + | # Calculer la racine carrée des abondances d'espèces réduites |
| - | </code> | + | (sqrt.scale.spe<-sqrt(scale.spe)) |
| + | # Comparer les résultats | ||
| + | sqrt.scale.spe | ||
| + | spe.hel | ||
| + | sqrt.scale.spe-spe.hel # ou: sqrt.scale.spe/spe.hel | ||
| - | 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 ? | + | # Chi-carré |
| + | # Premièrement calculer le total des abondances d'espèces par site | ||
| + | (site.totals<-apply(spe, 1, sum)) | ||
| + | # Ensuite calculer la racine carrée du total des abondances d'espèces | ||
| + | (sqrt.spe.totals<-sqrt(apply(spe, 2, sum))) | ||
| - | **Défi 1**: Solution | + | # Réduire les abondances d'espèces en les divisant par les totaux par sites et les totaux par espèces |
| + | scale.spe2<-spe | ||
| + | for (i in 1:nrow(spe)) { | ||
| + | for (j in 1:ncol(spe)) { | ||
| + | (scale.spe2[i,j]=scale.spe2[i,j]/(site.totals[i]*sqrt.spe.totals[j])) }} | ||
| - | <hidden> | + | #Ajuster les abondances en les multipliant par la racine carrée du total de la matrice des espèces |
| - | Votre code devrait ressembler à ceci: | + | (adjust.scale.spe2<-scale.spe2*sqrt(sum(rowSums(spe)))) |
| - | <code rsplus | RDA sur les données des acariens> | + | #Vérifier les résultats |
| - | #RDA avec toutes les variables environnementales | + | adjust.scale.spe2 |
| - | mite.spe.rda<-rda(mite.spe.hel~., data=mite.env) | + | spe.chi |
| - | + | adjust.scale.spe2-spe.chi # or: adjust.scale.spe2/spe.chi | |
| - | #Sélection des variables environnementales significatives | + | |
| - | ordiR2step(rda(mite.spe.hel~1, data=mite.env), | + | |
| - | scope= formula(mite.spe.rda), direction= "forward", R2scope=TRUE, pstep=1000) | + | |
| - | + | ||
| - | #Créez un nouveau tableau de données avec seulement les variables significatives | + | |
| - | mite.env.signif <- subset(mite.env, | + | |
| - | select = c("WatrCont", "Shrub", "Substrate", "Topo", "SubsDens")) | + | |
| - | + | ||
| - | #Refaire la RDA avec seulement les variables significatives et regardez le sommaire des résultats | + | |
| - | mite.spe.rda.signif=rda(mite.spe~., data=mite.env.signif) | + | |
| - | summary(mite.spe.rda.signif, display=NULL) | + | |
| - | + | ||
| - | #Calculez le R2 ajusté | + | |
| - | (R2adj <- RsquareAdj(mite.spe.rda.signif)$adj.r.squared) | + | |
| - | + | ||
| - | #Déterminez les axes significatifs de la RDA | + | |
| - | anova.cca(mite.spe.rda.signif, step=1000) | + | |
| - | anova.cca(mite.spe.rda.signif, step=1000, by="axis") | + | |
| - | + | ||
| - | #Représentation graphique de la RDA | + | |
| - | windows() | + | |
| - | plot(mite.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(mite.spe.rda.signif, display="sites", choices=c(1,2), scaling=1), | + | |
| - | pch=21, col="black", bg="steelblue", cex=1.2) | + | |
| - | text(scores(mite.spe.rda.signif, display="species", choices=c(1), scaling=1), | + | |
| - | scores(mite.spe.rda.signif, display="species", choices=c(2), scaling=1), | + | |
| - | labels=rownames(scores(mite.spe.rda.signif, display="species", scaling=1)), | + | |
| - | col="grey", cex=0.8) | + | |
| - | arrows(0,0, | + | |
| - | scores(mite.spe.rda.signif, display="bp", choices=c(1), scaling=1), | + | |
| - | scores(mite.spe.rda.signif, display="bp", choices=c(2), scaling=1), | + | |
| - | col="red") | + | |
| - | text(scores(mite.spe.rda.signif, display="bp", choices=c(1), scaling=1)+0.05, | + | |
| - | scores(mite.spe.rda.signif, display="bp", choices=c(2), scaling=1)+0.05, | + | |
| - | labels=rownames(scores(mite.spe.rda.signif, display="bp", choices=c(2), scaling=1)), | + | |
| - | col="red", cex=1) | + | |
| </code> | </code> | ||
| + | </hidden> | ||
| - | Cinq variables explicatives sont significatives : WatrCont, Shrub, Substrate, Topo, and SubsDens. La proportion de variance expliquée par ces variables est de 32.08% et le R2 ajusté de cette RDA s’élève à 19.20%. Bien que le premier axe d’ordination soit significatif (p<0.001), le modèle global de RDA n'est pas significatif (p=0.107). Trois groupes de sites apparaissent sur le triplot. Le premier groupe de sites posséde un contenu élevé en eau (coin supérieur droit) et a de fortes abondances en espèces 9 et 25. Les espèces 6, 12, 19 et 26 sont liés à un second groupe de sites caractérisés par une topographie de type hummock (coin inférieur gauche). Le dernier groupe de sites (coins supérieur gauche) ne semble pas abriter d’espèces particulières et montrent des conditions environnementales hétérogènes. | ||
| - | {{ :mite_rda1.png |}} | ||
| - | </hidden> | + | =====2.3 Groupement===== |
| + | Les matrices d’association nécessaires pour utiliser les méthodes de groupement. Le groupement n’est pas une méthode statistique en tant que telle puisqu’elle ne teste pas d’hypothèse, mais permet de déceler des structures dans les données en partitionnant soit les objets, soit les descripteurs. Les objets similaires sont agrégés en sous-groupes ce qui permet de mettre en relief des cassures (contrastes) entre les données. En tant que biologiste, il peut être intéressant de tenter de séparer une série de sites en groupes en fonction de leurs caractéristiques environnementales ou de leur composition en espèces. | ||
| + | Les résultats d’un groupement sont généralement représentés sous forme de dendrogramme (arbre), dans lequel les objets sont agrégés en groupes. Il y a plusieurs familles de méthodes de groupements, mais nous présenterons uniquement un aperçu de trois méthodes : le groupement agglomératif hiérarchique à liens simples, le groupement agglomératif hiérarchique à liens complets et la méthode de Ward. Pour plus de détails sur les différentes familles de méthodes de groupement, consulter Legendre et Legendre 2012 (chapitre 8). | ||
| - | =====2.2 RDA partielle===== | + | Dans les méthodes hiérarchiques, les éléments des petits ensembles se regroupent en groupes plus vastes de rang supérieur, et ainsi de suite (par exemple : espèces, genres, familles, ordre). Avant de faire le groupement, il faut créer une matrice d’association entre les objets. Une matrice de distances est le choix par défaut des fonctions de groupement dans R. La matrice d’association est premièrement classée en ordre croissant de distances. Ensuite, les groupes sont formés de manière hiérarchique selon les critères spécifiques à chaque méthode. |
| - | 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. | + | Prenons un exemple tout simple d'une matrice de distances euclidiennes entre 5 objets dont on a ordonné les distances en ordre croissant. |
| - | 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. | + | {{ :groupement.exemple.png |}} |
| - | Avec R, la RDA partielle est effectuée en utilisant le fonction rda et en ajoutant un argument correspondant aux co-variables : | + | Pour le groupement agglomératif à liens simples, les deux objets les plus proches se regroupent en premier. Ensuite, un deuxième groupe est formé à partir des deux objets les plus proches suivant (il se peut que ce soit deux objets différents, ou bien un objet et le groupe formé précédemment), et ainsi de suite. Cette méthode forme généralement de longues chaînes de groupes (dans l'exemple ci-haut, les objets 1 à 5 se regroupent successivement). À l’inverse, pour le groupement agglomératif à liens complet, un objet se regroupe à un autre objet/groupe seulement lorsqu’il est aussi lié à l’élément le plus éloigné de ce groupe. Ainsi, quand deux groupes fusionnent, tous les éléments des deux groupes sont liés à la distance considérée (ci-haut, le groupe 3-4 ne se lie au groupe 1-2 qu'à la distance à laquelle tous les autres éléments sont déjà liés). C’est pour cette raison que le groupement à liens complets forme généralement plusieurs petits groupes séparés et qu’elle peut être plus appropriée pour relever des contrastes ou des discontinuités dans les données. |
| + | Comparons ces deux méthodes en utilisant les données d’abondances de poissons de la rivière Doubs. | ||
| + | Les données d’abondances ont été au préalable transformées par la méthode Hellinger. Puisque les fonctions de groupement requièrent une matrice de distances, la première étape sera de générer une matrice de distances Hellinger. | ||
| - | <code rsplus | RDA partielle avec rda()> | + | <code rsplus | vegdist, distance de Hellinger> |
| - | #Divisez le tableau de données environnementales en deux: | + | spe.dhel<-vegdist(spe.hel,method="euclidean") #crée une matrice de distances Hellinger à partir des données d’abondance transformées |
| - | envtopo <- env[, c(1:3)] # Physiographie : tableau 1 | + | |
| - | names(envtopo) | + | #Pour voir la différence entre les deux types d’objets |
| - | envchem <- env[, c(4:10)] # Physico-chimie de l'eau : tableau 2 | + | head(spe.hel)# données d’abondances transformées Hellingerhead(spe.dhel)# matrice de distances de Hellinger entre les sites |
| - | 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) | + | |
| - | 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="grey", 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) | + | |
| - | 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="grey", 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) | + | |
| </code> | </code> | ||
| + | La plupart des méthodes de groupement sont disponible dans la fonction hclust() de la librairie stats | ||
| - | 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%. | + | <code rsplus | hclust, comparaison du groupement à liens simples et à liens complets> |
| + | #Faire le groupement à liens simples | ||
| + | spe.dhel.single<-hclust(spe.dhel, method="single") | ||
| + | plot(spe.dhel.single) | ||
| + | #Faire le groupement à liens complet | ||
| + | spe.dhel.complete<-hclust(spe.dhel, method="complete") | ||
| + | plot(spe.dhel.complete) | ||
| + | </code> | ||
| - | {{ :partial_rda.png |}} | + | {{ :clust_single.png |}}{{ :clust_complete.png |}} |
| + | Est-ce que les deux dendrogrammes sont très différents? | ||
| - | **Défi 2** | + | On remarque que pour le groupement à liens simple, plusieurs objets s’enchaînent (par exemple les sites 19, 29, 30, 20, 26, etc.) alors que des groupes plus distincts peuvent être observés dans le groupement à liens complets. |
| - | 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. | + | La méthode de Ward diffère légèrement des deux méthodes précédentes. Le critère utilisé est la méthode des moindres carrés (comme dans les modèles linéaires). Ainsi, des objets/groupes fusionnent de façon à minimise la variance intragroupes. Pour débuter, chaque objet est considéré comme un groupe. À chaque étape, la paire de groupes à fusionner est celle qui résulte à la plus petite augmentation de la somme des carrés des écarts intra-groupes. |
| - | **Défi 2**: Solution | + | La méthode de Ward est également disponible sous la fonction hclust(). Par contre, le dendrogramme produit par défaut montre les distances au carré. Afin de comparer ce dendrogramme à celui du groupement à liens simples et à liens complets, il faut calculer la racine carrée des distances. |
| - | <hidden> | + | <code rsplus | hclust, méthode de Ward> |
| + | #Faire le groupement de Ward | ||
| + | spe.dhel.ward<-hclust(spe.dhel, method="ward.D2") | ||
| + | plot(spe.dhel.ward) | ||
| - | Votre code pourrait ressembler à ceci: | + | #Refaire le dendrogramme en utilisant la racine carrée des distances |
| + | spe.dhel.ward$height<-sqrt(spe.dhel.ward$height) | ||
| + | plot(spe.dhel.ward) | ||
| + | plot(spe.dhel.ward, hang=-1) # hang=-1 permet d’afficher les objets sur la même ligne | ||
| + | </code> | ||
| - | <code rsplus | RDA partielle sur les données d'acariens> | + | {{ :clust_ward.png |}} |
| + | {{ :clust_wardfinal.png |}} | ||
| - | mite.spe.subs=rda(mite.spe.hel ~ Shrub + Topo | + | Les groupements générés par la méthode de Ward ont tendance à être plus sphériques et à contenir des quantités plus similaires d’objets. |
| - | + Condition(SubsDens + WatrCont + Substrate), data=mite.env) | + | |
| - | + | ||
| - | #Extraire les résultats | + | |
| - | summary(mite.spe.subs, display=NULL) | + | |
| - | (R2adj <- RsquareAdj(mite.spe.subs)$adj.r.squared) | + | |
| - | + | ||
| - | #Axes significatifs | + | |
| - | anova.cca(mite.spe.subs, step=1000) | + | |
| - | anova.cca(mite.spe.subs, step=1000, by="axis") | + | |
| - | + | ||
| - | #Triplot cadrage 1 | + | |
| - | windows(title="Partial RDA scaling 1") | + | |
| - | plot(mite.spe.subs, 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(mite.spe.subs, display="sites", choices=c(1,2), scaling=1), | + | |
| - | pch=21, col="black", bg="steelblue", cex=1.2) | + | |
| - | text(scores(mite.spe.subs, display="species", choices=c(1), scaling=1), | + | |
| - | scores(mite.spe.subs, display="species", choices=c(2), scaling=1), | + | |
| - | labels=rownames(scores(mite.spe.subs, display="species", scaling=1)), | + | |
| - | col="grey", cex=0.8) | + | |
| - | arrows(0,0, | + | |
| - | scores(mite.spe.subs, display="bp", choices=c(1), scaling=1), | + | |
| - | scores(mite.spe.subs, display="bp", choices=c(2), scaling=1), | + | |
| - | col="red") | + | |
| - | text(scores(mite.spe.subs, display="bp", choices=c(1), scaling=1)+0.05, | + | |
| - | scores(mite.spe.subs, display="bp", choices=c(2), scaling=1)+0.05, | + | |
| - | labels=rownames(scores(mite.spe.subs, display="bp", choices=c(2), scaling=1)), | + | |
| - | col="red", cex=1) | + | |
| - | + | ||
| - | #Triplot cadrage 2 | + | |
| - | windows(title="Partial RDA scaling 2") | + | |
| - | plot(mite.spe.subs, 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(mite.spe.subs, display="sites", choices=c(1,2), scaling=2), | + | |
| - | pch=21, col="black", bg="steelblue", cex=1.2) | + | |
| - | text(scores(mite.spe.subs, display="species", choices=c(1), scaling=2), | + | |
| - | scores(mite.spe.subs, display="species", choices=c(2), scaling=2), | + | |
| - | labels=rownames(scores(mite.spe.subs, display="species", scaling=2)), | + | |
| - | col="grey", cex=0.8) | + | |
| - | arrows(0,0, | + | |
| - | scores(mite.spe.subs, display="bp", choices=c(1), scaling=2), | + | |
| - | scores(mite.spe.subs, display="bp", choices=c(2), scaling=2), | + | |
| - | col="red") | + | |
| - | text(scores(mite.spe.subs, display="bp", choices=c(1), scaling=2)+0.05, | + | |
| - | scores(mite.spe.subs, display="bp", choices=c(2), scaling=2)+0.05, | + | |
| - | labels=rownames(scores(mite.spe.subs, display="bp", choices=c(2), scaling=2)), | + | |
| - | col="red", cex=1) | + | |
| - | </code> | + | |
| - | {{ :partial_rda_ch6.png |}} | + | Quelle méthode choisir? |
| - | La RDA partielle est significative (p<0.001) de même que les deux premiers axes d’ordination. Les variables environnementales expliquent 9.81% de la variation de l’abondance en poissons, les variables relatives au substrat expliquent 42.84% de la variation et la variance non expliquée est de 47.35%. Le R2 ajusté de cette RDA est de 8.33%. | + | Le choix de la bonne mesure d’association et de la bonne méthode de groupement dépend de l’objectif. Qu’est-ce qu’il est plus intéressant de démontrer : des gradients? des contrastes? Il est également important de tenir en compte les propriétés de la méthode utilisée dans l’interprétation des résultats. Si plus d’une méthode semble adéquate pour répondre à une question biologique, comparer les dendrogrammes serait une bonne option. Encore une fois, le groupement n’est pas une analyse statistique, mais il est possible de tester les résultats et d’identifier des partitions ayant un sens biologique. Il est également possible de déterminer le nombre de groupes optimal et de performer des tests statistiques sur les résultats. Les méthodes de groupement peuvent aussi être combinées à une ordination pour distinguer des groupes de sites. Ces avenues ne seront pas explorées dans cet atelier. Pour aller plus loin, consulter Borcard et al. 2011. |
| - | </hidden> | + | |
| - | =====2.3 Partitionnement de la variation par RDA partielles===== | + | ======3. Ordination sans contrainte====== |
| - | 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 : | + | Les analyses d'ordination sans contrainte permettent d'organiser des échantillons, des sites ou des espèces le long de gradients continus (ex. écologiques ou environnementaux). Les ordinations sans contrainte se différencient des analyses canoniques (voir plus loin dans cet atelier "analyses canoniques") par le fait que ces techniques ne tentent pas de définir une relation entre des variables dépendantes et indépendantes. |
| - | - 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. | + | L'ordination sans contrainte peut être utilisée pour: |
| + | - Évaluer les relations //au sein// d'un ensemble de variables (et non pas entre séries de variables). | ||
| + | - Trouver les éléments clés de variation au sein d'échantillons, de sites, d'espèces, etc. | ||
| + | - Réduire les dimensions d'un jeu de données multivariées sans perte importante d'informations. | ||
| + | - Créer de nouvelles variables pour une utilisation dans des analyses ultérieures (comme la régression). Les composantes principales des axes d'ordination sont en effet des combinaisons linéaires des variables d'origine. | ||
| + | [[http://www.umass.edu/landeco/teaching/multivariate/schedule/ordination1.pdf|Source]] | ||
| - | {{ :vennd_varpart.png |}} | + | =====3.1 Analyses en Composantes Principales (Principal Component Analysis, PCA)===== |
| - | 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) | + | L'Analyse en Composantes Principales (ou PCA) fur originellement décrite par Pearson (1901) bien qu'elle soit le plus souvent attribué à Hotelling (1933) qui l'a proposé indépendamment. Cette méthode ainsi que nombre de ces implications sont pour l'analyse de données sont présentées dans l'article fondateur de Rao (1964). La PCA est utilisée pour générer, à partir d'un large jeu de données, un nombre restreint de variables clefs qui permettent de représenter au maximum la variation présente dans le jeu de données. En d'autres termes, la PCA est utilisée pour générer des combinaisons de variables à partir d'un ensemble plus grand de variables tout en conservant la majorité de variation de l'ensemble des données. La PCA est une technique d'analyse puissante pour l'analyse dans descripteurs quantitatifs (tels que les abondances d'espèces), mais ne peut pas être appliquée aux données binaires (telles que l'absence/présence des espèces). |
| - | 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. | + | À partir d'un jeu de données contenant des variables à distribution normale, le premier axe de PCA (ou axe de composante principale) correspond à la droite qui traverse la plus grande dimension de l’ellipsoïde décrivant la distribution multi-normale des données. Les axes suivants traversent de façon similaire cet ellipsoïde selon l'ordre décroissant de ces dimensions. Ainsi, il est possible d'obtenir un maximum de p axes principaux à partir d'un jeu de données contenant p variables. |
| - | 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. | + | Pour cela, la PCA effectue un rotation du système d'axes originels défini par les variables de façon à ce que les axes d'ordination successifs soient orthogonaux entre eux et correspondent aux dimensions successives du maximum de variance observée dans le nuage de points (voir ci-dessus). Les nouvelles variables produites par la PCA sont non-corrélées entre elles (les axes d'ordination étant orthogonaux) et peuvent alors être utilisés dans d'autres types d'analyse telles que des régressions multiples (Gotelli et Ellison, 2004). Les composantes principales situent la position des objets dans le nouveau système de coordonnées calculé par la PCA. La PCA s'effectue sur une matrice d'association entre variables, et a pour caractéristique de préserver les distances euclidiennes et de détecter des relations linéaires, uniquement. En conséquence, les abondances brutes des espèces doivent être soumises à une pré-transformation (comme une transformation d'Hellinger) avant de réaliser une PCA. |
| + | //Pour faire une PCA, vous avez besoin :// | ||
| + | - Un ensemble de variables (sans distinction entre les variables indépendantes ou dépendantes, c-est-à-dire un ensemble d'espèces OU un ensemble de variables environnementales). | ||
| + | - De sites (objets, échantillons) dans lesquels sont mesurés les mêmes variables. | ||
| + | - En général, il est préférable d'avoir un plus grand nombre de sites que de variables dans le jeu de données (plus de lignes que de colonnes). | ||
| - | <code rsplus | varpart()> | + | La PCA est particulièrement utile pour des matrices de données contenant plus de deux variables, mais il est plus facile de décrire son fonctionnement avec un exemple bidimensionnel. Dans l'exemple suivant (d'après Clarke et Warwick 2001), la matrice contient les données d'abondance de deux espèces dans neuf sites: |
| - | ?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) | + | |
| - | </code> | + | |
| + | ^Site^Espèce 1^Espèce 2^ | ||
| + | ^A^6^2^ | ||
| + | ^B^0^0^ | ||
| + | ^C^5^8^ | ||
| + | ^D^7^6^ | ||
| + | ^E^11^6^ | ||
| + | ^F^10^10^ | ||
| + | ^G^15^8^ | ||
| + | ^H^18^14^ | ||
| + | ^I^14^14^ | ||
| - | Le sommaire ressemble à ceci: | + | La représentation des sites en deux dimensions devrait ressembler à ceci: |
| - | {{ :varpart_output.png |}} | + | {{:pcaex_1.png?500|}} |
| - | {{ :varpart_output_venn.png |}} | + | |
| - | 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. | + | Ce nuage de points est une ordination. Il présente la distribution des espèces entre les sites, mais vous pouvez imaginer qu'il est plus difficile de visualiser un tel graphique en présence de plus de deux espèces. Dans ce cas, l'objectif est de réduire le nombre de variables en composantes principales. Pour réduire les données bidimensionnelles précédentes à une dimension, une PCA peut être effectuée : |
| - | Il est également possible d’effectuer un partitionnement de la variation à partir de jeux de données explicatives contenant uniquement des variables significatives : | + | {{:pcaex_2.png?500|}} |
| + | Dans ce cas, la première composante principale est orientée dans le sens de la plus grande variation dans les points, ces points étant perpendiculaires à la ligne. | ||
| - | <code rsplus | varpart avec variables explicatives> | + | Une seconde composante principale est alors ajoutée perpendiculairement à la première: |
| - | #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] | + | |
| - | </code> | + | |
| + | {{:pcaex_3.png?500|}} | ||
| - | 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). | + | Dans le diagramme final, les deux axes de PCA sont pivotés et les axes sont maintenant les composantes principales (et non plus les espèces): |
| - | **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 ? | + | {{:pcaex_4.png?500|}} |
| - | + | ||
| - | **Défi 3**: Solution | + | |
| - | <hidden> | + | Pour les PCA avec plus de deux variables, les composantes principales sont ajoutées de la façon suivante (Clarke et Warwick 2001): |
| - | Votre code pourrait ressembler à ceci: | + | |
| - | <code rsplus | varpart avec variables significatives> | + | PC1 = axe qui maximise la variance des points qui sont projetées perpendiculairement à l'axe. |
| - | str(mite.env) | + | PC2 = axe perpendiculaire à PC1, mais dont la direction est à nouveau celle maximisant la variance lorsque les points y sont projetés perpendiculairement. |
| - | (mite.subs=mite.env[,c(1,2,3)]) #Premier jeu de données | + | PC3 et ainsi de suite: perpendiculaire aux deux premiers axes dont la direction est à nouveau celle maximisant la variance lorsque les points y sont projetés perpendiculairement. |
| - | (mite.other=mite.env[,c(4,5)]) #Deuxième jeu de données | + | |
| - | + | ||
| - | #RDA sur mite.subs | + | |
| - | rda.mite.subs <- rda(mite.spe.hel~., data=mite.subs) | + | |
| - | R2a.all.chem <- RsquareAdj(rda.mite.subs)$adj.r.squared | + | |
| - | + | ||
| - | #Sélection progressive pour mite.subs | + | |
| - | ordiR2step(rda(mite.spe.hel~1, data=mite.subs), | + | |
| - | scope= formula(rda.mite.subs), direction= "forward", R2scope=TRUE, pstep=1000) | + | |
| - | names(mite.subs) | + | |
| - | (mite.subs.pars <- mite.subs[, c(2, 3)]) | + | |
| - | + | ||
| - | #RDA sur mite.other | + | |
| - | rda.mite.other <- rda(mite.spe.hel~., data=mite.other) | + | |
| - | R2a.all.chem <- RsquareAdj(rda.mite.other)$adj.r.squared | + | |
| - | + | ||
| - | #Sélection progressive sur mite.other | + | |
| - | ordiR2step(rda(mite.spe.hel~1, data=mite.other), | + | |
| - | scope= formula(rda.mite.other), direction= "forward", R2scope=TRUE, pstep=1000) | + | |
| - | names(mite.other) | + | |
| - | (mite.other.pars <- mite.other[, c(1,2)]) | + | |
| - | + | ||
| - | #Partitionnement de la variation | + | |
| - | (mite.spe.part <- varpart(mite.spe.hel, ~WatrCont+Substrate, ~Shrub+Topo, | + | |
| - | data=mite.env)) | + | |
| - | windows(title="Variation partitioning - parsimonious subsets") | + | |
| - | plot(mite.spe.part, digits=2) | + | |
| - | + | ||
| - | # Test de significativité des fractions | + | |
| - | anova.cca(rda(mite.spe.hel~ WatrCont+Substrate, data=mite.env), step=1000) # Test of fractions [a+b] | + | |
| - | anova.cca(rda(mite.spe.hel~Shrub+Topo, data=mite.env), step=1000) # Test of fractions [b+c] | + | |
| - | (env.pars <- cbind(mite.env[,c(2,3,4,5)])) | + | |
| - | anova.cca(rda(mite.spe.hel~ WatrCont+Substrate+Shrub+Topo, data=env.pars), step=1000) # Test of fractions [a+b+c] | + | |
| - | anova.cca(rda(mite.spe.hel~WatrCont+Substrate + Condition(Shrub+Topo), data=env.pars), step=1000) # Test of fraction [a] | + | |
| - | anova.cca(rda(mite.spe.hel~Shrub+Topo+ Condition(WatrCont+Substrate ), data=env.pars), step=1000) # Test of fraction [c] | + | |
| - | </code> | + | |
| - | Dans ce cas, les variables relatives au substrat expliquent 14.00% de la variation en espèces tandis que les autres variables expliquent 9.1% de la variation. L’interaction de ces deux types de variables explique 16.90% de la variation. Toutes ces fractions sont significatives (p<0.001). | + | Lorsqu'il y a plus de deux dimensions, la PCA produit un nouvel espace dans lequel tous les axes de PCA sont orthogonaux (ce qui signifie que la corrélation entre chaque combinaison de deux axes est nulle), et où les axes de PCA sont ordonnés selon la proportion de variance des données d'origine qu'ils expliquent. |
| - | </hidden> | + | Les données "spe" comprennent 27 espèces de poissons. Pour simplifier cette diversité à un petit nombre de variables ou pour identifier différents groupes de sites associés à des espèces particulières, une PCA peut être effectuée. |
| + | Exécuter une PCA sur les données d'abondance d'espèces soumises à la transformation d'Hellinger : | ||
| + | <code rsplus | PCA avec rda() de (vegan)> | ||
| + | #Exécuter la PCA avec la fonction rda()- cette fonction calcule à la fois des PCA et des RDA | ||
| + | spe.h.pca<-rda(spe.hel) | ||
| - | ======3. Arbre de régression multivarié====== | + | #Extraire les résultats |
| + | summary(spe.h.pca) | ||
| + | </code> | ||
| - | 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. | + | Résultats: |
| - | Avantages du MRT par rapport à la RDA: | + | {{:pca_outputfr_1.png?800|}} |
| - | * n'assume pas de relation linéaire entre les abondances d'espèces et les variables environnementales (quantitatives ou qualitatives), | + | {{:pca_outputfr_2.png?800|}} |
| - | * la méthode est robuste en présence de valeurs manquantes, | + | {{:pca_outputfr_3.png?800|}} |
| - | * 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: | + | **Interprétation des résultats de PCA** |
| + | Cette sortie R contient les valeurs propres ou «eigenvalue» de la PCA. La valeur propre est la valeur de la variation ramenée à la longueur d'un vecteur, et correspond à la quantité de variation expliquée par chaque axe d'ordination de la PCA. Comme vous pouvez le voir, la fonction summary fournit de nombreuses informations. Parmi les résultats, la proportion de variance des données expliquée par les variables sans contraintes est une information importante. Dans cet exemple, la variance totale des sites expliquée par les espèces est de 0,5 (50%). Le résumé vous indique également quelle proportion de la variance totale expliquée est répartie entre chaque composantes principales de la PCA: le premier axe de PCA explique 51.33% de la variation tandis que le second axe explique 12.78%. Vous pouvez également extraire certaines parties des résultats: | ||
| + | <code rsplus | Résultats de PCA> | ||
| + | summary(spe.h.pca, display=NULL) # seulement les valeurs propres | ||
| + | eigen(cov(spe.hel)) # vous pouvez aussi trouver les valeurs propres par cette ligne de code | ||
| + | </code> | ||
| - | Feuille: Groupe terminal de sites | + | Les scores (c'est-à-dire les coordonnées) des sites ou des espèces peuvent également être extraits d'une PCA. Ces scores permettent, par exemple, d'utiliser une composante principale comme une variable dans une autre analyse, ou de faire des graphiques supplémentaires. Par exemple, vous pouvez obtenir par PCA une variable unique issue du jeu de données "spe" puis l'utiliser pour la corréler par régression à une autre variable, ou déterminer un gradient spatial. Pour extraire les scores d'une PCA, utiliser la fonction scores (): |
| + | <code rsplus | scores()> | ||
| + | spe.scores<-scores(spe.h.pca, display="species", choices=c(1,2)) # scores des espèces selon les premier et deuxième axes | ||
| + | site.scores<-scores(spe.h.pca, display="sites", choices=c(1,2)) # scores des sites selon les premier et deuxième axes | ||
| + | #Remarque: si vous ne spécifiez pas le nombre de composantes principales à l'aide de choices = c (1,2) | ||
| + | #(ou choices = c (1: 2)), les scores selon toutes les composantes principales seront extraits. | ||
| + | </code> | ||
| - | Noeud: Point où les données se divisent en deux groupes. Chaque noeud est caractérisé par une valeur d'une variable explicative. | + | La PCA des données d'abondances de poissons produit autant de composantes principales qu'il y d'espèces (i.e. de colonnes dans le jeu de données), soit 27 composantes principales. Le nombre de variables à traiter n'est donc pas directement réduit par la PCA. Pour réduire le nombre de variables, il est alors nécessaire de déterminer quelles composantes principales sont significatives et doivent être conservées, par exemple à l'aide du critère de Kaiser-Guttman. Ce critère compare la variance expliquée par chaque composante principale à la moyenne de la variance expliquée par l'ensemble des composantes principales. Un histogramme illustrant la significativité des différentes composantes principale peut ensuite être tracé à l'aide du code ci-dessous : |
| + | <code rsplus | Les axes> | ||
| + | # Identification des axes significatifs de la PCA à l'aide du critère de Kaiser-Guttman | ||
| + | ev<-spe.h.pca$CA$eig | ||
| + | ev[ev>mean(ev)] | ||
| + | n<-length(ev) | ||
| + | bsm<-data.frame(j=seq(1:n), p=0) | ||
| + | bsm$p[1]=1/n | ||
| + | for (i in 2:n) { | ||
| + | bsm$p[i]=bsm$p[i-1]+(1/(n=1-i))} | ||
| + | bsm$p=100*bsm$p/n | ||
| + | bsm | ||
| + | barplot(ev, main="valeurs propres", col="grey", las=2) | ||
| + | abline(h=mean(ev), col="red") | ||
| + | legend("topright", "moyenne des valeurs propres", lwd=1, col=2, bty="n") | ||
| + | </code> | ||
| - | Branche: Chaque lignée formée par un noeud. | + | {{:pca_sigaxes_sp.png?500|}} |
| + | Cet histogramme montre que la proportion de la variance expliquée par chaque composante chute en-dessous de la proportion moyenne expliquée par l'ensemble des composantes après le sixième axe d"ordination PC6. En consultant le résumé de nouveau, on peut constater que la proportion cumulée de la variance expliquée par les cinq premières composantes principales est de 85%. | ||
| - | 1- Partitionnement des données sous contrainte (construction de l'arbre) | + | La PCA n'est pas seulement appropriée pour les données de composition d'espèces, mais peut également être exécutée sur des variables environnementales standardisées: |
| + | <code rsplus | PCA pour les variables environnementales> | ||
| + | #Exécuter la PCA | ||
| + | env.pca<-rda(env.z) # ou rda(env, scale=TRUE) | ||
| - | 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. | + | #Extraction des résultats |
| + | summary(env.pca) | ||
| + | summary(env.pca, scaling=2) | ||
| + | </code> | ||
| + | «Scaling» réfère à quelle portion de la PCA est redimensionnée aux valeurs propres. Scaling = 2 signifie que les scores des espèces sont mises à l'échelle des valeurs propres, alors que scaling = 1 signifie que les scores des sites sont mises à l'échelle des valeurs propres. Scaling = 3 signifie qu'à la fois les scores des espèces et des sites sont mis symétriquement à l'échelle de la racine carrée des valeurs propres. En scaling = 1,les distances euclidiennes entre les sites (lignes de la matrice de données) sont conservées tandis qu'en scaling = 2 les corrélations entre espèces (les colonnes de la matrice de données) sont conservées. Cela implique que lorsque vous regardez un biplot de PCA en Scaling = 2, l'angle entre les descripteurs représente leur corrélation. | ||
| - | 2- Validation croisée et élagage de l'arbre | + | <code rsplus | Axes significatifs> |
| + | ev<-env.pca$CA$eig | ||
| + | ev[ev>mean(ev)] | ||
| + | n<-length(ev) | ||
| + | bsm<-data.frame(j=seq(1:n), p=0) | ||
| + | bsm$p[1]=1/n | ||
| + | for (i in 2:n) { | ||
| + | bsm$p[i]=bsm$p[i-1]+(1/(n=1-i))} | ||
| + | bsm$p=100*bsm$p/n | ||
| + | bsm | ||
| + | barplot(ev, main="valeurs propres", col="grey", las=2) | ||
| + | abline(h=mean(ev), col="red") | ||
| + | legend("topright", "moyenne des valeurs propres", lwd=1, col=2, bty="n") | ||
| + | </code> | ||
| + | Comparer cet histogramme de valeurs propres avec celui que vous avez créé pour la PCA sur les abondances d'espèces. | ||
| - | 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. | + | Bien que beaucoup d'informations puissent être extraites d'une PCA par la fonction summary de la PCA, l'interpétation et la communication des résultats est souvent facilitée en traçant un biplot. Sur un biplot de PCA, l'axe des x correspond à la première composante principale et l'axe des y à la deuxième composante principale. |
| - | Créons un arbre de régression multivarié sur les données de la rivière Doubs. | + | La fonction plot () permet de tracer des biplot sur lequels les sites figurent en chiffres noirs et les espèces sont représentées en rouge. Le biplot de la //PCA des abondances d'espèces// peut être appelé comme suit: |
| + | <code rsplus | Graphique PCA des abondances d'espèces > | ||
| + | plot(spe.h.pca) | ||
| + | </code> | ||
| - | <code rsplus | mvpart()> | + | {{:spe_PCA1.png?800|}} |
| - | ?mvpart | + | |
| - | + | ||
| - | #Preparer les données: enlever la variable “distance from source” | + | |
| - | env <- subset(env, select = -das) | + | |
| - | # Créer l'arbre | + | La construction de biplots de PCA s'articule en trois étapes: |
| - | doubs.mrt<-mvpart(data.matrix(spe.hel)~.,env, | + | <code rsplus | Biplot de PCA> |
| - | legend=FALSE,margin=0.01,cp=0,xv="pick",xval=nrow(spe.hel),xvmult=100,which=4) | + | plot(spe.h.pca, type=”n”) # Produit une figure vierge |
| + | points(spe.h.pca, dis=”sp”, col=”blue”) # ajoute les points correspondant aux espèces | ||
| + | #utilizer text() plutôt que points() si vous préférez que les codes des espèces s'affichent (nom des colonnes) | ||
| + | points(spe.h.pca, dis=”sites”, col=”red”) # ajoute les points correspondant aux sites | ||
| </code> | </code> | ||
| + | {{:spe_PCA2.png?800|}} | ||
| - | 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). | + | Pour créer de plus beaux biplots, essayez ce code: |
| + | <code rsplus | Graphique PCA des abondances d'espèces v.2> | ||
| + | #Scaling 1 | ||
| + | windows() | ||
| + | plot(spe.h.pca) | ||
| + | windows() | ||
| + | biplot(spe.h.pca) | ||
| + | windows() | ||
| + | # scaling 1 = distance biplot : | ||
| + | # distances entre les objets est une approximation de leur distance euclidienne | ||
| + | # les angles entre les descripteurs ne réflètent PAS leur corrélation | ||
| + | plot(spe.h.pca, scaling=1, type="none", | ||
| + | xlab<-c("PC1 (%)", round((spe.h.pca$CA$eig[1]/sum(spe.h.pca$CA$eig))*100,2)), | ||
| + | ylab<-c("PC2 (%)", round((spe.h.pca$CA$eig[2]/sum(spe.h.pca$CA$eig))*100,2))) | ||
| + | points(scores(spe.h.pca, display="sites", choices=c(1,2), scaling=1), | ||
| + | pch=21, col="black", bg="steelblue", cex=1.2) | ||
| + | text(scores(spe.h.pca, display="species", choices=c(1), scaling=1), | ||
| + | scores(spe.h.pca, display="species", choices=c(2), scaling=1), | ||
| + | labels=rownames(scores(spe.h.pca, display="species", scaling=1)), | ||
| + | col="red", cex=0.8) | ||
| + | </code> | ||
| + | Le code ci-dessus a produit trois biplots mais le dernier est le plus attrayant : | ||
| + | {{:spe_PCA3.png?800|}} | ||
| - | {{ :cross_validation.png |}} | + | Sur ce graphique, les sites scores sont indiqués par des points bleus et les noms d'espèces sont en rouge. Il est également possible de représenter les sites par leurs noms. |
| - | 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. | + | Comment interpréter ce type de graphique ? |
| + | Ce biplot permet d'observer qu'il existe certains groupes de sites homogènes du point de vue de la composition de leur communautés de poissons. On y voit également que l'espèce «ABL» n'a pas la même prévalence dans la majorité des sites que les autres espèces plus proches du centre du graphique. | ||
| - | 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. | + | Les biplots ne doivent pas seulement être interprétés en termes de proximité, mais également d'angles. Deux variables séparées d'un angle de 90 degrés ne sont pas corrélées. Deux variables très rapprochés sont fortement corrélées. Deux variables aux directions opposées sont corrélées négativement. |
| + | Maintenant regardons le biplot de la //PCA environnement//: | ||
| + | <code rsplus | Graphique PCA des variables environnementales> | ||
| + | #Scaling 2 | ||
| + | windows() | ||
| + | plot(env.pca) | ||
| + | windows() | ||
| + | # scaling 2 = graphique de corrélations : | ||
| + | # les distances entre les objets ne sont PAS des approximations de leur distance euclidienne | ||
| + | # les angles entres les descripteurs reflètent leur corrélation | ||
| + | plot(env.pca, scaling=2, type="none", | ||
| + | xlab<-c("PC1 (%)", round((env.pca$CA$eig[1]/sum(env.pca$CA$eig))*100,2)), | ||
| + | ylab<-c("PC2 (%)", round((env.pca$CA$eig[2]/sum(env.pca$CA$eig))*100,2)), | ||
| + | xlim<-c(-1,1), ylim=c(-1,1)) | ||
| + | points(scores(env.pca, display="sites", choices=c(1,2), scaling=2), | ||
| + | pch=21, col="black", bg="darkgreen", cex=1.2) | ||
| + | text(scores(env.pca, display="species", choices=c(1), scaling=2), | ||
| + | scores(env.pca, display="species", choices=c(2), scaling=2), | ||
| + | labels<-rownames(scores(env.pca, display="species", scaling=2)), | ||
| + | col="red", cex=0.8) | ||
| + | </code> | ||
| - | {{ :mrt_1se.png |}} | + | {{:env_PCA1.png?800|}} |
| - | 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. | + | Rappelez-vous qu'un biplot de PCA est en fait un nuage de points dans lequel les axes sont des combinaisons linéaires des variables d'origine. Il existe donc beaucoup de façons différentes de tracer un biplot. Par exemple, vous pouvez utiliser la fonction ggplot () et les compétences acquises de l'atelier 3 pour tracer votre graphique d'ordination dans ggplot. |
| - | 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). | + | **Utilisation des axes de PCA comme variable explicative composite** |
| + | Dans certains cas, l'utilisateur cherche à réduire un grand nombre de variables environnementales en un plus faible nombre de variables composites. Lorsque les axes de PCA représentent des gradients écologiques (i.e. lorsue les variables environnementales sont corrélées de façon cohérente avec les axes de PCA), l'utlisateur peut utiliser les scores des sites le long axes de PCA dans de nouvelles analyses (au lieu d'utiliser les variables environnementales brutes). En d'autres termes, étant donné que les scores des sites le long des axes de PCA représentent des combinaisons linéaires des descripteurs, ils peuvent être utilisés comme proxy des conditions écologiques dans de nouvelles analyses. | ||
| - | <code rsplus | Comparaison des arbres> | + | Dans l'exemple ci-dessus, le premier axe de PCA peut être identifié comme un gradient écologique allant des sites oligotrophes riches en oxygène aux sites eutrophes pauvres en oxygène: de gauche à droite, le premier groupe de sites montrent les plus hautes altitudes (alt) et pente (slope), et les plus faibles débit (deb) et distance à la source (das). Le second groupe de sites possède les plus hautes valeurs de concentration en oxygène (oxy) et les plus faibles concentrations en nitrates (nit). Un troisième groupe de sites montrent des valeurs intermédiaire pour l'ensemble de ces variables. |
| - | # 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 | + | Dans ce cas, si l'objectif est d'identifier si une espèce particulière est associée au gradient oligotrophe-eutrophe, il est possible de corréler l'abondance de cette espèce aux scores des sites le long du premier axe de PCA. Par exemple, si l'utilisateur veut identifier si l'espèce TRU est associée à des eaux oligotrophes ou eutrophes, il lui est possible d’utiliser le modèle linéaire suivant: |
| - | 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) | + | |
| - | </code> | + | |
| - | {{ :mrt_cvre.png |}} | + | <code rsplus | Graphique PCA des variables environnementales> |
| - | {{ :mrt_4.png |}} | + | Sites_scores_Env_Axis1<- scores(env.pca, display="sites", choices=c(1), scaling=2) |
| + | spe$ANG | ||
| + | plot( Sites_scores_Env_Axis1, spe$TRU) | ||
| + | summary(lm(spe$TRU~Sites_scores_Env_Axis1)) | ||
| + | abline(lm(spe$TRU~Sites_scores_Env_Axis1)) | ||
| + | </code> | ||
| - | 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. | + | Ce modèle simple montre que l'abondance de l'espèce TRU est significativement liée aux scores des sites le long du premier axe de PCA (t = -5.30, p = 1.35e-05, adj-R2 = 49.22%), c’est-à-dire qu'elle dépend d'un gradient oligotrophe-eutrophe. dans ce cas l'espèce TRU préfère donc les eaux oligotrophes. |
| - | Le sommaire de la fonction contient beaucoup d'informations. | + | **Défi 3** |
| + | Exécuter une PCA sur l'abondance des espèces de mites. Quels sont les axes significatifs ? Quels sont groupes de sites pouvez-vous identifier? Quelles espèces sont liées à chaque groupe de sites? | ||
| + | <code rsplus | Données mites> | ||
| + | mite.spe<-data(mite) # données disponibles dans vegan | ||
| + | </code> | ||
| + | **Défi - Solution** | ||
| + | <hidden> | ||
| + | Votre code ressemble certainement à celui-ci: | ||
| + | <code rsplus | Mite PCA> | ||
| - | <code rsplus | sommaire du MRT> | + | # Transformation de Hellinger |
| - | summary(doubs.mrt) | + | mite.spe.hel<-decostand(mite.spe, method="hellinger") |
| - | </code> | + | mite.spe.h.pca<-rda(mite.spe.hel) |
| + | |||
| + | # Quels sont les axes significatifs? | ||
| + | ev<-mite.spe.h.pca$CA$eig | ||
| + | ev[ev>mean(ev)] | ||
| + | n<-length(ev) | ||
| + | bsm<-data.frame(j=seq(1:n), p=0) | ||
| + | bsm$p[1]=1/n | ||
| + | for (i in 2:n) { | ||
| + | bsm$p[i]=bsm$p[i-1]+(1/(n=1-i))} | ||
| + | bsm$p=100*bsm$p/n | ||
| + | bsm | ||
| + | barplot(ev, main="Valeurs propres", col="grey", las=2) | ||
| + | abline(h=mean(ev), col="red") | ||
| + | legend("topright", "Moyenne des valeurs propres", lwd=1, col=2, bty="n") | ||
| - | {{ :doubs_mrt_summary.png |}} | + | # Résultats |
| + | summary(mite.spe.h.pca, display=NULL) | ||
| + | windows() | ||
| - | 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. | + | # Représentation graphique de la PCA |
| + | plot(mite.spe.h.pca, scaling=1, type="none", | ||
| + | xlab=c("PC1 (%)", round((mite.spe.h.pca$CA$eig[1]/sum(mite.spe.h.pca$CA$eig))*100,2)), | ||
| + | ylab=c("PC2 (%)", round((mite.spe.h.pca$CA$eig[2]/sum(mite.spe.h.pca$CA$eig))*100,2))) | ||
| + | points(scores(mite.spe.h.pca, display="sites", choices=c(1,2), scaling=1), | ||
| + | pch=21, col="black", bg="steelblue", cex=1.2) | ||
| + | text(scores(mite.spe.h.pca, display="species", choices=c(1), scaling=1), | ||
| + | scores(mite.spe.h.pca, display="species", choices=c(2), scaling=1), | ||
| + | labels=rownames(scores(mite.spe.h.pca, display="species", scaling=1)), | ||
| + | col="red", cex=0.8) | ||
| + | </code> | ||
| + | Votre graphique ressemblera à ceci. | ||
| + | {{:mite_pca.png?800|}} | ||
| - | <code rsplus | Trouver des espèces discriminantes et indicatrices> | + | Bien que les sites soient tous de composition semblable (aucun groupe distinct de sites n'apparait sur le biplot), certaines espèces semblent souvent être présentes ensemble, par exemple Spec 01, Spec 10, Spec 14 et Spec 15. |
| - | # Trouver les espèces discriminantes de l'arbre de régression | + | </hidden> |
| - | doubs.mrt.wrap<-MRT(doubs.mrt,percent=10,species=colnames(spe.hel)) | + | |
| - | summary(doubs.mrt.wrap) | + | |
| - | # Extraire les p-values des valeurs indval | + | =====3.2. Analyse des Correpondances (Correspondence Analysis, CA)===== |
| - | doubs.mrt.indval<-indval(spe.hel,doubs.mrt$where) | + | |
| - | doubs.mrt.indval$pval | + | |
| - | # Extraire les espèces indicatrices à chaque noeud et leur valeur indval | + | L'une des hypothèses clefs de la PCA postule que les espèces sont liées les unes aux autres de façon linéaire, et qu'elles répondent de façon linéaire aux gradients écologiques. Ce n'est cependant pas nécessairement le cas dans les données écologiques (e.g. beaucoup d'espèces montrent en effet un distribution unimodale le long des gradients environnementaux). Utiliser une PCA avec des données contenant des espèces à distribution unimodale, ou un grand nombre de zéros (absence des espèces), peut conduire à un phénomène statistique appelé "horseshoe effect" (ou effet fer à cheval) se produisant le long de gradients écologiques. Dans de tels cas, l'Analyse des Correspondances (CA) permet de mieux représenter les données (voir Legendre et Legendre pour plus d'informations). Comme la CA préserve les distances de Chi2 entre objets (tandis que la PCA préserve les distances euclidiennes), cette technique est, en effet, plus appropriée pour ordonner les jeux de données contenant des espèces à distribution unimodale, et a, pendant longtemps, était l'une des techniques les plus employées pour analyser les données d'absence-présence ou d'abondances d'espèces. Lors d'une CA, les données brutes sont d'abord transformées en une matrice Q des contributions cellule-par-cellule à la statistique Chi2 de Pearson, puis la matrice résultante est soumise à une décomposition en valeurs singulières afin de calculer les valeurs propres et vecteurs propres de l'ordination. |
| - | doubs.mrt.indval$maxcls[which(doubs.mrt.indval$pval<=0.05)] | + | |
| - | doubs.mrt.indval$indcls[which(doubs.mrt.indval$pval<=0.05)] | + | |
| - | </code> | + | |
| - | {{ ::doubs_mrt_discriminant.png |}} | + | Le résultat d'une CA représente donc une ordination dans laquelle les distances de Chi2 entre objets sont préservées (au lieu de la distance euclidienne dans une PCA), le distance de Chi2 n'étant pas influencée par la présence de double-zéros. Ainsi, la CA constitue une méthode d’ordination puissante pour l'analyse des abondances brutes d'espèces (i.e. sans pré-transformation). Contrairement à la PCA, la CA peut être appliquée sur des données quantitatives ou binaires (telles que des abondance ou absence-présence d'espèces). Comme dans une PCA, le critère de Kaiser-Guttman peut être utilisé pour identifier les axes significatifs d'une CA, et les scores des objets le long des axes d'ordination peuvent être extraits pour être utlisés dans des régressions multiples par exemple. |
| - | {{ ::doubs_mrt_finalpart.png |}} | + | |
| - | 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. | + | Exécuter une CA sur les données d'abondance d'espèces: |
| - | {{ ::doubs_mrt_indval.png |}} | + | <code rsplus | CA par cca() de vegan> |
| + | #Effectuer une CA à l'aide de la fonction cca() (NB: cca() est utilisée à la fois pour les CA et CCA) | ||
| + | spe.ca <- cca(spe) | ||
| + | |||
| + | # Identifier les axes significatifs | ||
| + | ev<-spe.ca$CA$eig | ||
| + | ev[ev>mean(ev)] | ||
| + | n=length(ev) | ||
| + | barplot(ev, main="Eigenvalues", col="grey", las=2) | ||
| + | abline(h=mean(ev), col="red") | ||
| + | legend("topright", "Average eigenvalue", lwd=1, col=2, bty="n") | ||
| + | </code> | ||
| - | 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). | + | {{ :ca_guttmankaiser.png |}} |
| - | **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’après cet histogramme, à partir du sixième axe d’ordination CA6, la proportion de variance expliquée diminue sous la proportion moyenne expliquée par l'ensemble des axes. La sortie R de la CA ci-dessous montre également que les cinq premiers axes d'ordination explique une proportion cumulée de variance expliquée de 84.63%. |
| - | **Défi 4** - Solution | + | <code rsplus | Extraire les résultats d'une CA> |
| + | summary(spe.h.pca) | ||
| + | summary(spe.h.pca, diplay=NULL) | ||
| + | </code> | ||
| - | <hidden> | + | {{ :ca_summary.png |}} |
| - | <code rsplus | Créé un MRT avec les données sur les acariens> | + | Les résultats d'une CA sont présentés sous R de la même façon que ceux d'une PCA. On y observe que le premier axe CA1 explique 51.50% de la variation de l'abondance des espèces tandis que le second axe CA2 explique 12.37% de la variation. |
| - | mite.mrt<-mvpart(data.matrix(mite.spe.hel)~.,mite.env, | + | |
| - | legend=FALSE,margin=0.01,cp=0,xv="pick", | + | |
| - | xval=nrow(mite.spe.hel),xvmult=100,which=4) | + | |
| - | summary(mite.mrt) | + | |
| - | mite.mrt.wrap<-MRT(mite.mrt,percent=10,species=colnames(mite.spe.hel)) | + | <code rsplus | Construction des biplots> |
| - | summary(mite.mrt.wrap) | + | par(mfrow=c(1,2)) |
| + | #### scaling 1 | ||
| + | plot(spe.ca, scaling=1, type="none", main='CA - biplot scaling 1', xlab=c("CA1 (%)", round((spe.ca$CA$eig[1]/sum(spe.ca$CA$eig))*100,2)), | ||
| + | ylab=c("CA2 (%)", round((spe.ca$CA$eig[2]/sum(spe.ca$CA$eig))*100,2))) | ||
| - | mite.mrt.indval<-indval(mite.spe.hel,mite.mrt$where) | + | points(scores(spe.ca, display="sites", choices=c(1,2), scaling=1), pch=21, col="black", bg="steelblue", cex=1.2) |
| - | mite.mrt.indval$pval | + | |
| - | mite.mrt.indval$maxcls[which(mite.mrt.indval$pval<=0.05)] | + | text(scores(spe.ca, display="species", choices=c(1), scaling=1), |
| - | mite.mrt.indval$indcls[which(mite.mrt.indval$pval<=0.05)] | + | scores(spe.ca, display="species", choices=c(2), scaling=1), |
| + | labels=rownames(scores(spe.ca, display="species", scaling=1)),col="red", cex=0.8) | ||
| + | |||
| + | #### scaling 2 | ||
| + | plot(spe.ca, scaling=1, type="none", main='CA - biplot scaling 2', xlab=c("CA1 (%)", round((spe.ca$CA$eig[1]/sum(spe.ca$CA$eig))*100,2)), | ||
| + | ylab=c("CA2 (%)", round((spe.ca$CA$eig[2]/sum(spe.ca$CA$eig))*100,2)), ylim=c(-2,3)) | ||
| + | |||
| + | points(scores(spe.ca, display="sites", choices=c(1,2), scaling=2), pch=21, col="black", bg="steelblue", cex=1.2) | ||
| + | text(scores(spe.ca, display="species", choices=c(1), scaling=2), | ||
| + | scores(spe.ca, display="species", choices=c(2), scaling=2), | ||
| + | labels=rownames(scores(spe.ca, display="species", scaling=2)),col="red", cex=0.8) | ||
| </code> | </code> | ||
| - | 25.6% de la variation dans la communauté d'acariens entre les sites est expliquée par le partitionnement des sites en fonction de la quantité d'eau dans le substrat (à 385.1 mg/L). LCIL est une espèces discriminante pour les sites à haute teneur en eau et sa valeur IndVal et de 0.715. | + | {{ :ca_biplot.png |}} |
| - | </hidden> | + | Ces biplots montrent qu'un groupe de sites (à gauche) possède des communautés similaires de poissons caractérisées par de nombreuses espèces dont GAR, TAN, PER, ROT, PSO et CAR; dans le coin supérieur droit, un second groupe de sites se caractérisent par les espèces LOC, VAI et TRU; le dernier groupe de sites dans le coin inférieur droit montrent des communautés abondantes en BLA, CHA et OMB. |
| + | **Défi 4** | ||
| + | Exécuter une CA sur les données d'abondance des //espèces d'acariens// (données mite). Quels sont les axes importants? Quels groupes de sites pouvez-vous identifier? Quelles espèces sont liées à chaque groupe de sites? | ||
| - | ======4. Analyse discriminante linéaire====== | + | **Défi 4 - Solution** |
| + | <hidden> | ||
| + | Votre code devrait s'apparenter à celui-ci: | ||
| - | 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. | + | <code rsplus | CA sur les données mite> |
| + | # CA | ||
| + | mite.spe.ca<-cca(mite.spe) | ||
| - | 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. | + | # Quels sont les axes importants? |
| + | ev<-mite.spe.ca$CA$eig | ||
| + | ev[ev>mean(ev)] | ||
| + | n=length(ev) | ||
| + | barplot(ev, main="Eigenvalues", col="grey", las=2) | ||
| + | abline(h=mean(ev), col="red") | ||
| + | legend("topright", "Average eigenvalue", lwd=1, col=2, bty="n") | ||
| - | 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. | + | # Résultats |
| + | summary(mite.spe.ca, display=NULL) | ||
| - | 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. | + | # Biplot |
| - | | + | windows() |
| + | plot(mite.spe.ca, scaling=1, type="none", | ||
| + | xlab=c("PC1 (%)", round((mite.spe.ca$CA$eig[1]/sum(mite.spe.ca$CA$eig))*100,2)), | ||
| + | ylab=c("PC2 (%)", round((mite.spe.ca$CA$eig[2]/sum(mite.spe.ca$CA$eig))*100,2))) | ||
| + | points(scores(mite.spe.ca, display="sites", choices=c(1,2), scaling=1), | ||
| + | pch=21, col="black", bg="steelblue", cex=1.2) | ||
| + | text(scores(mite.spe.ca, display="species", choices=c(1), scaling=1), | ||
| + | scores(mite.spe.ca, display="species", choices=c(2), scaling=1), | ||
| + | labels=rownames(scores(mite.spe.ca, display="species", scaling=1)), | ||
| + | col="red", cex=0.8) | ||
| + | </code> | ||
| + | |||
| + | Et votre biplot devrait ressembler à celui-ci: | ||
| + | |||
| + | {{ :ca_mite_biplot.png |}} | ||
| + | |||
| + | </hidden> | ||
| - | <code rsplus | Charger les données spatiales et classifier les abondances de poissons par site en fonction de la latitude> | + | =====3.3 Analyse en coordonnées principales (Principal Coordinates Analysis, PCoA)===== |
| - | #Charger les données spatiales pour déterminer les groupes | + | La PCA, comme la CA, impose une préservation des distances entre objets: la distance euclidienne dans le cas de la PCA, et la distance de Chi2 dans la CA. Si l'objectif est d'ordonner les objets sur la base d'une autre mesure de distance plus appropriée au problème, la PCoA constitue une technique de choix. Dans une PCA, les données sont pivotées de façon à ce que la première composante principale (correspondant à une combinaison linéaire des descripteurs) explique la plus forte proportion de variation possible; la contribution de chaque descripteur (espèces ou variables environnementales) à chaque composante principale peut alors être évaluée d'après son score. |
| - | spa <- read.csv ('http://www.davidzeleny.net/anadat-r/data-download/DoubsSpa.csv', row.names = 1) | + | La PCoA est une seconde méthode d'ordination sans contrainte dans laquelle les points sont ajoutés les uns après les autres à l'espace d'ordination en utilisant la distance euclidienne //ou n'importe quelle mesure de distance (dissimilarité) métrique vous choisissez//. Un premier point est ainsi placé dans l'espèce d'ordination, puis un second point placé à la valeur de distance du premier, puis un troisième et ainsi de suite en ajoutant autant d'axes (de dimensions) que nécessaire. |
| - | spa <- spa[,-8] | + | Il est parfois difficile de choisir entre effectuer une PCA ou une PCoA. La PCA permet toutefois de réduire des données multivariables en un faible nombre de dimensions tandis que la PCoA est utile pour visualiser les distances entre sites (ou objets). La PCoA est aussi particulièrement adaptées pour les jeux de données présentant plus de colonnes que de lignes. Par exemple, si des centaines d'espèces ont été observées dans un petit nombree de quadrats, une approche basée sur une PCoA utilisant la distance de Bray-Curtis (voir ci-dessous) peut être plus adaptée. |
| - | #Visualiser la matrice de données | + | PCoA avec DoubsSpe (transformé Hellinger): |
| - | View (spa) | + | <code rsplus | PCoA> |
| + | # En utilisant la fonction cmdscale() | ||
| + | ?cmdscale | ||
| + | cmdscale(dist(spe.hel), k=(nrow(spe)-1), eig=TRUE) | ||
| - | #Ajouter les numéros de site | + | # En utilisant la fonction pcoa() |
| - | numbers<-(1:30) | + | ?pcoa |
| - | numbers<-numbers[!numbers%in%8] | + | spe.h.pcoa<-pcoa(dist(spe.hel)) |
| - | spa$site<-numbers | + | |
| - | #Faire des groupes en fonction de la latitudey<82=group1, 82<y<156=group2, y>156=group3 | + | # Extraction des résultats |
| - | spa.group<-ddply(.data=spa, .variables=.(x, y, site), .fun= summarise, group = if(y <= 82) 1 else if (y <= 156) 2 else 3) | + | spe.h.pcoa |
| - | #Ordonner par site | + | # Représentation graphique |
| - | spa.group<-spa.group[with(spa.group, order(site)), ] | + | biplot.pcoa(spe .h.pcoa, spe.hel, dir.axis2=-1) |
| </code> | </code> | ||
| + | Les résultats de cette PCoA sont: | ||
| - | 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. | + | {{:pcoa_outputfr_1.png?800|}} |
| - | 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. | + | {{:pcoa_outputfr_2.png?800|}} |
| + | Et le graphique: | ||
| - | <code rsplus | Faire l'analyse discriminante linéaire (LDA)> | + | {{:pcoa_spe.png?500|}} |
| - | #faire la LDA | + | |
| - | LDA<-lda(env,spa.group[,4]) | + | |
| - | #classification des objets en fonction de la LDA | + | Vous pouvez aussi exécuter cette PCoA avec une autre mesure de distance (ex. Bray-Curtis): |
| - | spe.class <- predict(LDA)$class | + | |
| - | #probabilités que les objets appartiennent à chaque groupe a posteriori | + | <code rsplus | PCoA avec Bray-Curtis> |
| - | spe.post <- predict(LDA)$posterior | + | spe.bray.pcoa<-pcoa(spe.db) # il s'agit de la matrice de distances de Bray-Curtis qu'on a générée plus tôt |
| + | spe.bray.pcoa | ||
| + | biplot.pcoa(spe.bray.pcoa, spe.hel, dir.axis2=-1) | ||
| + | # Le choix d'une mesure de distance est très important car ça influence les résultats! | ||
| + | </code> | ||
| - | #tableau des classifications a priori et prédites | + | **Défi 5** |
| - | spe.table <- table(spa.group[,4], spe.class) | + | Exécuter une PCoA sur les données d'abondance des //espèces d'acariens// transformées Hellinger (données mite). Quels sont les axes importants? Quels groupes de sites pouvez-vous identifier? Quelles espèces sont liées à chaque groupe de sites? Comment les résultats de cette PCoA se comparent-ils avec ceux de la PCA? |
| - | #proportion de classification correcte | + | **Défi 5 - Solution** |
| - | diag(prop.table(spe.table, 1)) | + | <hidden> |
| + | <code rsplus | PCoA avec les données d'acariens> | ||
| + | mite.spe.h.pcoa<-pcoa(dist(mite.spe.hel)) | ||
| + | mite.spe.h.pcoa | ||
| + | windows() | ||
| + | biplot.pcoa(mite.spe.h.pcoa, mite.spe.hel, dir.axis2=-1) | ||
| </code> | </code> | ||
| - | {{ :lda_spetable.png?300 |}} | + | Représentation graphique: |
| + | {{:pcoa_mite.png?500|}} | ||
| - | 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? | + | Les espèces 16 et 31 sont plus éloignées des autres espèces en termes de distance, et donc leur distribution entre les sites est très différente de celle des autres espèces d'acariens. Les sites dont les étiquettes se chevauchent sont de bons exemples de sites à forte similarité en termes de communautés d'acariens. |
| + | </hidden> | ||
| - | 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. | + | =====3.4. Positionnement multidimensionnel non-métrique (Nonmetric Multidimensional Scaling, NMDS)===== |
| - | Chargez le fichier classifyme.csv. Celui-ci contient des données fictives de 5 nouveaux sites. | + | Les méthodes d'ordination non contrainte présentées ci-dessus permettent d'organiser les objets (ex. les sites) caractérisés par des descripteurs (ex. les espèces) dans un espace comprenant l'ensemble des dimensions décrites par l’ellipsoïde représentant le nuage des points de données. En d'autres termes, la PCA, la CA et la PCoA calculent un grand nombre d'axes d'ordination (nombre proportionnel au nombre de descripteurs) représentant la variation des descripteurs entre sites et préservant les distances entre objets (distance euclidienne dans une PCA, distance de Chi2 dans une CA et distance définie par l'utilisateur dans une PCoA). L'utilisateur peut ensuite sélectionner les axes d'intérêt (généralement les deux premiers axes d'ordination) pour représenter les objets dans un biplot. Le biplot produit représente ainsi correctement les distances entre objets (ex. la similarité des sites), mais ne permet pas de représenter l'ensemble des dimensions de la variation dans l'espace d'ordinations (étant donnée que l'Axe 3, l'Axe 4,..., l'Axe n'apparaissent pas sur le biplot, mais contribuent tout de même à expliquer la variation entre objets). |
| + | Dans certains cas, la priorité n'est pas de préserver la distance exacte entre les objets, mais au contraire de représenter aussi fidèlement que possible les relations entre objets selon un petit nombre d'axes (généralement deux ou trois) spécifiés par l'utilisateur. Dans de tels cas, le positionnement multidimensionnel non-métrique (NMDS) est la solution. Si l'utilisateur définit un nombre d'axe égal à deux, le biplot produit par le NMDS correspond à la meilleure solution graphique pour représenter en deux dimensions la similarité entre objets (les objets dissimilaires étant les plus éloignées, et les objets similaires étant les plus proches). De plus, le NMDS permet à l'utilisateur de choisir la mesure de distance qu'il souhaite pour ordonner les objets. | ||
| - | <code rsplus | Charger les données fictives et prédire le groupement des nouvelles données> | + | Afin de trouver la meilleure représentation des objets, le NMDS applique une procédure itérative qui vise à positionner les objets dans le nombre spécifié de dimensions de façon à minimiser une fonction de stress (variant de 0 à 1) qui mesure la qualité de l'ajustement de la distance entre objets dans l'espace d'ordination. Ainsi, plus la valeur du stress sera faible, plus la représentation des objets dans l'espace d'ordination sera exacte. Un second moyen d'évaluer l'exactitude d'un NMDS consiste à construire un diagramme de Shepard qui représente les distances entre objets sur le biplot d'ordination en fonction de leurs distances réelles. Le R2 obtenu à partir de la régression entre ces deux types de distance mesure la qualité de l'ajustement du NMDS. |
| - | #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 | + | <code rsplus | NMDS sur les données d’abondance spe avec une distance de Bray-Curtis et k=2 axes d'ordination> |
| - | group.new<-predict.group$class | + | # NMDS |
| - | </code> | + | spe.nmds<-metaMDS(spe, distance='bray', k=2) |
| + | |||
| + | ### Extraction des résultats | ||
| + | spe.nmds | ||
| - | {{ :lda_newgroups.png?200 |}} | + | ### Évaluation de la qualité de l'ajustement et construction du diagramme de Shepard |
| + | spe.nmds$stress | ||
| + | stressplot(spe.nmds, main='Shepard plot') | ||
| - | Nos nouveaux sites, en ordre, ont été classifiés dans les groupes 1, 1, 1, 3 et 3 respectivement. | + | # Construction du biplot |
| + | windows() | ||
| + | plot(spe.nmds, type="none", main=paste('NMDS/Bray - Stress=', round(spe.nmds$stress, 3)), | ||
| + | xlab=c("NMDS1"), | ||
| + | ylab=c("NMDS2")) | ||
| + | points(scores(spe.nmds, display="sites", choices=c(1,2)), | ||
| + | pch=21, col="black", bg="steelblue", cex=1.2) | ||
| + | text(scores(spe.nmds, display="species", choices=c(1)), | ||
| + | scores(spe.nmds, display="species", choices=c(2)), | ||
| + | labels=rownames(scores(spe.nmds, display="species")), | ||
| + | col="red", cex=0.8) | ||
| + | </code> | ||
| + | {{ :shepard_plot.png |}} | ||
| - | **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? | + | Le diagramme de Shepard identifie une forte corrélation entre les distances observées et les distances de l'ordination (R2 > 0.95), et donc une bonne qualité de l'ajustement du NMDS. |
| - | **Défi 5**: Solution | + | {{ :nmds_biplot.png |}} |
| - | <hidden> | + | Le biplot du NMDS identifie un groupe de sites caractérisés par les espèces BLA, TRU, VAI, LOC, CHA et OMB, tandis que les autres espèces caractérisent un groupe de sites situés dans le coin supérieur droit du biplot. Quatre sites situés dans le coin inférieur droit sont fortement différents des autres. |
| - | <code rsplus | LDA on mite data> | + | |
| - | mite.xy$site<-seq(1:70) | + | |
| - | (max(mite.xy[,2])-min(mite.xy[,2]))/4 | + | |
| - | mite.xy.group<-ddply(.data=mite.xy, .variables=.(x, y, site), .fun= summarise, group = if(y <= 2.5) 1 else if (y <= 4.9) 2 else if (y <= 7.3) 3 else 4) | ||
| - | mite.xy.group<-mite.xy.group[with(mite.xy.group, order(site)), ] | ||
| - | LDA.mite<-lda(mite.env[,1:2],mite.xy.group[,4]) | + | **Défi 6** |
| - | mite.class <- predict(LDA.mite)$class | + | Exécuter un NMDS sur les données d'abondance des //espèces d'acariens// (données mite) en deux dimensions à partir de distances de Bray-Curtis. Évaluer la qualité de l'ajustement et interpréter le biplot. |
| - | mite.post <- predict(LDA.mite)$posterior | + | |
| - | mite.table <- table(mite.xy.group[,4], mite.class) | + | |
| - | diag(prop.table(mite.table, 1)) | + | |
| - | </code> | + | |
| - | </hidden> | + | **Défi 5 - Solution** |
| + | <hidden> | ||
| + | <code rsplus> | ||
| + | ### NMDS | ||
| + | mite.spe.nmds<-metaMDS(mite.spe, distance='bray', k=2) | ||
| + | ### Extraction des résultats | ||
| + | mite.spe.nmds | ||
| - | ======5. Autres méthodes d'ordination utiles====== | + | ### Évaluation de la qualité de l'ajustement |
| + | mite.spe.nmds$stress | ||
| + | stressplot(mite.spe.nmds, main='Shepard plot') | ||
| - | <code rsplus | autres méthodes> | + | ### Construction du biplot |
| - | ?cca # Analyse canonique des correspondances (Constrained Correspondence Analysis, CCA) | + | windows() |
| - | + | plot(mite.spe.nmds, type="none", main=paste('NMDS/Bray - Stress=', round(mite.spe.nmds$stress, 3)), | |
| - | # Méthode d’analyse canonique similaire à la RDA préservant les distance de Chi-carré entre objets (au lieu des | + | xlab=c("NMDS1"), |
| - | # distances euclidiennes dans le cas d’une RDA). Cette méthode est mieux adaptée à l’étude de longs gradients | + | ylab=c("NMDS2")) |
| - | # que la RDA. | + | points(scores(mite.spe.nmds, display="sites", choices=c(1,2)), |
| - | + | pch=21, col="black", bg="steelblue", cex=1.2) | |
| - | + | text(scores(mite.spe.nmds, display="species", choices=c(1)), | |
| - | + | scores(mite.spe.nmds, display="species", choices=c(2)), | |
| - | ?CCorA # Analyse canonique de corrélations (Canonical Correlation Analysis, CCorA) | + | labels=rownames(scores(mite.spe.nmds, display="species")), |
| - | + | col="red", cex=0.8) | |
| - | # 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 | + | |
| </code> | </code> | ||
| + | {{ :nmds_mite_shepard.png |}} | ||
| - | ======Références====== | + | La corrélation entre distance observée et distance d'ordination (R2 > 0.91) et la valeur de stress relativement faible identifient une bonne qualité de l'ajustement du NMDS. |
| - | 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. | + | {{ :nmds_mite_biplot.png |}} |
| - | Borcard, Gillet & Legendre (2011). Numerical Ecology with R. Springer New York. | + | Aucun groupe de sites ne peut être précisément identifié à partir du biplot, ce qui montre que la plupart des espèces sont présentes dans la plupart des sites, i.e. peu de sites présentent des communautés distinctes. |
| - | + | </hidden> | |
| - | 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. | + | =====3.5. En résumé===== |
| - | Gotelli & Ellison (2004). A Primer of Ecological Statistics. Sinaeuer Associates Inc., Sunderland MA. | + | L’ordination constitue une puissante méthode d'analyse pour étudier les relations entre objets caractérisés par différents descripteurs (ex. des sites décrits par leurs communautés biologiques, ou leurs variables environnementales), mais de nombreuse méthodes d'ordination existent. Ces méthodes diffèrent principalement par le type de distance qu'elles préservent, le type de variables qu'elles autorisent, et le nombre de dimensions de l'espace d'ordination. Pour mieux guider votre choix de la méthode d'ordination à utiliser, le tableau ci-dessous identifie les caractéristiques de chacune des quatre méthodes d'ordination présentées lors de cet atelier. |
| - | Legendre & Legendre (2012). Numerical Ecology 3rd edition. Elsevier Science BV, Amsterdam. | + | {{ :resume_ordination.jpg |}} |
| - | Poulin, Andersen & Rochefort (2013) A new approach for tracking vegetation change after restoration: a case study with peatlands. Restoration Ecology, 21, 363-371. | + | Lors du prochaine atelier, vous verrez comment identifier les relations entre variables environnementales et communautés biologiques décrivant un même ensemble de sites, à l'aide des méthodes d'analyses canoniques. |
