Differences

This shows you the differences between two versions of the page.

--- r_atelier10 [2015/09/08 14:45]
zofia.taranu created
+++ r_atelier10 [2021/10/14 03:52] (current)
lsherin
@@ Line 1: / Line 1: @@
+<WRAP group>
+<WRAP centeralign>
+<WRAP important>
+<wrap em> __AVIS IMPORTANT__ </wrap>
+<wrap em> Depuis l'automne 2021, ce wiki a été discontinué et n'est plus activement développé. </wrap>
+<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-10/|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>
+<wrap em> Merci de votre compréhension, </wrap>
+<wrap em> Vos coordonnateurs de la série d’ateliers R du CSBQ. </wrap>
+</WRAP>
+</WRAP>
+<WRAP clear></WRAP>
 ======= Série d'ateliers en R du CSBQ =======
@@ Line 5: / Line 22: @@
 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.
-====== Atelier 10: Programmation en R ======
+//Le contenu de cet atelier a été révisé par plusieurs membres du CSBQ. Si vous souhaitez y apporter des modifications, veuillez SVP contacter les coordonnateurs actuels de la série, listés sur la page d'accueil//
-Développé par: Johanna Bradie, Sylvain Christin, Ben Haller, Guillaume Larocque
-**Résumé:** Cet atelier vise à vous apprendre les bases de la programmation en R. Vous apprendrez à utiliser des structures de contrôle (boucles for, if, while) afin d'éviter la répétition de code, de faciliter l'organisation et d'effectuer des simulations. Vous apprendrez également à écrire vos propres fonctions et quelques astuces pour programmer plus efficacement. La dernière partie de l'atelier portera sur des librairies de R qui peuvent être très utiles pour les participants, mais qui n'ont pas été couvertes dans la série d'ateliers en R du CSBQ.
+<wrap em>AVIS IMPORTANT: MISES À JOUR MAJEURES</wrap>
-Lien vers le Prezi associé: [[https://prezi.com/vjbz_n-rweof/csbq-atelier-r-8/|Prezi]]
+**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/workshop10|la page GitHub]] des ateliers R du CSBQ.
-Téléchargez le script R pour cet atelier:
+Le matériel disponible inclut;
-  * [[http://qcbs.ca/wiki/_media/qcbs_workshop8_r_code.r | script R]]
+  - La [[https://qcbsrworkshops.github.io/workshop10/pres-fr/workshop10-pres-fr.html|présentation Rmarkdown]] pour cet atelier;
+  - Un [[https://qcbsrworkshops.github.io/workshop10/book-fr/workshop10-script-fr.R|script R]] qui suit la présentation - //*en construction*//.
+  - [[https://qcbsrworkshops.github.io/workshop10/book-fr/index.html|Le matériel écrit]] qui accompagne la présentation en format bookdown - //*en construction*//.
-===== Objectifs d'apprentissage =====
-  - Structures de contrôle
+====== Atelier 10: Analyses multivariées avancées ======
-  - Écriture de fonctions en R
-  - Réduire le temps d’exécution des codes
-  - Paquets utilises pour les biologistes
-===== 1. Structures de contrôle =====
+Développé par: Monica Granados, Emmanuelle Chrétien, Bérenger Bourgeois, Amanda Winegardner and Xavier Giroux-Bougard. (Le matériel des scripts R est adapté de: Borcard, Gillet & Legendre (2011). //Numerical Ecology with R//. Springer New York.)
-Les structures de contrôle vous permettent d'exécuter un même groupe de commandes plusieurs fois selon des conditions diverses. Dans cette section, vous apprendrez comment:
-  * Exécuter des commandes **conditionnellement** en utilisant: if, if/else
-  * Exécuter les mêmes commandes **plusieurs fois** en utilisant: des boucles, for, while, repeat
-  * Modifier l'exécution des boucles en utilisant: break, next.
-==== Commandes if et if/else ====
+**Résumé:** Durant cet atelier, vous apprendrez à réaliser des analyses multivariées avancées sur des données de communauté. Cet atelier se concentre sur les méthodes sous contraintes, telles que l'analyse canonique de redondances (RDA), l'arbre de régression multivarié (MRT) et l'analyse discriminante linéaire (LDA) afin
+d'explorer comment les variables environnementales peuvent expliquer les patrons de composition en espèces à travers différents sites.
-Les commandes if et if/else sont utiles pour:
-  * Vérifier s'il y a des problèmes ou le non-respect de conditions.
-  * Traiter différentes lignes de vos données de façons différentes.
-  * Vérifier l'existence d'un fichier ou d'une variable.
-**Syntaxe**
+**Lien vers la nouvelle [[https://qcbsrworkshops.github.io/workshop10/workshop10-fr/workshop10-fr.html|présentation Rmarkdown]]**
-<code rsplus>
+Lien vers la présentation Prezi associée : [[https://prezi.com/iuaffmpxhie7/
-if (condition) {
+| Prezi]]
-  expression      # N'importe quelle commande que vous voulez exécuter.
-}
-if (condition) {
-  expression
-} else {
-  expression
-}
-</code>
-Dans ce qui suit, nous allons utiliser les **opérateurs logiques**.
+Téléchargez le code R, les librairies et données requises pour cet atelier:
-Notez que pour tester si deux valeurs sont égales, vous devez utiliser "==".
+  *  [[http://qcbs.ca/wiki/_media/multivar2_f.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]]
-| == | égal à |
-| != | pas égal à |
-| !x | non x |
-| < | plus petit que |
-| < = | plus petit que ou égal à |
-| > | plus grand que |
-| >= | plus grand que ou égal à |
-| x & y | x ET y |
-| x%%|%%y| x OU y |
-| isTRUE(x) | est-ce que X est vrai? |
+Assurez-vous d'importer les librairies suivantes dans R Studio (procédure fournie dans le code R):
+  *  [[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
-Par exemple,
-<code rsplus>
-if ((2+2) == 4) {
-  print("Les maths, c'est logique!.")
-}
-if ((2+1) == 4) {
+<code rsplus | Installer et importer les librairies>
-  print("Les maths, c'est logique!.")
+install.packages("vegan")
-}
+install.packages("mvpart")
-</code>
+install.packages("labdsv")
+install.packages("plyr")
+install.packages("MASS")
-Les accolades { } sont utilisées pour indiquer à R que d'autre chose s'en vient. Quand on utilise une accolade, R attend l'accolade de fermeture pour exécuter la commande en entier. Si les accolades ne sont pas utilisées, R peut ne pas se comporter tel qu'attendu. Par exemple, essayez:
+# Pour les deux librairies suivantes, importez le fichier d'archivage donné sur la page wiki.
+# Pour ce faire, cliquez sur l'onglet "Packages" de la section en bas à droite sur R Studio
+# Cliquez sur "Install Packages"
+# Choisissez d'installer depuis "Package Archive file" et importez les deux fichiers d'archivage
+install.packages("MVPARTwrap")
+install.packages("rdaTest")
-<code rsplus>
+library(vegan)
-if ((2 + 1) == 4) print("Les maths, c'est logique!.")
+library(mvpart)
-else print("Houston, on a un problème.")
+library(MVPARTwrap)
+library(rdaTest)
+library(labdsv)
+library(plyr)
+library(MASS)
 </code>
-La commande else ne fonctionne pas puisque R évalue la première ligne sans savoir que votre commande n'est pas complète.
-À la place, utilisez:
-<code rsplus>
-if ((2 + 2) == 4) {
-  print("Les maths, c'est logique!.")  # R n'évalue pas encore cette expression puisque l'accolade n'est pas fermée.
-} else {
-  print("Houston, on a problem.")
-}  # Comme toutes les accolades sont fermées, R va évaluer les commandes en entier.
-</code>
-Notez que if et if/else testent une condition sur une valeur unique. Si vous voulez tester un vecteur de conditions, et obtenir un vecteur en résultat, utilisez la fonction ifelse.
+======Introduction aux analyses multivariées avancées======
-Par exemple,
+L'atelier précédent donnait un aperçu des analyses multivariées de base:
-<code rsplus>
+    * comment choisir les mesures de distances et les transformations appropriées selon le type de données
-a <- 1:10
+    * groupement hiérarchique
-ifelse(a > 5, "oui", "non")
+    * ordinations sans contraintes (analyses en composantes principales, analyse en coordonnées principales, analyse de correspondances, positionnement multidimensionnel non-métrique)
-</code>
-Vous pouvez aussi utiliser la commande ifelse dans une autre fonction afin d'effectuer une opération seulement selon certaines conditions:
+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.
-Par exemple,
+======Introduction aux données======
-<code rsplus>
-a <- (-4):5
-sqrt(ifelse(a >= 0, a, NA))
-</code>
+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 les données>
-**Exercice 1**
+#Matrice d'abondances d'espèces: “DoubsSpe.csv”
+spe<- read.csv(file.choose(), row.names=1)
+spe<- spe[-8,] #Pas d'espèces dans le site 8, supprimer site 8. Exécuter cette ligne une seule fois
-<code rsplus>
+#Matrice de données environnementales: “DoubsEnv.csv”
-Minou <- "chat"
+env<- read.csv(file.choose(), row.names=1)
-Pitou <- "chien"
+env<- env[-8,] #Supprimer le site 8 puisqu'on l'a supprimé de la matrice d'abondance. N'exécuter qu'une seule fois.
-Filou <- "chat"
-animaux <- c(Minou, Pitou, Filou)
 </code>
-. Utilisez une commande "if" pour afficher "meow" si Animal a la valeur "chat".
-++++ Exercice 1.1 : Réponse|
-<code rsplus>
-if(Minou == "chat") {
-  print("meow")
-}
-</code>
-++++
-. Utilisez une commande if/else pour afficher "woof" si Animal a la valeur "chien" et "meow" sinon. Essayez en d'abord la valeur Animal="chien" et ensuite "lion".
+======1. Exploration et préparation des données======
+=====1.1 Données d'abondances d'espèces =====
-++++ Exercice 1.2 : Réponse|
+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>
-if(Animal == "chien") {
-  print("woof")
-} else {
-  print("meow")
-}
-</code>
-++++
-. Utilisez la commande ifelse pour afficher "woof" si les animaux sont des chiens et "meow" pour les chats.
-++++ Exercice 1.3 : Réponse|
+<code rsplus | Explorez DoubsSpe>
-<code rsplus>
+names(spe) #voir les noms des colonnes spe
-ifelse(animaux == "chien", "woof", "meow")
+dim(spe) #dimensions de spe; nombre de lignes et de colonnes
+str(spe) #la structure interne de la matrice
+head(spe) #les premières lignes
+summary(spe) #statistiques descriptives
 </code>
-++++
------
-==== Boucles ====
-Les boucles sont utiles pour:
+Regardez la distribution des espèces.
-  * faire quelque chose pour chaque élément d'un objet.
-  * faire quelque chose jusqu'à la fin du traitement de données.
-  * faire quelque chose pour chaque fichier dans un répertoire.
-  * faire quelque chose qui peut échouer, jusqu'à ce que ça marche.
-  * faire des calculs itératifs jusqu'à convergence.
-===boucles for===
+<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
-La boucle **for** est le type de boucle la plus commune. Utilisez là pour exécuter un bloc de commande un nombre fixe de fois.
+barplot(ab, las=1, xlab="Abundance class", ylab="Frequency", col=grey(5:0/5))
-**Syntaxe**
-<code rsplus>
-for (variable in séquence) {
-  expression
-}
 </code>
-L'exécution une fois d'un groupe de commandes dans un boucle s'appelle une itération.
+{{:spe_barplot.png?300|}}
+Il y a une grande fréquence de zéros dans les données d'abondance.
-Par exemple:
+Calculez le nombre d'absences.
-<code rsplus>
+<code rsplus | Absences dans les données de poissons>
-for (i in 1:5) {
+sum(spe==0)
-  print(i)
-}
 </code>
-Dans cet exemple, notre séquence contient 5 éléments (1, 2, 3, 4, 5). R va ensuite évaluer l'expression 5 fois. La variable i prendra successivement les valeurs de 1 à 5 pendant l'exécution de la boucle. Pendant l'itération 1, R remplacera chaque valeur de i dans la boucle par 1, ainsi de suite.
+Regardez la proportion de zéros dans les données d'abondances de poissons.
+<code rsplus | Proportion de zéros>
+sum(spe==0)/(nrow(spe)*ncol(spe))
+</code>
-La lettre "i" peut être remplacée par n'importe quelle nom de variable et la séquence peut être à peut prêt n'importe quoi, même une liste de vecteurs.
+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().
-Essayez:
+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).
-<code rsplus>
-for (m in 4:10) {
-  print(m * 2)
-}
-for (a in c("Bonjour", "programmeurs", "en R")) {
-  print(a)
-}
-for (z in 1:30) {
-  a <- rnorm(n = 1, mean = 5, sd = 2) # obtenir une valeur aléatoire provenant d'une distribution normale avec une moyenne de 5 et un écart type de 2.
-  print(a)
-}
-elements <- list(1:3, 4:10)
+<code rsplus | transformation de Hellinger>
-for (element in elements) {
+spe.hel <- decostand(spe, method="hellinger") # vous pouvez également utiliser method="hell"
-  print(element)
-}
 </code>
-Les boucles sont souvent utilisées pour faire des opérations successives sur un jeu de données. Nous utiliserons les boucles pour évaluer des fonctions sur le jeu de données CO2 qui vient avec R. Ce jeu de données contient des concentrations et valeurs d'absorption par des plantes situées au Québec et au Mississippi qui ont été exposées à des traitements de refroidissement ("chilled") ou sans refroidissement (non-chilled). Notez que c'est le même jeu de données utilisé pour l'atelier 2. Le code ci-dessous fournit quelques exemples d'utilisation de boucles avec ces données.
-<code rsplus>
+=====1.2 Données environnementales=====
-data(CO2)
+Explorez les données environnementales et détectez les colinéarités.
-for (i in 1:length(CO2[,1])) { # pour chaque ligne du jeu de donnée CO2
-  print(CO2$conc[i]) #affichez les concentrations de CO2
-}
-for (i in 1:length(CO2[,1])) { # pour chaque ligne du jeu de donnée CO2
+<code rsplus | Exploration de la matrice Doubs Env>
-  if(CO2$Type[i] == "Quebec") { # si le type est Quebec
+names(env)
-    print(CO2$conc[i]) #affichez les concentrations de CO2 }
+dim(env)
-  }
+str(env)
-}
+head(env)
+summary(env)
-# Truc 1 : pour obtenir le nombre de lignes d'un jeu de donnée, on peut utiliser la fonction nrow().
+pairs(env, main="Bivariate Plots of the Environmental Data" )
-for (i in 1:nrow(CO2)) { # pour chaque ligne du jeu de donnée CO2
-  print(CO2$conc[i]) # affichez les concentrations de CO2
-}
-# Truc 2 : Si on veut faire des opérations seulement sur les valeurs d'une seule colonne, on peut s'en servir comme source pour l'itération.
-for (i in CO2$conc) { # pour chacune des valeurs de concentration de CO2
-  print(i) # afficher cette valeur
-}
 </code>
-La partie "expression" de la boucle peut contenir plusieurs lignes de commandes différents.
+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>
+<code rsplus | Decostand>
-for (i in 4:5) { # pour i de 4 à 5
+env.z <- decostand(env, method="standardize")
-  print(colnames(CO2)[i])
+apply(env.z, 2, mean) # les données sont centrées (moyenne~0)
-  print(mean(CO2[,i])) # afficher les moyennes de cette colonne
+apply(env.z, 2, sd)   # les données sont réduites (écart type =1)
-}
 </code>
-Notez que certaines de ces opérations peuvent se faire plus efficacement en utilisant la fonction apply(), mais ce n'est pas vrai pour les commandes plus complexes comprenant plusieurs lignes d'expressions.
-=== boucles while et repeat ===
-les //boucles while// et //repeat// opèrent de façon similaire aux //boucles for//. Une fois que vous avez compris le fonctionnement des boucles "for", vous devriez comprendre facilement les autres sortes de boucles. Dans la prochaine section, vous verrez des exemples d'utilisation des boucles while et repeat. Notez que dans plusieurs situations, vous pouvez accomplir la même tâche de plusieurs façons différentes, soit en utilisant boucle for, while ou repeat, tel que démontré ci-dessous.
+======2. Analyses canoniques======
-=== boucles nichées ===
+Décrites à l’origine par Rao (1964), les analyses canoniques rassemblent de nombreuses méthodes statistiques combinant les concepts d’ordination et de régression et partageant un but commun, à savoir identifier les relations entre un ensemble de variables réponse (matrice Y décrivant généralement la composition en espèces de différentes communautés) et un ensemble de variables explicatives (matrice X contenant le plus souvent des variables environnementales). Les analyses canoniques permettent de tester des hypothèses écologiques relatives aux déterminants environnementaux de la composition en espèces des communautés. Parmi l’ensemble de méthodes d’analyses canoniques existantes, nous insisterons ici principalement sur l’analyse canonique de redondance (RDA).
+=====2.1 Analyse canonique de redondances (ACR ou RDA)=====
-Dans certains cas, vous voudrez peut-être utiliser des boucles nichées une dans l'autre. Dans ce cas, il est important d'utiliser un nom de variable d'itération différent pour chaque boucle (ici on utilise i et n).
+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.
-<code rsplus>
+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 :
-for (i in 1:5) {
-  for (n in 1:5) {
+Yfit = X[X'X]-1 [X'Y]
-    print (i*n)
-  }
-}
-</code>
------
+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.
-**Exercice 2**
-. Vous avez pris conscience que votre outil pour mesurer le l'absorption de CO2 n'était pas bien calibré aux sites situés au Québec, et toutes les mesures sont deux unités trop élevées. Utilisez une boucle pour corriger les mesures pour tous les sites aux Québec.
+{{ :constrained_ord_diagram.png |}}
-++++ Exercice 2 : Réponse|
-<code rsplus>
-for (i in 1:length(CO2[,1])) {
-  if(CO2$Type[i] == "Quebec") {
-    CO2$uptake[i] <- CO2$uptake[i] - 2
-  }
-}
-</code>
+Une RDA est un processus en deux étapes (d’après Legendre et Legendre, 2012)
-Faites attention de ne pas rouler ce code plus d'une fois puisque que les valeurs vont graduellement diminuer.
-++++
------
-Assurez-vous de bien recharger le jeu de données pour ainsi travailler avec les données originales pour le reste de l'exercice:
-<file rsplus>
-data(CO2)
-</file>
-\\
-==== Modifications aux boucles ====
-Normalement, les itérations s'exécutent successivement jusqu'à la dernière. Pour changer ce comportement, on peut utiliser **break** pour rompre complètement l'exécution de la boucle, ou **next** pour arrêter l'exécution de l'itération courante et passer à la suivante.
+Plusieurs valeurs statistiques peuvent être extraites de la RDA, en particulier :
-Par exemple,
+- Le R² mesure la force de la relation canonique entre Y et X en calculant la proportion de la variation de Y expliquée par les variables X,
+- Le R² ajusté mesure également la force de la relation entre Y et X, mais applique une correction du R² afin de tenir compte du nombre de variables explicatives,
+- La statistique F correspond à un test global de la significativité de la RDA en comparant le modèle étudié à un modèle nul.
-<file rsplus>
+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.
-# Afficher les concentrations de CO2 pour les traitements "chilled" et garder le compte du nombre total de traitements.
+Avec R, la RDA peut être calculée en utilisant la fonction rda du package vegan, comme suit:
-count <- 0 # la valeur de count est mise à zéro de façon à ensuite pouvoir modifier la valeur dans la boucle.
-for (i in 1:length(CO2[,1])) {
+<code rsplus | effectuer une RDA avec rda() de vegan>
-  if (CO2$Treatment[i] == "nonchilled") next #Passer à l'itération suivante si c'est "nonchilled"
+#Préparer les données
-  count <- count + 1
+env.z <- subset(env.z, select = -das) # Enlever la variable "distance from the source"
-  print(CO2$conc[i])
-}
+#Faire la RDA
-print(count) # La valeur de count est maintenant égale à 42.
+?rda
+spe.rda <- rda(spe.hel~., data=env.z)
-# Ceci pourrait être écrit de façon équivalente en utilisant une boucle repeat:
+### Extraire les résultats
+summary(spe.rda, display=NULL)
+#Les résultats peuvent être vus en utilisant summary():
+summary(spe.rda, display=NULL) #display = NULL optional
+</code>
-count <- 0
-i <- 0
-repeat {
-      i <- i + 1
-      if (CO2$Treatment[i] == "nonchilled") next  # sauter cette itération
-      count <- count + 1
-      print(CO2$conc[i])
-      if (i == length(CO2[,1])) break     # arrêter l'itération
-    }
-print(count)
+Le sommaire ressemble à ceci:
-### On pourrait aussi l'écrire avec une boucle while
-i <- 0
+{{ :rda_output_1.png |}}
-count <- 0
-while (i < length(CO2[,1]))
-{
-  i <- i + 1
-  if (CO2$Treatment[i] == "nonchilled") next  # sauter cette itération
-  count <- count + 1
-  print(CO2$conc[i])
-}
-print(count)
-### On pourrait aussi l'écrire plus simplement avec un opérateur logique servant d'index:
+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.
-print(CO2$conc[CO2$Treatment=="chilled"])
+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):
-print(sum(CO2$Treatment=="chilled"))
-</file>
------
+<code rsplus | ordiR2step() pour la sélection progressive>
-**Exercice 3**
+?ordiR2step
+ordiR2step(rda(spe.hel~1, data=env.z), scope= formula(spe.rda), direction= "forward", R2scope=TRUE, pstep=1000)
+env.signif <- subset(env.z, select = c("alt", "oxy", "dbo"))
+</code>
-. Vous venez de vous rendre compte que votre outil pour mesurer la concentration ne fonctionne pas correctement. Sur les sites au Mississippi, les concentrations de moins de 300 était correctes, mais les concentrations de plus de 300 étaient surestimées par 20 unités. Utilisez une boucle pour corriger ces mesures.
-++++ Exercice 3 : Réponse|
+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 :
-<file rsplus>
-for (i in 1:length(CO2[,1])) {
-  if(CO2$Type[i] == "Mississippi") {
-    if(CO2$conc[i] < 300) next
-    CO2$conc[i] <- CO2$conc[i] - 20
-  }
-}
-# Note : on pourrait également utiliser une seule boucle if, ce qui est plus claire
-for (i in 1:nrow(CO2)) {
-  if(CO2$Type[i] == "Mississippi" && CO2$conc[i] >= 300) {
-    CO2$conc[i] <- CO2$conc[i] - 20
-  }
-}
-# Ou utiliser un index logique
+<code rsplus | RDA avec seulement les variables significatives>
-ind=CO2$Type == "Mississippi" & CO2$conc >= 300 # cet index a des valeurs de "TRUE" quand les conditions sont vraies et FALSE sinon
+spe.rda.signif <- rda(spe.hel~., data=env.signif)
-CO2[ind,'conc']= CO2$conc[ind] - 20
+summary(spe.rda.signif, display=NULL)
+</code>
-</file>
-++++
------
-Assurez-vous de travailler avec les données originales pour le reste de l'exercice
+Les variables explicatives expliquent désormais 59% de la variance de Y.
-<file rsplus>
-data(CO2)
-</file>
+Le R² ajusté de cette RDA est calculé à partir de la fonction RsquareAdj :
-**Utiliser les structures de contrôle pour faire un graphique complexe**
-Nous avons l'idée de créer un graphique à partir des données de concentration et absorption où chaque point est associé à un type (Québec ou Mississippi) et un traitement ("chilled" et "nonchilled") et nous voulons représenter ces points différemment sur le graphique.
+<code rsplus | R2 ajusté>
+(R2adj <- RsquareAdj(spe.rda.signif)$adj.r.squared)
+</code>
-Vous pouvez en lire davantage sur le "typesetting" mathématique avec ?plotmath, et voir les options de couleurs, grosseurs, rotations, etc. avec ?par.
-<file rsplus>
+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):
-head(CO2) # Voir les données
-unique(CO2$Type)
-unique(CO2$Treatment)
-# Créer le graphique dans lequel chaque type et traitement a une couleur différente
+<code rsplus | anova.cca pour tester la significativité des axes>
+?anova.cca
+anova.cca(spe.rda.signif, step=1000)
+anova.cca(spe.rda.signif, step=1000, by="axis")
+#Le modèle global est significatif ainsi que les trois axes canoniques.
+</code>
-plot(x=CO2$conc, y=CO2$uptake, type="n", cex.lab=1.4, xlab="CO2 concentration", ylab="CO2 uptake") # Type "n" dit à R de ne pas créer le graphique
-for (i in 1:length(CO2[,1])) {
-  if (CO2$Type[i] == "Quebec" & CO2$Treatment[i] == "nonchilled") {
-    points(CO2$conc[i], CO2$uptake[i], col="red",type="p")
-  }
-  if (CO2$Type[i] == "Quebec" & CO2$Treatment[i] == "chilled") {
-    points(CO2$conc[i], CO2$uptake[i], col="blue")
-  }
-  if (CO2$Type[i] == "Mississippi" & CO2$Treatment[i] == "nonchilled") {
-    points(CO2$conc[i], CO2$uptake[i], col="orange")
-  }
-  if (CO2$Type[i] == "Mississippi" & CO2$Treatment[i] == "chilled") {
-    points(CO2$conc[i], CO2$uptake[i], col="green")
-  }
-}
-</file>
+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.
-Notez, qu'il y a d'autres façons de créer des graphiques complexes. ggplot peut, qui a été utiliser dans l'exercice 4 peut également être utilisé dans ce cas.
+Pour obtenir un triplot scaling 1 d’une RDA, le code suivant peut être employé:
------
-**Exercice 4**
-. Créez un graphique montrant les concentrations en fonction de l'absorption et où chaque plante a des points de différents couleurs. Essayez de le faire avec une boucle nichée!
+<code rsplus | représentations graphiques des RDAs>
+#Pour faire des graphiques rapidement et facilement
+windows()
+plot(spe.rda.signif, scaling=1, main="Triplot RDA (scaling 1)")
+windows()
+plot(spe.rda.signif, scaling=2, main="Triplot RDA (scaling 2)")
+#Code avancé pour graphiques plus attrayants
+#cadrage 1
+windows()
+plot(spe.rda.signif, scaling=1, main="Triplot RDA - scaling 1", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
+points(scores(spe.rda.signif, display="sites", choices=c(1,2), scaling=1),
+       pch=21, col="black", bg="steelblue", cex=1.2)
+arrows(0,0,
+       scores(spe.rda.signif, display="species", choices=c(1), scaling=1),
+       scores(spe.rda.signif, display="species", choices=c(2), scaling=1),
+       col="black",length=0)
+text(scores(spe.rda.signif, display="species", choices=c(1), scaling=1),
+     scores(spe.rda.signif, display="species", choices=c(2), scaling=1),
+     labels=rownames(scores(spe.rda.signif, display="species", scaling=1)),
+     col="black", cex=0.8)
+arrows(0,0,
+       scores(spe.rda.signif, display="bp", choices=c(1), scaling=1),
+       scores(spe.rda.signif, display="bp", choices=c(2), scaling=1),
+       col="red")
+text(scores(spe.rda.signif, display="bp", choices=c(1), scaling=1)+0.05,
+     scores(spe.rda.signif, display="bp", choices=c(2), scaling=1)+0.05,
+     labels=rownames(scores(spe.rda.signif, display="bp", choices=c(2), scaling=1)),
+     col="red", cex=1)
+#cadrage 2
+windows()
+plot(spe.rda.signif, scaling=2, main="Triplot RDA - scaling 2", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
+points(scores(spe.rda.signif, display="sites", choices=c(1,2), scaling=2),
+       pch=21, col="black", bg="steelblue", cex=1.2)
+arrows(0,0,
+       scores(spe.rda.signif, display="species", choices=c(1), scaling=2)*2,
+       scores(spe.rda.signif, display="species", choices=c(2), scaling=2)*2,
+       col="black",length=0)
+text(scores(spe.rda.signif, display="species", choices=c(1), scaling=2)*2.1,
+     scores(spe.rda.signif, display="species", choices=c(2), scaling=2)*2.1,
+     labels=rownames(scores(spe.rda.signif, display="species", scaling=2)),
+     col="black", cex=0.8)
+arrows(0,0,
+       scores(spe.rda.signif, display="bp", choices=c(1), scaling=2),
+       scores(spe.rda.signif, display="bp", choices=c(2), scaling=2),
+       col="red")
+text(scores(spe.rda.signif, display="bp", choices=c(1), scaling=2)+0.05,
+     scores(spe.rda.signif, display="bp", choices=c(2), scaling=2)+0.05,
+     labels=rownames(scores(spe.rda.signif, display="bp", choices=c(2), scaling=2)),
+     col="red", cex=1)
+</code>
-++++ Exercice 4 : Réponse|
-<file rsplus>
-plot(x=CO2$conc, y=CO2$uptake, type="n", cex.lab=1.4,xlab="CO2 concentration", ylab="CO2 uptake")
-plants <- unique(CO2$Plant)
+Les triplots finaux:
-for (i in 1:length(CO2[,1])){
+{{ :doubs_rda1.png |}}
-  for (p in 1:length(plants)) {
+{{ :doubs_rda2.png |}}
-    if (CO2$Plant[i] == plants[p]) {
-      points(CO2$conc[i], CO2$uptake[i], col=p, type="p")
-    }
-  }
-}
-</file>
-++++
------
-===== 2. Créer ses fonctions =====
-==== Pourquoi créer ses fonctions? ====
+**Défi 1**: Effectuer la RDA de l’abondance des espèces du tableau d'acariens en fonction des variables environnementales mite.env
-La plupart du travail lourd dans R est effectué par les fonctions. Elles sont utiles pour:
-  * répéter une même tâche mais en changeant ses paramètres
-  * rendre votre code plus lisible
-  * rendre votre code plus facile à modifier et à maintenir
-  * partager du code entre différentes analyses
-  * partager votre code avec d'autres personnes
-  * modifier les fonctionalités par défaut de R
+<code rsplus | Chargez les données d'acariens>
+#Ces données sont disponibles dans la librairie vegan
+data(mite)
+mite.spe<-mite
+mite.spe.hel <- decostand(mite.spe, method="hellinger")
-Mais qu'est ce qu'une fonction au juste? Une fonction, c'est essentiellement une boîte noire qui transforme des données. Elle prend en entrée des valeurs - appelées arguments -, utilise du code R pour les traiter et renvoie optionnellement une valeur de retour.
+data(mite.env)
+</code>
+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 ?
-{{::fonction_schema.png|200}}
-==== Comment écrire des fonctions? ====
+**Défi 1**: Solution
-Voici la syntaxe basique d'une fonction:
+<hidden>
+Votre code devrait ressembler à ceci:
-<code rsplus>
+<code rsplus | RDA sur les données des acariens>
-nom_de_la_fonction <- function(argument1, argument2, ...) {
+#RDA avec toutes les variables environnementales
-  expression...  # Ce que l'on veut que la fonction fasse
+mite.spe.rda<-rda(mite.spe.hel~., data=mite.env)
-  return(value)  # Optionnel. Si l'on veut acceder au resultat de la fonction
-}
+#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>
-=== Arguments ===
+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.
-Les arguments sont les données fournies en entrée à votre fonction. Il s'agit de l'information dont votre fonction a besoin pour marcher correctement. Une fonction peut avoir entre 0 et une infinité d'arguments.
+{{ :mite_rda1.png |}}
-D'un point de vue technique, les arguments sont des variables comme les autres et s'utilisent donc de la même façon. La seule différence est qu'elles sont uniquement disponible à l'interieur de votre fonction. Leur valeur sera déterminée au moment où votre fonction sera appelée.
-Commençons avec une fonction très basique que nous appelerons **print_number** qui va prendre un nombre en entrée et l'afficher.
+</hidden>
-<code rsplus>
-print_number <- function(number) {
-  print(number)
-}
-</code>
-Maintenant, pour l'utiliser, nous l'utilisons comme n'importe quelle autre fonction.
-<code rsplus>
+=====2.2 RDA partielle=====
-print_number(2)
-print_number(231)
-</code>
-Nous pouvons utiliser plus d'un argument. Par exemple, créons une fonction qui prend un premier nombre (number1), l'additionne à un second (number2), multiplie le résultat par un troisième nombre (number3) et enfin affiche le résultat.
+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.
-<code rsplus>
+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.
-operations <- function(number1, number2, number3) {
-  result <- (number1 + number2) * number3
-  print(result)
-}
-operations(1, 2, 3)
+Avec R, la RDA partielle est effectuée en utilisant le fonction rda et en ajoutant un argument correspondant aux co-variables :
-operations(17, 23, 2)
-</code>
-La partie expression de notre fonction - son corps - peut virtuellement contenir n'importe quoi. On peut y retrouver des simples expression R, des boucles, des conditions if/else, même d'autres fonctions.
------
+<code rsplus | RDA partielle avec rda()>
-**Défi 5**
+#Divisez le tableau de données environnementales en deux:
+envtopo <- env[, c(1:3)] # Physiographie : tableau 1
-En utilisant ce que vous avez vu précédemment sur les structures de contrôle, créez une fonction appelée **print_animal** qui prend un animal en argument et donne les résultats suivants:
+names(envtopo)
+envchem <- env[, c(4:10)] # Physico-chimie de l'eau : tableau 2
-<code rsplus>
+names(envchem)
-Scruffy <- "dog"
-Paws <- "cat"
+#FAire la RDA partielle
+spechem.physio=rda(spe.hel, envchem, envtopo)
-print_animal(Scruffy)
+summary(spechem.physio, display=NULL)
-[1] "woof"
+#ou
+spechem.physio2=rda(spe.hel ~ pH + dur + pho + nit + amm + oxy + dbo
-print_animal(Paws)
+                    + Condition(alt + pen + deb), data=env)
-[1] "meow"
+#Extraire les résultats
+summary(spechem.physio, display=NULL)
+#Calcul du R2 ajusté de la RDA partielle
+(R2adj <- RsquareAdj(spechem.physio)$adj.r.squared)
+#Test de la significativité des axes
+anova.cca(spechem.physio, step=1000)
+anova.cca(spechem.physio2, step=1000, by="axis")
+#Faire les triplots
+#Cadrage 1
+windows(title="Partial RDA scaling 1")
+plot(spechem.physio, scaling=1, main="Triplot partial RDA - scaling 1", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
+points(scores(spechem.physio, display="sites", choices=c(1,2), scaling=1),
+       pch=21, col="black", bg="steelblue", cex=1.2)
+arrows(0,0,
+      scores(spechem.physio, display="species", choices=c(1), scaling=1),
+      scores(spechem.physio, display="species", choices=c(2), scaling=1),
+      col="black")
+text(scores(spechem.physio, display="species", choices=c(1), scaling=1),
+     scores(spechem.physio, display="species", choices=c(2), scaling=1),
+     labels=rownames(scores(spechem.physio, display="species", scaling=1)),
+     col="black", cex=0.8)
+arrows(0,0,
+      scores(spechem.physio, display="bp", choices=c(1), scaling=1),
+      scores(spechem.physio, display="bp", choices=c(2), scaling=1),
+      col="red")
+text(scores(spechem.physio, display="bp", choices=c(1), scaling=1)+0.05,
+     scores(spechem.physio, display="bp", choices=c(2), scaling=1)+0.05,
+     labels=rownames(scores(spechem.physio, display="bp", choices=c(2), scaling=1)),
+     col="red", cex=1)
+#Cadrage 2
+windows(title="Partial RDA scaling 2")
+plot(spechem.physio, scaling=2, main="Triplot partial RDA - scaling 2", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
+points(scores(spechem.physio, display="sites", choices=c(1,2), scaling=2),
+       pch=21, col="black", bg="steelblue", cex=1.2)
+arrows(0,0,
+      scores(spechem.physio, display="species", choices=c(1), scaling=2),
+      scores(spechem.physio, display="species", choices=c(2), scaling=2),
+      col="black")
+text(scores(spechem.physio, display="species", choices=c(1), scaling=2),
+     scores(spechem.physio, display="species", choices=c(2), scaling=2),
+     labels=rownames(scores(spechem.physio, display="species", scaling=2)),
+     col="black", cex=0.8)
+arrows(0,0,
+       scores(spechem.physio, display="bp", choices=c(1), scaling=2),
+       scores(spechem.physio, display="bp", choices=c(2), scaling=2),
+       col="red")
+text(scores(spechem.physio, display="bp", choices=c(1), scaling=2)+0.05,
+     scores(spechem.physio, display="bp", choices=c(2), scaling=2)+0.05,
+     labels=rownames(scores(spechem.physio, display="bp", choices=c(2), scaling=2)),
+     col="red", cex=1)
 </code>
-++++ Défi 5 : Réponse|
-<file rsplus| Défi 5 Réponse>
-print_animal <- function(animal) {
-  if (animal == "dog") {
-    print("woof")
-  } else if (animal == "cat") {
-    print("meow")
-  }
-}
-</file>
-++++
------
-Les arguments peuvent également être optionnels, auquel cas on peut leur donner une valeur par défaut. Ceci peut s'avérer utile si l'on prévoit d'utiliser fréquemment une fonction avec les mêmes paramètres pour éviter d'avoir à les réécrire à chaque fois, mais si l'on veut tout de même garder la possibilité de changer leur valeur si nécessaire.
+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>
-operations <- function(number1, number2, number3=3) {
-  result <- (number1 + number2) * number3
-  print(result)
-}
-operations(1, 2, 3) # devient équivalent a
+{{ :doubs_rdapart2.png |}}
-operations(1, 2)
-operations(1, 2, 2) # on peut toujours changer la valeur de number3 si necessaire
-</code>
-R fournit également l'argument spécial "**...**". Ceci vous permet de dire à R que votre fonction peut accepter un nombre indéfini d'arguments. ceci est utile pour deux raisons principalement:
+**Défi 2**
-  * Passer des arguments à une autre fonction utilisée dans votre fonction. Ceci vous permet d'utiliser tous les arguments d'autres fonctions sans avoir à les définir un à une lors de la création de votre fonction. Par exemple, créons une fonction à partir de notre exemple précédent où nous traçons l'absorption de CO2 en fonction de la concentration. Nous allons tracer nos graphes avec deux couleurs différentes selon la région. Les paramètres de plot() et points() seront passés via "...".
+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.
-<code rsplus>
+**Défi 2**: Solution
-plot.CO2 <- function(CO2, ...) {
+<hidden>
-  plot(x=CO2$conc, y=CO2$uptake, type="n", ...)  # On utilise ... pour passer les arguments a plot().
-  for (i in 1:length(CO2[,1])){
+Votre code pourrait ressembler à ceci:
-     if (CO2$Type[i] == "Quebec") {
-       points(CO2$conc[i], CO2$uptake[i], col="red", type="p", ...) # idem pour points()
-     } else if (CO2$Type[i] == "Mississippi") {
-       points(CO2$conc[i], CO2$uptake[i], col="blue", type="p", ...) # idem pour points()
-     }
-  }
-}
-plot.CO2(CO2, cex.lab=1.4, xlab="CO2 concentration", ylab="CO2 uptake")
+<code rsplus | RDA partielle sur les données d'acariens>
-plot.CO2(CO2, cex.lab=1.4, xlab="CO2 concentration", ylab="CO2 uptake", pch=20)
+mite.spe.subs=rda(mite.spe.hel ~ Shrub + Topo
+                  + 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)
+arrows(0,0,
+       scores(mite.spe.subs, display="species", choices=c(1), scaling=1),
+       scores(mite.spe.subs, display="species", choices=c(2), scaling=1),
+       col="black")
+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="black", 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)
+arrows(0,0,
+       scores(mite.spe.subs, display="species", choices=c(1), scaling=2),
+       scores(mite.spe.subs, display="species", choices=c(2), scaling=2),
+       col="black")
+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="black", 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>
-> NOTE : Il est très important de garder en tête que lorsque l'on passe des arguments en utilisant "...", il est nécessaire de spécifier le nom de ces derniers pour éviter toute confusion.
+{{ :partial_rda_ch6.png |}}
-  * Autoriser l'utilisateur à entrer un nombre indéfini d'arguments. La valeur de chaque argument devra alors être récupérée manuellement. Cela se fait en transformant "..." en liste et en itérant dessus. Par exemple, créons une fonction somme qui accepte un nombre indéfini d'arguments.
+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%.
+</hidden>
+=====2.3 Partitionnement de la variation par RDA partielles=====
-<code rsplus>
+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 :
-sum2 <- function(...){
-  args <- list(...)
-  result <- 0
-  for (i in args)  {
-    result <- result + i
-  }
-  return (result)
-}
-sum2(2, 3)
+- 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.
-sum2(2, 4, 5, 7688, 1)
-</code>
-=== Valeur de retour ===
+{{ :vennd_varpart.png |}}
-Comme montré dans les exemples précédents, si nous voulons être capable de sauvegarder le résultat de notre fonction pour pouvoir l'utiliser plus tard, nous devons le renvoyer à la fin de la fonction en utilisant return(). Il est important de garder en tête qu'une seule valeur de retour peut être renvoyée par une fonction. Si vous désirez renvoyer plus d'un objet, vous devrez utiliser des objets composites tels que des listes ou des dataframes.
+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)
-Par ailleurs, il est important de noter que l'execution de la fonction se termine dès qu'elle atteint le mot clé return().
-<code rsplus>
+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.
-returntest <- function(a, b) {
+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.
-  return (a) # La fonction s'arrete la
-  a <- a + b # Pas interprete
-  return (a + b) # Pas interprete
-}
-returntest(2, 3) # Par defaut R affiche la valeur de retour de la fonction
-c <- returntest(2, 3) # Pour la sauvegarder, ne pas oublier de l'assigner a une variable
+<code rsplus | varpart()>
-c
+?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>
------
-**Défi 6**
-En utilisant ce que vous avez appris jusqu'ici sur les fonctions et les structures de contrôle, créez une fonction **bigsum** qui prend deux arguments **a** et **b** et :
+Le sommaire ressemble à ceci:
-  * retourne 0 si la somme de a et b est strictement inférieure à 50
+{{ :varpart_output.png |}}
-  * retourne la somme de a et b sinon
+{{ :varpart_output_venn.png |}}
-++++ Défi 6 : Réponse|
+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.
-<file rsplus| Défi6 Réponse>
-bigsum <- function(a, b) {
-  result <- a + b
-  if (result < 50) {
-    return(0)
-  } else
-    return (result)
-  }
-}
-</file>
-++++
------
-==== Accessibilé des variables ====
-Lorsque l'on travaille avec des structures de contrôle et des fonctions, il est essentiel de toujours savoir où sont nos variables, si elles sont définies et accessibles.
+Il est également possible d’effectuer un partitionnement de la variation à partir de jeux de données explicatives contenant uniquement des variables significatives :
-Voici quelques conseils à garder en tête.
-  * Les variables définies à l'intérieur d'une fonction ne sont pas accessibles en dehors.
-  * Les variables définies en dehors d'une fonction sont accessibles à l'intérieur. Cependant, ce n'est **JAMAIS** une bonne idée de les utiliser à l'intérieur, car votre fonction pourrait arrêter de marcher si la variable n'existe pas.
-<code rsplus>
+<code rsplus | varpart avec variables explicatives>
+#RDA des variables physico-chimiques
+spe.chem <- rda(spe.hel~., data=envchem)
+#Sélection des variables significatives
+R2a.all.chem <- RsquareAdj(spe.chem)$adj.r.squared
+ordiR2step(rda(spe.hel~1, data=envchem),
+           scope= formula(spe.chem), direction= "forward", R2scope=TRUE, pstep=1000)
+names(envchem)
+(envchem.pars <- envchem[, c( 4, 6, 7 )])
+#RDA avec les autres variables significatives
+spe.topo <- rda(spe.hel~., data=envtopo)
+R2a.all.topo <- RsquareAdj(spe.topo)$adj.r.squared
+ordiR2step(rda(spe.hel~1, data=envtopo),
+           scope= formula(spe.topo), direction= "forward", R2scope=TRUE, pstep=1000)
+names(envtopo)
+envtopo.pars <- envtopo[, c(1,2)]
+#Varpart
+spe.part <- varpart(spe.hel, envchem.pars, envtopo.pars)
+windows(title="Variation partitioning - parsimonious subsets")
+plot(spe.part, digits=2)
+#Tests de significativité des fractions
+anova.cca(rda(spe.hel, envchem.pars), step=1000) # Test of fractions [a+b]
+anova.cca(rda(spe.hel, envtopo.pars), step=1000) # Test of fractions [b+c]
+env.pars <- cbind(envchem.pars, envtopo.pars)
+anova.cca(rda(spe.hel, env.pars), step=1000) # Test of fractions [a+b+c]
+anova.cca(rda(spe.hel, envchem.pars, envtopo.pars), step=1000) # Test of fraction [a]
+anova.cca(rda(spe.hel, envtopo.pars, envchem.pars), step=1000) # Test of fraction [c]
+</code>
-rm(list=ls()) # supprimons toutes nos variables pour eviter toute confusion
-var1 <- 3     # var1 est definie en dehors de la fonction
+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).
-vartest <- function() {
-  a <- 4      # a est definie a l'interieur
-  print(a)    # affiche a
-  print(var1) # affiche var1
-}
-a             # affiche a. Ca ne marche pas, a est seulement visible dans la fonction
-vartest()     # appeller vartest() affiche a et var1
-rm(var1)      # supprime var1
-vartest()     # la fonction ne marche plus, var1 n'existe plus
-</code>
+**Défi 3** Effectuer le partitionnement de la variation de l’abondance des espèces du tableau mite en utilisant un premier jeu de variables explicatives relatives au substrat (SubsDens, WaterCont and Substrate) et un second jeu de variables explicatives contenant les autres variables significatives (Shrud and Topo) Quelle est la proportion de variance expliquée par chaque jeu de données ? Quelles sont les fractions significatives ?
+**Défi 3**: Solution
-À la place, utilisez des arguments!!Dans une fonction, le nom des arguments remplacera le nom des autres variables.
+<hidden>
+Votre code pourrait ressembler à ceci:
-<code rsplus>
+<code rsplus | varpart avec variables significatives>
-var1 <- 3     # var1 est definie en dehors de la fonction
+str(mite.env)
-vartest <- function(var1) {
+(mite.subs=mite.env[,c(1,2,3)]) #Premier jeu de données
-  print(var1) # affiche var1
+(mite.other=mite.env[,c(4,5)]) #Deuxième jeu de données
-}
-vartest(8)     # Dans notre fonction, var1 est maintenant notre argument et prend sa valeur
+#RDA sur mite.subs
-var1          # var1 a toujours la meme valeur
+rda.mite.subs <- rda(mite.spe.hel~., data=mite.subs)
-</code>
+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>
-Faites très attention lorsque vous créez des variables à l'intérieur d'une condition car la variable pourrait ne jamais être créée et causer des erreurs
+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).
-<code rsplus>
+</hidden>
-a <- 3
-if (a > 5) {
-  b <- 2
-}
-a + b    # Erreur! b n'est pas cree
-</code>
-En général, il est recommandé de définir les variables en dehors des conditions et de les modifier leur valeur à l'intérieur pour éviter tout problème
-<code rsplus>
+======3. Arbre de régression multivarié======
-a <- 3
-b <- 0
-if (a > 5) {
-  b <- 2
-}
-a + b
-</code>
-===== Bonnes pratiques et comment accélerer vore code =====
-==== Bonnes pratiques ====
+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.
-Voici quelques conseils de programmation qui peuvent vous faciliter la vie, vous aider à avoir un code plus lisible et qui rendent le partage et la réutilisation de votre code bien moins difficile. Avoir un code facile à lire aide à réduire le temps que vous passeriez à essayer de le comprendre, donc ce n'est jamais du temps perdu.
+Avantages du MRT par rapport à la RDA:
+    * n'assume pas de relation linéaire entre les abondances d'espèces et les variables environnementales (quantitatives ou qualitatives),
+    * la méthode est robuste en présence de valeurs manquantes,
+    * la méthode est robuste en présence de colinéarité entre les descripteurs,
+    * les MRTs ne nécessitent pas la transformation des variables explicatives (les valeurs brutes peuvent être utilisées).
+    * l'arbre obtenu est facile à interpréter pour un public scientifique ou néophyte.
+Le MRT divise les données en groupes ayant des compositions en espèce semblables et caractérisés par des variables environnementales. La méthode implique deux volets s'effectuant en parallèle: 1) la construction de l'arbre et 2) la sélection de la partition finale optimale par validation croisée. La fonction mvpart() de la librairie mvpart {} permet de créer les arbres de régression.
-=== Gardez un code beau et propre ===
+Vocabulaire lié aux MRTs:
-L'une des choses qui aide le plus lors que l'on lit un code informatique, c'est d'avoir un code bien formatté, bien espacé et bien indenté. Certains standards de programmation existent pour vous aider à obtenir une plus grande consistance, mais cela dépend souvent ultimement des préférences de chacun. Voici quelques trucs pour vous aider:
+Feuille: Groupe terminal de sites
-  * Mettez des espaces avant et après vos opérateurs
-  * Utilisez toujours le même opérateur d'assignation. `<-` est souvent préférable, `=` est ok mais ne changez pas tout le temps entre les deux
-  * Utilisez des crochets pour encadrer vos structures de contrôle, même si c'est juste pour une ligne. Chaque ligne de code à l'intérieur des crochets devrait être indentée d'au moins deux espaces. Les crochets de fermeture occupent généralement leur propre ligne, sauf s'ils précèdent une condition **else**. Ces pratiques adent grandement lorsque l'on veut déterminer ou l'on se trouve, en particulier si l'on a beaucoup de conditions/boucles imbriquées les unes dans les autres.
-  * Définissez chaque variable sur sa propre ligne
-Voici un exemple de code difficile à lire
+Noeud: Point où les données se divisent en deux groupes. Chaque noeud est caractérisé par une valeur d'une variable explicative.
-<code rsplus>
-a<-4;b=3
-if(a<b){
-if(a==0)print("a zero") } else {
-if(b==0){print("b zero")} else print(b)}
-</code>
-Voici une version plus aisée. Elle prend plus d'espace mais le flot du code est plus facile à voir.
+Branche: Chaque lignée formée par un noeud.
-<code rsplus>
-a <- 4
-b <- 3
-if(a < b){
-  if(a == 0) {
-    print("a zero")
-  }
-} else {
-  if(b == 0){
-    print("b zero")
-  } else {
-    print(b)
-  }
-}
-</code>
-Certains guides de style peuvent être trouvés sur internet. En voici un exemple: [[https://google-styleguide.googlecode.com/svn/trunk/Rguide.xml|https://google-styleguide.googlecode.com/svn/trunk/Rguide.xml]]
+- Partitionnement des données sous contrainte (construction de l'arbre)
-=== Utilisez des fonctions si possible ===
+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.
-Maintenant que vous savez comment créer une fonction, n'hésitez pas a les utiliser. Dès que vous aperceve une portion de code qui est répétée plus que deux fois dans votre code, vous devriez vous dire "Hmmm... est ce que ça ne serait pas mieux d'écrire une fonction à la place?". Si seulement une portion de ce code change, essayez de penser à des façons d'utiliser des arguments dans une fonction à la place. Ceci vous aidera à réduire le nombre d'erreurs réalisées en faisant des copier/coller, réduira le temps passé a les corriger et facilitera les modifications futures éventuelles.
-Par exemple, modifions l'exemple du défi 3 et supposons que toutes les absorptions de CO2 du Mississipi étaient surestimées de 20 et que celles du Québec étaient sous-estimées de 50. Nous pourrions écrire ceci.
-<code rsplus>
+- Validation croisée et élagage de l'arbre
-for (i in 1:length(CO2[,1])) {
-  if(CO2$Type[i] == "Mississippi") {
-    CO2$conc[i] <- CO2$conc[i] - 20
-  }
-}
-for (i in 1:length(CO2[,1])) {
-  if(CO2$Type[i] == "Quebec") {
-    CO2$conc[i] <- CO2$conc[i] + 50
-  }
-}
-</code>
-Ou alors nous pourrions faire ceci à la place.
+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.
-<code rsplus>
+Créons un arbre de régression multivarié sur les données de la rivière Doubs.
-recalibrate <- function(CO2, type, bias) {
-  for (i in 1:nrow(CO2)) {
-    if(CO2$Type[i] == type) {
-      CO2$conc[i] <- CO2$conc[i] + bias
-    }
-  }
-  # On doit retourner notre nouveau jeu de donnees car l'original n'est pas modifie
-  return (CO2)
-}
-newCO2 <- recalibrate(CO2, "Mississipi", -20)
-# Notez que nous recalibrons ici notre variable newCO2 parce que le CO2 original n'est pas modifie
-newCO2 <- recalibrate(newCO2, "Quebec", +50)
-</code>
+<code rsplus | mvpart()>
+?mvpart
-Et maintenant, nous réalisons que ce que l'on a modifié dans nos exemples précédents n'était l'absorption, mais la concentration... Maintenant on doit changer toutes les occurences de //CO2\$conc[i]// par //CO2\$uptake[i]//. Dans le premier cas, cela veut dire que l'on a à le changer 4 fois, contre seulement deux fois dans notre fonction! (Bon ok, vous pouvez vous dire ici que cela ne vaut pas vraiment le coup, qu'avec un simple rechercher/remplacer vous le faites super vite et effectivement vous auriez raison. Mais ce n'est qu'un simple exemple! Imaginez si vous aviez à le remplacer 10 fois au lieu de 2. Un bon programmeur est un programmeur paresseuz. Et aussi, admettez le, ca en jette plus avec une fonction...)
+#Preparer les données: enlever la variable “distance from source”
+env <- subset(env, select = -das)
-=== Donnez des noms qui ont du sens à vos variables et fonctions ===
+# Créer l'arbre
+doubs.mrt<-mvpart(data.matrix(spe.hel)~.,env,
+legend=FALSE,margin=0.01,cp=0,xv="pick",xval=nrow(spe.hel),xvmult=100,which=4)
+</code>
-Ceci aide à voir au premier coup d'oeil qui fait quoi. Soyez encore plus prudents dans le choix du nom de vos arguments quand vous créez une fonction car c'est ce que les utilisateurs voient. Toutefois il est parfois judicieux de choisir des noms courts pour éviter d'avoir à les taper tout le temps et ainsi éviter les fautes de frappe, donc un bon équilibre doit être choisi.
-Voici ce à quoi notre exemple précédent pourrait ressembler avec des noms vagues. Comprendre ce que cette fonction fait demande maintenant un peu plus d'efforts.
+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).
-<code rsplus>
-rc <- function(c, t, b) {
-  for (i in 1:nrow(c)) {
-    if(c$Type[i] == t) {
-      c$uptake[i] <- c$uptake[i] + b
-    }
-  }
-  return (c)
-}
-</code>
-=== Commentaires ===
+{{ :cross_validation.png |}}
-Même avec des noms évidents, ce n'est jamais une mauvaise chose d'ajouter des commentaires pour décrire tout ce que votre code fait, que ce soit le but de la fonction, comment utiliser ses arguments ou une description détaillée de la fonction étape par étape.
+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.
-<code rsplus>
+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.
-## Recalibre le jeu de donnees CO2 en modifiant l'absorption de CO2
-## d'une valeur fixe selon la region
-# Arguments
-# CO2: le jeu de donnees CO2
-# type: le type de donnees qui doivent etre recalibrees. Valeurs: "Mississippi" ou "Quebec"
-# bias: la quantite a ajouter a l'absorption. Utilisez des valeurs negative pour les surestimations
-recalibrate <- function(CO2, type, bias) {
-  for (i in 1:nrow(CO2)) {
-    if(CO2$Type[i] == type) {
-      CO2$uptake[i] <- CO2$uptake[i] + bias
-    }
-  }
-  # On doit retourner notre nouveau jeu de donnees car l'original n'est pas modifie
-  return (CO2)
-}
-</code>
+{{ :mrt_1se.png |}}
-===== 3. Accélerer votre code =====
+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.
-Ici sont présentés quelques trucs pour programmer de manière plus efficace avec R et vous aider à obtenir de meilleures performances et un code plus rapide. Toutefois, avant d'optimiser votre code, il est important de s'assurer que vous avez d'abord un code qui fonction. Un code lent qui fonctionne sera toujours meilleur qu'un code rapide qui ne marche pas. Par ailleurs, parfois, cela ne sert à rien d'optimiser. Passer 2 heures à réécrire des lignes de code pour gagner quelques secondes à l'execution n'est pas forcément la solution la plus efficace.
+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.
-==== Avant de commencer : évaluer nos performances ====
+Nous pouvons comparer cet arbre avec la solution minimisant le CVRE (contenant 10 feuilles), ou avec une solution intermédiaire (ex. avec 4 feuilles).
-Si l'on veut optimiser notre code, la première étape est de savoir combien de temps chaque tâche prend.
-La façon la plus simple de le faire est d'utiliser la fonction system.time(//expression//)
+<code rsplus | Comparaison des arbres>
+# En utilisant le critère CVRE
+doubs.mrt.cvre<-mvpart(data.matrix(spe.hel)~.,env,
+legend=FALSE,margin=0.01,cp=0,xv="pick",xval=nrow(spe.hel),xvmult=100,which=4)
-<code rsplus>
+# En choisissant un arbre avec 4 feuilles
-system.time({
+doubs.mrt.4<-mvpart(data.matrix(spe.hel)~.,env,
-a <- 0
+legend=FALSE,margin=0.01,cp=0,xv="pick",xval=nrow(spe.hel),xvmult=100,which=4)
-  for (i in 1:1000) {
-    a <-  a + i
-  }
-})
 </code>
-Notez que la plupart du temps, R travaille vraiment rapidement et il faut avoir des tâches qui demandent vraiment beaucoup de puissance ou le temps risque de ne même pas être enregistré. C'est pourquoi il est recommandé de répéter plusieurs fois la tâche que l'on veut évaluer ou alors de travailler sur de gros jeux de données.
+{{ :mrt_cvre.png |}}
+{{ :mrt_4.png |}}
-<code rsplus>
+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.
-system.time(replicate(1000, {
-  a <- 0
+Le sommaire de la fonction contient beaucoup d'informations.
-  for (i in 1:1000) {
-    a <-  a + i
-  }
-}))
-</code>
-Un autre outil simple et utile est la fonction Rprof(). L'avantage principal de Rprof() est qu'elle enregistre de l'infromation sur le temps passé dans chaque fonction dans un fichier auquel on peut accéder plus tard. Voici comment l'utiliser.
-<code rsplus>
+<code rsplus | sommaire du MRT>
-Rprof("profile.txt")  # on peut changer profile.txt par le nom de fichier desire
+summary(doubs.mrt)
-for (i in 1:1000) {
-    a <- 0
-    for (i in 1:1000) {
-      a <-  a + i
-    }
-  }
-Rprof()               # Ceci termine le profilage
-summaryRprof("profile.txt")  # Utilisez le nom de fichier enregistre precedemment pour afficher le resume
 </code>
-Enfin, si vous voulez comparer l'efficacite de plusieurs fonctions côte à côte, un très bon outil est le package microbenchmark
+{{ :doubs_mrt_summary.png |}}
-<code rsplus>
+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.
-install.packages("microbenchmark")
-library(microbenchmark)
-f1 <- function() {
-  a <- 0
-  for (i in 1:1000) {
-    a <-  a + i
-  }
-}
-microbenchmark(f1(), times=1000) # l'argument times nous permet de determiner le nombre d'iterations voulues
+<code rsplus | Trouver des espèces discriminantes et indicatrices>
-</code>
+# Trouver les espèces discriminantes de l'arbre de régression
+doubs.mrt.wrap<-MRT(doubs.mrt,percent=10,species=colnames(spe.hel))
+summary(doubs.mrt.wrap)
-==== Première étape : réfléchir un peu! ====
+# Extraire les p-values des valeurs indval
+doubs.mrt.indval<-indval(spe.hel,doubs.mrt$where)
+doubs.mrt.indval$pval
-Si vous regardez attentivement votre code, souvent vous realiserez qu'il existe d'autres façons plus simples, plus efficaces de faire ce que vous désirez et que certaines opérations peuvent facilement être supprimées pour gagner du temps.
+# Extraire les espèces indicatrices à chaque noeud et leur valeur indval
+doubs.mrt.indval$maxcls[which(doubs.mrt.indval$pval<=0.05)]
+doubs.mrt.indval$indcls[which(doubs.mrt.indval$pval<=0.05)]
+</code>
-Par exemple, créons une fonction qui prend un nombre **a**. Nous allons ajouter **a** à chaque nombre de 1 à 100, et si **a** est inférieur à 5, alors nous ajouterons 2*a à la place. Ensuite, nous additionnerons ensemble tous les elements de la sequence.
+{{ ::doubs_mrt_discriminant.png |}}
+{{ ::doubs_mrt_finalpart.png |}}
-Voici une façon de le faire
+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.
-<code rsplus>
+{{ ::doubs_mrt_indval.png |}}
-f2 <- function(a) {
-  # initialisation du résultat
-  result <- 0
-  # on itere sur la sequence de 1 à 100
-  for (i in 1:100) {
-    if (a < 5) {
-      # a est < 5, on ajoute 2*a a la sequence. On met le tout dans result
-      result <- result + i + (2*a)
-    } else {
-      # a est >= 5, on n'ajoute que a
-      result <- result + i + a
-    }
-  }
-  return(result)
-}
-f2(4)
-</code>
-Notre fonction fait ce que l'on désire et est une solution tout à fait acceptable. Cependant, nous avons plein d'étapes inutiles dans notre code. Par exemple, nous n'avons pas besoin d'effectuer notre condition à chaque itération, car le résultat sera toujours le même.
+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).
-On peut donc la sortir de la boucle.
-<code rsplus>
+**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?
-f3 <- function(a) {
-  # initialisation du résultat
-  result <- 0
-  # On verifie si a < 5, si oui, a vaut maintenant 2*a
-  if (a < 5) {
-   a <- 2 * a
-  }
-  # nous n'avons meme pas besoin de tester l'alternative puisque a reste identique
-  # on itere sur la sequence de 1 à 100
-  for (i in 1:100) {
-      result <- result + i + a
-  }
-  return(result)
-}
-f3(4)
+**Défi 4** - Solution
-microbenchmark(f2(4),
-               f3(4), times=1000)
-</code>
+<hidden>
-Nous avons effectué seulement une simple modification mais nous avons ici réussi à accélerer notre code d'environ 40% (les résultats peuvent varier selon les ordinateurs). De plus, notre code est plus facile à lire et à comprendre. Parfois, on peut gagner à la fois vitesse et lisibilité juste en réflechissant à la place de nos conditions et à ce qu'elles testent.
+<code rsplus | Créé un MRT avec les données sur les acariens>
+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)
-Mais en utilisant les forces de R, on peut faire encore mieux!
+mite.mrt.wrap<-MRT(mite.mrt,percent=10,species=colnames(mite.spe.hel))
+summary(mite.mrt.wrap)
-<code rsplus>
+mite.mrt.indval<-indval(mite.spe.hel,mite.mrt$where)
-f4 <- function(a, n) {
+mite.mrt.indval$pval
-  result <- 0
-  if (a < 5) {
-    a <- a + 1
-  }
-  result <- sum(1:n + a)
-  return(result)
-}
-f4(4)
-microbenchmark(f3(4),
-               f4(4), times=1000)
+mite.mrt.indval$maxcls[which(mite.mrt.indval$pval<=0.05)]
+mite.mrt.indval$indcls[which(mite.mrt.indval$pval<=0.05)]
 </code>
-Wow, ici notre modification est beaucoup plus efficace... Mais qu'est ce qui s'est passé exactement? Ceci nous amène à notre prochain point.
+.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.
-==== Vectorisation ====
+</hidden>
-Cette partie est un rappel de choses que vous avez probablement déjà vues lors des premiers ateliers. Cependant, ces notions sont souvent oubliées et les mauvaises performances de R peuvent fréquemment être attribuées à une mauvaise vectorisation.
-R est conçu pour travailler avec les vecteurs et par conséquent, de nombreuses fonctions sont optimisées pour la vectorisation. Pour comprendre ceci, il est d'abord important de comprendre comment R fonctionne. R est un langage interprété, ce qui veut dire que lorsque vous executez votre code R, en réalité vous envoyez vos instructions à des fonctions programmées dans un autre langage (le langage C). Ceci ralentit l'execution de vos programmes puisque le code R doit d'abord être décodé puis envoyé à d'autres fonctions. Lorsque vous créez une boucle, vous devez décoder chaque itération puis la transférer.
-Les fonctions vectorisées d'un autre côté sont des fonctions qui travaillent directement avec des vecteurs. Elles executent elles aussi une boucle sur votre vecteur, mais la grosse différence est qu'elles le font directement en C, ce qui est beaucoup plus rapide. la fonction sum() est un exemple de fonction vectorisée.
-L'un des plus gros challenges de R est d'apprendre à penser et à programmer avec des vecteurs et non avec des éléments simples. Par exemple, la plupart des opérations de base peuvent être faites sur des vecteurs.
-<code rsplus>
+======4. Analyse discriminante linéaire======
-v1 <- 1:5
-v2 <- 2:6
-v3 <- 1:3
-v1 + 2      # Addition sur un vecteur : ajoute 2 a tous les elements
-v1 + v2     # Ajoute chaque element de v2 a v1
-v1 + v3     # v1 et v3 ne sont pas de la meme taille, on recommence a additionner a partir du debut de v3
-sum(v1)     # Additionne tous les elements de v1 ensemble
-sum(v1, v2) # Fait la somme de tous les elements de v1 et v2
-mean(v1)    # Fait la moyenne de v1
-mean(c(v1, v2)) # Moyenne des elements de v1 et v2. Contrairement a sum(), on doit les combiner avant
-</code>
-=== Extraire des sous-ensembles ===
+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.
-Pour vectoriser efficacement, il est également important d'être capable d'extraire des valeurs de nos données rapidement.
+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.
-R offre des outils de sélection de sous-ensembles pour appliquer un traitement sur des élements spécifiques de vecteurs ou de dataframes qui sont parfois plus efficaces et plus faciles à écrire que des boucles et des conditions.
-L'extraction de sous-ensembles est faite vis les opérateurs **[** et **$** (pour un dataframe). Nous pouvons insérer directement nos condition dans la partie **[]** pour extraire rapidement des valeurs de nos données. Il est également possible d'utiliser la fonction **which()** pour tester une condition. which() retourne les indexs des élements qui remplissent la condition.
+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.
-<code rsplus>
+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.
-v1 <- 1:10
-v1[7]      # Extrait la 7eme valeur
-v1[v1 > 5] # Extrait les valeurs > 5 seulement
-v1[which(v1 > 5)]  # pareil que precedemment
-</code>
-Dans les dataframe, **$** permet d'acceder à une colonne par nom. Nous pouvons également le faire en fournissant directement le nom de la colonne.
+<code rsplus | Charger les données spatiales et classifier les abondances de poissons par site en fonction de la latitude>
-<code rsplus>
+#Charger les données spatiales pour déterminer les groupes
-CO2 <- read.csv("co2_good.csv")
+spa <- read.csv ('http://www.davidzeleny.net/anadat-r/data-download/DoubsSpa.csv', row.names = 1)
-CO2$Type  # Affiche la colonne Type
+spa <- spa[,-8]
-CO2[, "Type"] # Idem
-CO2[CO2$Type == "Quebec", ] #Extrait toutes les lignes de CO2 dont le Type est "Quebec"
-</code>
+#Visualiser la matrice de données
+View (spa)
------
+#Ajouter les numéros de site
-**Défi 7**
+numbers<-(1:30)
+numbers<-numbers[!numbers%in%8]
+spa$site<-numbers
-Créez une nouvelle fonction recalibrate2(), qui est une réécriture de la fonction recalibrate vue précedemment, en utilisant des techniques de vectorisation et d'extraction de sous-ensembles. La nouvelle fonction ne devrait pas faire plus de 3 lignes.\\
+#Faire des groupes en fonction de la latitudey<82=group1, 82<y<156=group2, y>156=group3
-Rappel:
+spa.group<-ddply(.data=spa, .variables=.(x, y, site), .fun= summarise, group = if(y <= 82) 1 else if (y <= 156) 2 else 3)
-<code rsplus>
+#Ordonner par site
-recalibrate <- function(CO2, type, bias) {
+spa.group<-spa.group[with(spa.group, order(site)), ]
-  for (i in 1:nrow(CO2)) {
-    if(CO2$Type[i] == type) {
-      CO2$uptake[i] <- CO2$uptake[i] + bias
-    }
-  }
-  return (CO2)
-}
 </code>
-++++ Défi 7 : Réponse|
-<code rsplus>
-recalibrate2 <- function(CO2, type, bias) {
-  # D'abord recuperons les index des donnees avec le bon type
-  # Astuce de reflexion : puisque nous utilisons les index deux fois ci dessous, plutot que d'utiliser which()
-  # deux fois, sauvegardons le resultat du premier appel!
-  idx <- which(CO2$Type == type)
-  # Modifions uniquement les donnees concernees en utilisant les index.
-  CO2$uptake[idx] <- CO2$uptake[idx] + bias
-  return (CO2)
-}
-# Verifions que tous les résultats sont les memes
+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.
-all.equal(recalibrate(CO2, "Quebec", 20), recalibrate2(CO2, "Quebec", 20))
-# Verifions que c'est effectivement plus rapide
+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.
-microbenchmark(recalibrate(CO2, "Quebec", 20),
-               recalibrate2(CO2, "Quebec", 20))
-</code>
-++++
------
-==== Les objets qui grossissent ====
-Vectoriser est une bonne chose, mais cela s'avère parfois difficile et il se peut que cela vous prenne plus de temps que d'écrire une simple boucle. Parfois, les boucles sont absolument nécessaires et il ne faut pas se restreindre de les utiliser. Cependant, dans votre boucle, si vous désirez avoir des performances décentes, vous devrez faire attention aux objets qui grossissent. C'est à dire les objets qui deviennent de plus en plus gros à chaque itération.
+<code rsplus | Faire l'analyse discriminante linéaire (LDA)>
-Illustrons ceci simplement en créant une fonction qui itère sur une séquence et crée un vecteur avec. Nous comparerons deux façons de faire: en laissant notre objet grandir ou en préallouant notre objet résultat et en le modifiant à chaque itération.
+#faire la LDA
+LDA<-lda(env,spa.group[,4])
-<code rsplus>
+#classification des objets en fonction de la LDA
-growing <- function(n) {
+spe.class <- predict(LDA)$class
-  # on declare notre objet resultat
-  result <- NULL
-  for (i in 1:n) {
-    # on cree notre resultat en le faisant grandir a chaque iteration
-    result <- c(result, i)
-  }
-  return(result)
-}
-growing2 <- function(n) {
+#probabilités que les objets appartiennent à chaque groupe a posteriori
-  # on declare notre resultat : ici on cree un vecteur de taille n avec des 0 dedans
+spe.post <- predict(LDA)$posterior
-  result <- numeric(n)
-  for (i in 1:n) {
-    # maintenant on modifie juste la valeur au lieu de recreer le vecteur
-    result[i] <- i
-  }
-  return(result)
-}
-</code>
-Maintenant comparons leurs vitesses respectives
+#tableau des classifications a priori et prédites
+spe.table <- table(spa.group[,4], spe.class)
-<code rsplus>
+#proportion de classification correcte
-system.time({
+diag(prop.table(spe.table, 1))
-  growing(10000)
-})
-system.time({
-  growing2(10000)
-})
 </code>
-Avec un vecteur de 10000 éléments, les vitesses sont encore comparables et prennent moins d'une seconde. Maintenant utilisons 50000 éléments
+{{ :lda_spetable.png?300 |}}
-<code rsplus>
-system.time({
-  growing(50000)
-})
-system.time({
-  growing2(50000)
-})
-</code>
-En multipliant le nombre d'objet par seulement 5, cela nous prend maintenant plusieurs secondes pour créer le vecteur par itération alors que modifier un vecteur prédéfini est toujours quasiment instantané. Qu'est ce qui s'est passé ici? La raison est que lorsque vous appelez une fonction, les arguments sont tout d'abord copiés avant d'être passés à la fonction. Donc lorsque vous écrivez // result <- c(result, i) //, à chaque fois **result** est copié avant d'être passé à c(). Au fur et à mesure que result grossit, à chaque itération cela prend de plus en plus de temps à le copier. Plus l'objet final est gros, plus cela prendra de temps. C'est pourquoi il est toujours préférable de créer votre objet résultat avant votre boucle si vous savez la taille qu'il aura.
+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?
-Ceci est particulièrement valide lorsque l'on travaille avec des dataframes et des fonctions telles que **rbind()** et **cbind()**.
+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.
-Malheuresement, définir au préalable des dataframes ne marche pas si bien que ça et il existe une meilleure façon, bien que plus compliquée. Il faut alors stocker chaque ligne (ou colonne) dans une liste préallouée et ensuite appeler **rbind()**  (ou **cbind()**) sur tous les éléments d'un coup via la fonction **do.call()**. La fonction **do.call()** vous permet d'executer une fonction donnée sur une liste d'arguments. De cette façon, **rbind()** est appelée une seule fois, a la fin, ce qui élimine le problème de copier l'objet au fur et à mesure qu'il grossit.
-<code rsplus>
+Chargez le fichier classifyme.csv. Celui-ci contient des données fictives de 5 nouveaux sites.
-growingdf <- function(n, row) {
-  # predefinission notre dataframe
-  df <- data.frame(numeric(n), character(n), stringsAsFactors=FALSE)
-  for (i in 1:n) {
-    # remplacons la ieme ligne par row
-    df[i,] <- row
-  }
-  return(df)
-}
-growingdf2 <- function(n, row) {
+<code rsplus | Charger les données fictives et prédire le groupement des nouvelles données>
-  # Voici la facon d'allouer une liste a n elements
+#prédire la classification des nouvelles données
-  df <- vector("list", n)
+#charger les nouvelles données
-  for (i in 1:n) {
+classify.me<-read.csv("classifyme.csv", header = T)
-    # on place row dans le ieme element
-    df[[i]] <- row
-  }
-  return(do.call(rbind, df))
-}
-# Stockons notre ligne dans une liste puisque nous avons des elements differents (un nombre et une chaine
+#prédire le groupement des nouvelles données
-# de caracteres)
+predict.group<-predict(LDA, newdata=classify.me)
-row <- list(1, "Hello World")
-microbenchmark(growingdf(5000, row),
+#donner la classification pour chaque site
-               growingdf2(5000, row),
+group.new<-predict.group$class
-               times=10)
 </code>
+{{ :lda_newgroups.png?200 |}}
-==== La famille apply ====
+Nos nouveaux sites, en ordre, ont été classifiés dans les groupes 1, 1, 1, 3 et 3 respectivement.
-Afin d'éviter le problème des objets qui grossisent dans les boucles et pour faciliter l'application de fonctions sur des objets tels que des dataframes, R nous offre ce que nous appelerons les fonctions apply (parce qu'elles possèdent toutes apply dans leur nom...). Il s'agit d'un  groupe de fonctions qui vont exécuter une autre fonction sur un objet d'une type particulier. Leur utilisation diffère seulement selon le type d'objet sur lequel on applique la fonction ou du type de la valeur de retour.
+**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?
-Les fonctions de la famille apply ne sont pas toujours le meilleur choix d'un point de vue performance car elles vont souvent cacher une boucle écrite en R dans leur code. Cependant, elles peuvent grandement réduire le temps de programmation par le confort d'utilisation qu'elles fournissent.
+**Défi 5**: Solution
-L'une des plus populaires est simplement apply(), qui execute une fonction sur les lignes ou les colonnes d'une dataframe ou d'une matrice. La fonction prend 3 arguments principaux:
+<hidden>
-  * l'objet sur lequel on veut appliquer notre fonction
+<code rsplus | LDA on mite data>
-  * le sous-ensemble sur lequel on veut appliquer la fonction. 1 est pour les lignes, 2 pour les colonnes
+mite.xy$site<-seq(1:70)
-  * la fonction à appliquer
+(max(mite.xy[,2])-min(mite.xy[,2]))/4
-  * les arguments éventuels à passer à la fonction fournie
-<code rsplus>
+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)
-df <- data.frame(1:100, 101:200)
+mite.xy.group<-mite.xy.group[with(mite.xy.group, order(site)), ]
-# Somme sur les lignes
-apply(df, 1, sum)
-# Moyenne sur les colonnes
-apply(df, 2, mean)
-# on peut egalement fournir des arguments supplementaire a la fonction
-apply(df, 2, mean, na.rm=TRUE)
-# On peut egalement definir directement une fonction. le premiere argument de cette
-# fonction sera obligatoirement ce sur quoi on veut iterer. Ici chaque ligne est consideree
-# comme un vecteur de nombre comme montre par la fonction str()
-apply(df, 1, function(x){str(x)})
-# On peut egalement ajouter d'autres arguments
-apply(df, 1, function(x, y){x[2] - x[1] + y}, y=5)
-</code>
-Toutes les fonctions apply fonctionnent sur le même modèle. D'un point de vue performance, les plus intéressantes sont sans doute lapply et vapply car ce sont des fontions primitives écrites en C. lapply retourne une liste de la même longueur que notre objet original. vapply vous autorise à spécifier le format de la valeur de retour de votre fonction. Ce peut être un vecteur ou un tableau.
+LDA.mite<-lda(mite.env[,1:2],mite.xy.group[,4])
+mite.class <- predict(LDA.mite)$class
+mite.post <- predict(LDA.mite)$posterior
+mite.table <- table(mite.xy.group[,4], mite.class)
+diag(prop.table(mite.table, 1))
+</code>
-<code rsplus>
+</hidden>
-a <- list(1:100, 101:200)
-# appliquons mean a tous les elements
-lapply(a, mean)  # on obtient une liste en retour
-unlist(lapply(a, mean)) # utilisez unlist pour avoir un vecteur a la place
-vapply(a, mean, 0) # on dit donc a vapply que notre resultat sera un simple nombre
-</code>
-======= 4. Brève introduction à quelques paquets utiles dans R =======
+======5. Autres méthodes d'ordination utiles======
-===== Knitr =====
+<code rsplus | autres méthodes>
-Knitr est un paquet qui peut être utilisé pour générer des rapports de façon dynamique ou des pages web à partir de code en R. Le code est évalué au moment de créer le rapport.
+?cca # Analyse canonique des correspondances (Constrained Correspondence Analysis, CCA)
-Le code peut être écrit facilement en RStudion en utilisant le langage Markdown. :
-++++ Exemple de code Markdown |
-<code rsplus>
----
-title: "Interesting packages for R / QCBS Workshop on R programming"
-output: html_document
----
-### Data table
+# Méthode d’analyse canonique similaire à la RDA préservant les distance de Chi-carré entre objets (au lieu des
-A package to facilitate and to improve the efficiency of certain operations in R.
+# distances euclidiennes dans le cas d’une RDA). Cette méthode est mieux adaptée à l’étude de longs gradients
-```{r}
+# que la RDA.
-library(data.table)
-mydf=data.frame(a=rep(LETTERS,each=1e5),b=rnorm(26*1e5))
-mydt=data.table(mydf)
-setkey(mydt,a) # We set the column that will be used as a key for the data table
-```
-Returns all rows with column a (the key) equal to F
-````{r}
-mydt['F']
-````
-Gives the mean value of column b for each letter in column a.
-```{r}
-mydt[,mean(b),by=a]
-```
-methods/html/Compare.html">Compare
+?CCorA # Analyse canonique de corrélations (Canonical Correlation Analysis, CCorA)
-```{r}
-system.time(t1<-mydt[,mean(b),by=a])
-```
-### With tapply()
+# Cette analyse canonique diffère d’une RDA puisque les deux matrices étudiées sont considérées symétriques
-```{r}
+# tandis que dans une RDA la matrice Y dépend toujours de X. Cette méthode est principalement utilisée pour
-system.time(t2<-tapply(mydf$b,mydf$a,mean))
+# tester la significativité des corrélations entre deux jeux de données multivariés et explorer la structure des
-```
+# données présentes.
-### With reshape2
+help(coinertia, package=ade4) # Coinertia Analysis
-```{r message=FALSE}
-library(reshape2)
-meltdf=melt(mydf)
-system.time(t3<-dcast(meltdf,a~variable,mean))
-```
-### With plyr
+help(coinertia, package=ade4) # Analyse de co-inertie (Coinertia Analysis, CoIA)
-```{r}
-library(plyr)
-system.time(t4<-ddply(mydf,.(a),summarize,mean(b)))
-```
-### With dplyr
+# Méthode d’analyse canonique symétrique permet de comparer des jeux de données jouant des rôles
-```{r message=FALSE}
+# équivalents lors de l’analyse. Cette méthode calcule un espace d’ordination commun dans lequel les objets et les
-library(dplyr)
+# variables des deux jeux de données sont projetées et comparées. Par rapport à une CCorA, la CoIA n’impose pas
-ti1<-proc.time()
+# de contraintes vis-à-vis du nombre de variables des deux jeux de données, et permet donc de comparer des
-groups <- group_by(mydf, a)
+# 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
-t4b <- summarise(groups, total = mean(b))
+# données appariés.
-eltime<-proc.time()-ti1
-eltime
-```
-### With sqldf
-```{r}
-library(sqldf)
-system.time(t5<-sqldf('SELECT a, avg(b) FROM mydf GROUP BY a'))
-```
-### With a for loop
+help(mfa, package=ade4) # Analyse factorielle multiple (Multiple Factorial Analysis, MFA)
-```{r}
-ti1<-proc.time()
-t6<-data.frame(letter=unique(mydf$a),mean=rep(0,26))
-for (i in t6$letter ){
-  t6[t6$letter==i,2]=mean(mydf[mydf$a==i,2])
-}
-eltime<-proc.time()-ti1
-eltime
-```
+# 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.
-### With a parallelized FOR loop
-```{r}
-library(foreach)
-library(doMC)
-registerDoMC(4) #Four-core processor
-ti1<-proc.time()
-t7<-data.frame(letter=unique(mydf$a),mean=rep(0,26))
-t7[,2] <- foreach(i=t7$letter, .combine='c') %dopar% {
- mean(mydf[mydf$a==i,2])
-}
-eltime<-proc.time()-ti1
-eltime
-```
+# Les packages AEM et PCNM permettent d’effectuer diverses formes d’analyses spatiales :
+https://r-forge.r-project.org/R/?group_id=195
-### RgoogleMaps
-```{r message=FALSE}
-library(RgoogleMaps)
-myhome=getGeoCode('Olympic Stadium, Montreal');
-mymap<-GetMap(center=myhome, zoom=14)
-PlotOnStaticMap(mymap,lat=myhome['lat'],lon=myhome['lon'],cex=5,pch=10,lwd=3,col=c('red'));
-```
-### Taxize
-```{r message=FALSE, warning=FALSE}
-library(taxize)
-spp<-tax_name(query=c("american beaver"),get="species", db = 'ncbi')
-spp
-fam<-tax_name(query=c("american beaver"),get="family", db = 'ncbi')
-fam
-correctname <- tnrs(c("fraxinus americanus"))
-correctname
-cla<-classification("acer rubrum", db = 'itis')
-cla
-```
-### spocc
-```{r message=FALSE}
-library(spocc)
-occ_data <- occ(query = 'Acer nigrum', from = 'gbif')
-mapggplot(occ_data)
-```
-### Combine spocc and RgoogleMaps
-```{r message=FALSE, warning=FALSE}
-occ_data <- occ(query = 'Puma concolor', from = 'gbif')
-occ_data_df=occ2df(occ_data)
-occ_data_df<-subset(occ_data_df,!methods/html/is.html">is.na">is.na(latitude) & latitude!=0)
-mymap<-GetMap(center=c(mean(occ_data_df$latitude),mean(occ_data_df$longitude)), zoom=2)
-PlotOnStaticMap(mymap,lat=occ_data_df$latitude,lon=occ_data_df$longitude,cex=1,pch=16,lwd=3,col=c('red'));
-```
-### geonames
-```{r message=FALSE}
-library(geonames)
-options(geonamesUsername="glaroc")
-res<-GNsearch(q="Mont Saint-Hilaire")
-res[,c('toponymName','fclName')]
-dc<-GNcities(45.4, -73.55, 45.7, -73.6, lang = "en", maxRows = 10)
-dc[,c('toponymName')]
-```
 </code>
-++++
-[[http://qcbs.ca/wp-content/uploads/2014/12/QCBS_atelier8_knitr.html|Voir la page web résultante]].
-===== Data Table =====
+======Références======
-[[http://cran.r-project.org/web/packages/data.table/index.html|Data table]] est un paquet très utile qui peut faciliter et améliorer l'efficacité de certaines opérations en R. Les tables de données (data tables) sont très similaires aux data frames. Vous pouvez même les construire à partir de data frames.
-[[http://cran.r-project.org/web/packages/data.table/vignettes/datatable-intro.pdf|Introduction to Data table (PDF)]]
+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.
-<code rsplus>
+Borcard, Gillet & Legendre (2011). Numerical Ecology with R. Springer New York.
-install.packages('data.table')
-library(data.table)
-</code>
-Créez un très long jeu de données avec une colonne contenant des lettres et une colonne avec des nombres aléatoires.
+Breiman, L., J. H. Friedman, et al. (1984). Classification and Regression Trees. Belmont, California, USA, Wadsworth International Group.
-<file rsplus>
-mydf<-data.frame(a=rep(LETTERS,each=1e5),b=rnorm(26*1e5))
-</file>
-Convertir le data frame en table de données.
+Budyko, M.I. (1969) The effect of solar radiation variations on the climate of the Earth. Tellus, 21(5), 611-619.
-<file rsplus>
-mydt<-data.table(mydf)
-</file>
-Une clé doit être attribuée à chaque table de données. Cette clé doit être une (ou plus) colonne provenant de la table. Cette clé est à la base de l'organisation de la table de données.
+Clarke & Warwick (2001). Change in Marine Communities: An Approach to Statistical Analysis and Interpretation 2nd edition. Primer-E Ltd.
-<file rsplus>
-setkey(mydt,a)
-</file>
-Une fois que la clé est attribuée, nous pouvons facilement extraire, par exemple, toutes les lignes contentant la clé (colonne a) égale à F.
+De'ath, G. (2002). Multivariate regression trees : a new technique for modeling species-environment relationships. Ecology, 83(4), 1105–1117.
-<file rsplus>
-mydt['F']
-</file>
-Donne la valeur moyenne pour la colonne b, pour chaque lettre de la colonne a.
+Gotelli & Ellison (2004). A Primer of Ecological Statistics. Sinaeuer Associates Inc., Sunderland MA.
-<file rsplus>
-mydt[,mean(b),by=a]
-</file>
-Comparons maintenant la performance de Data table avec les autres méthodes nous permettant d'effectuer la même tâche.
+Legendre & Legendre (2012). Numerical Ecology 3rd edition. Elsevier Science BV, Amsterdam.
-<file rsplus>
-system.time(t1<-mydt[,mean(b),by=a])
-</file>
-**Avec tapply()**
+Poulin, Andersen & Rochefort (2013) A new approach for tracking vegetation change after restoration: a case study with peatlands. Restoration Ecology, 21, 363-371.
-<file rsplus>
-system.time(t2<-tapply(mydf$b,mydf$a,mean))
-</file>
-**Avec [[http://cran.r-project.org/web/packages/reshape2/index.html|reshape2]]**
-**NOTE**: plyr et reshape2 ont été traités dans l'[[r_atelier4|atelier 4]].
-<file rsplus>
-library(reshape2)
-meltdf<-melt(mydf)
-system.time(t3<-dcast(meltdf,a~variable,mean))
-</file>
-**Avec [[http://cran.r-project.org/web/packages/plyr/index.html|plyr]]**
-, un suite d'outils qui peuvent être utiliser pour séparer des données en blocs homogènes, appliquer une fonction sur chaque bloc, et remettre les blocs ensemble.
-<file rsplus>
-library(plyr)
-system.time(t4<-ddply(mydf,.(a),summarize,mean(b)))
-</file>
-**Avec [[http://cran.rstudio.com/web/packages/dplyr/vignettes/introduction.html|dplyr]]**
-, un nouvelle version de plyr qui est plus rapide et adaptée spécialement aux data frames.
-<file rsplus>
-library(dplyr)
-ti1<-proc.time()
-groups <- group_by(mydf, a)
-t5 <- summarise(groups, total = mean(b))
-eltime<-proc.time()-ti1
-</file>
-**Avec [[http://cran.r-project.org/web/packages/sqldf/index.html|sqldf]]**. Ce paquet permet d'écrire des requêtes de types SQL (Structured Query Language) sur des data frames.
-<file rsplus>
-library(sqldf)
-system.time(t6<-sqldf('SELECT a, avg(b) FROM mydf GROUP BY a'))
-</file>
-**Avec une boucle FOR**
-<file rsplus>
-ti1<-proc.time()
-# Initialisation d'un data frame vide avec deux colonnes et 26 lignes.
-t7<-data.frame(letter=unique(mydf$a),mean=rep(0,26))
-for (i in t6$letter ){
-  t7[t7$letter==i,2]=mean(mydf[mydf$a==i,2])
-}
-eltime<-proc.time()-ti1
-eltime
-</file>
-**Avec une boucle FOR parallèlisée**
-On pour utiliser les paquets [[http://cran.r-project.org/web/packages/foreach/index.html|foreach]] et [[http://cran.r-project.org/web/packages/doMC/index.html|doMC]] pour executer des sections de code en parallèle sur des ordinateurs avec plusieurs coeurs. C'est particulièrement utile pour accélérer des calculs dans des boucles FOR dans lesquelles chaque itération roule indépendemment des autres. Notez que doMC peut ne pas fonctionne sous Windows. Ça devrait cependant fonctionner sous Linux ou Mac OSX.
-<file rsplus>
-library(foreach)
-library(doMC)
-registerDoMC(4) #Processeur quatre-coeurs
-ti1<-proc.time()
-t8<-data.frame(letter=unique(mydf$a),mean=rep(0,26))
-t8[,2] <- foreach(i=t8$letter, .combine='c') %dopar% {
- mean(mydf[mydf$a==i,2])
-}
-eltime<-proc.time()-ti1
-eltime
-</file>
-===== RgoogleMaps! =====
-Le paquet [[http://cran.r-project.org/web/packages/RgoogleMaps/index.html|RgoogleMaps]] vous permet d'afficher très simplement des images de Google Maps ou Google Satellite  allows to very simply show Google maps or Google Satellite dans R, centrées sur la localisation de votre choix. Vous pouvez également superposer des données spatiales relativement aisément sur ces cates. La fonction getGeocode() permet de transformer une entrée de recherche pour un code postal ou un nom de lieu en coordonnées latitude, longitudes en utilisant les services Google.
-<file rsplus>
-library(RgoogleMaps)
-myhome=getGeoCode('Olympic stadium, Montreal');
-mymap<-GetMap(center=myhome, zoom=14)
-PlotOnStaticMap(mymap,lat=myhome['lat'],lon=myhome['lon'],cex=5,pch=10,lwd=3,col=c('red'));
-</file>
-===== Taxize =====
-Le projet [[http://ropensci.org/|rOpenSci]] supporte le développement d'un nombre important de paquets R pour faciliter l'accès à des sources de données en ligne. Parmi celles-ci, figure le paquet [[http://cran.r-project.org/web/packages/taxize/index.html|Taxize]], qui peut être utilisé pour extraire l'information taxonomique provenant de différents bases de données. On peut extraire, par exemple, des synonymes, des hiérarchies taxonomiques, les noms communs, et plus, de plus d'une dizaine de sources.
-<file rsplus>
-library(taxize)
-spp<-tax_name(query=c("american beaver"),get="species")
-fam<-tax_name(query=c("american beaver"),get="family")
-correctname <- tnrs(c("fraxinus americanus"))
-cla<-classification("acer rubrum", db = 'itis')
-</file>
-===== Spocc =====
-Un autre paquet de rOpenSci qui est très utile est [[http://cran.r-project.org/web/packages/spocc/index.html|Spocc]]. Il permet d'effectuer des recherches de données d'occurrence d'espèces provenant de plusieurs sources, dont [[http://www.gbif.org/|Global Biodiversity Information Facility]], une immense base de données mondiale contenant des centaines de millions d'occurrences venant de données de terrain ou de collections.
-<file rsplus>
-library(spocc)
-occ_data <- occ(query = 'Acer nigrum', from = 'gbif')
-mapggplot(occ_data)
-</file>
-Combinez spocc et RgoogleMaps
-<file rsplus>
-occ_data <- occ(query = 'Puma concolor', from = 'gbif')
-occ_data_df=occ2df(occ_data)
-occ_data_df<-subset(occ_data_df,!is.na(latitude) & latitude!=0)
-mymap<-GetMap(center=c(mean(occ_data_df$latitude),mean(occ_data_df$longitude)), zoom=2)
-PlotOnStaticMap(mymap,lat=occ_data_df$latitude,lon=occ_data_df$longitude,cex=1,pch=16,lwd=3,col=c('red'));
-</file>
-===== geonames =====
-[[http://geonames.r-forge.r-project.org/|Geonames]] connecte R à [[http://geonames.org|Geonames.org]], une base de données de noms de lieux et de toponymes.
-<file rsplus>
-library(geonames)
-options(geonamesUsername="glaroc")
-# Trouver les noms de lieux qui contiennent le terme "Mont Saint-Hilaire"
-res<-GNsearch(q="Mont Saint-Hilaire")
-res[,c('toponymName','fclName')]
-#Extraire toutes les villes dans un rectangle définit par deux coins géographiques.
-dc<-GNcities(45.4, -73.55, 45.7, -73.6, lang = "en", maxRows = 10)
-dc[,c('toponymName')]
-</file>