Differences

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

Link to this comparison view

Both sides previous revision Previous revision
Next revision
Previous revision
r_atelier8 [2016/03/18 08:00]
emmanuelle.chretien [3.5. En résumé]
r_atelier8 [2021/10/13 23:52] (current)
lsherin
Line 1: Line 1:
-======= Ateliers R du CSBQ =======+<WRAP group> 
 +<WRAP centeralign>​ 
 +<WRAP important>​ 
 +<wrap em> __AVIS IMPORTANT__ </​wrap> ​
  
-[[http://​qcbs.ca/fr/​|{{:​logo_text.png?​nolink&​500|}}]]+<wrap em> Depuis l'​automne 2021, ce wiki a été discontinué et n'est plus activement développé</wrap>
  
-Cette série ​de [[r|10 ateliers]] ​guide les participants ​à travers les étapes requises ​afin de maîtriser le logiciel R pour une grande variété ​d’analyses statistiques pertinentes en recherche en biologie ​et en écologie. Ces ateliers en libre accès ont été créés par des membres du CSBQ à la fois pour les membres du CSBQ et pour la grande communauté d’utilisateurs de R.+<wrap em> Tout le matériel mis à jour et les annonces pour la série ​d'​ateliers R du CSBQ se trouvent maintenant sur le [[https://r.qcbs.ca/​fr/​workshops/​r-workshop-08/​|site web de la série d'ateliers ​R du CSBQ]]. Veuillez mettre ​à jour vos signets en conséquence ​afin d'​éviter les documents périmés ​et/ou les liens brisés</​wrap>​
  
-====== Atelier 8: Analyses multivariées ======+<wrap em> Merci de votre compréhension,​ </​wrap>​
  
-Développé par : Bérenger Bourgeois, Xavier Giroux-Bougard,​ Amanda Winegardner,​ Emmanuelle Chrétien et Monica Granados. +<wrap em> Vos coordonnateurs de la série d’ateliers ​du CSBQ</wrap>
-(Le script est en parti issu de: Borcard, Gillet & Legendre (2011). //Numerical Ecology with R//. Springer New York). ​+
  
-**Résumé:​** Apprenez les bases des analyses multivariées : choisissez les mesures de distance et transformations appropriées pour réaliser des analyses multivariées,​ et apprenez à effectuer des groupements,​ des Analyses en Composantes Principales,​ des Analyses de Correspondance,​ des Analyses en Coordonnées Principales et des Positionnements Multidimensionnels Non-métriques,​ afin de trouver les patrons de diversité des communautés biologiques !+</​WRAP>​ 
 +</​WRAP>​ 
 +<WRAP clear></​WRAP>​
  
-Lien vers la présentation Prezi associée : [[https://​prezi.com/​puivxp0qd-tz/​csbq-atelier-r-8-good-version/​| Prezi]]+======= Ateliers R du CSBQ =======
  
 +[[http://​qcbs.ca/​fr/​|{{:​logo_text.png?​nolink&​500|}}]]
  
-**Téléchargez le script R et les données pour cet atelier:​** +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 écologieCes 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.
-  * [[http://​qcbs.ca/​wiki/​_media/​code_atelier8.r| script R]] +
-  * [[http://​qcbs.ca/​wiki/​_media/​DoubsEnv.csv|DoubsEnv]] +
-  * [[http://​qcbs.ca/​wiki/​_media/​DoubsSpe.csv|DoubsSpe]] +
-  * [[http://​qcbs.ca/​wiki/​_media/​coldiss.R|Coldiss (fonction R)]]+
  
-**Télechargez les paquets R pour cet atelier:** +//Le contenu de cet atelier ​a été révisé par plusieurs membres du CSBQSi vous souhaitez y apporter des modifications,​ veuillez SVP contacter les coordonnateurs actuels de la série, listés sur la page d'​accueil//
-  * [[http://​cran.r-project.org/web/packages/​vegan/​index.html|vegan]] +
-  * [[http://​cran.r-project.org/​web/​packages/​gclus/​index.html|gclus]] +
-  * [[http://​cran.r-project.org/​web/​packages/​ape/​index.html|ape]]+
  
-<code rsplus | Chargez les paquets et les fonctions nécessaires> +<wrap em>AVIS IMPORTANT: MISES À JOUR MAJEURES</wrap>
-install.packages("​vegan"​) +
-install.packages("​gclus"​) +
-install.packages("​ape"​) +
-library(vegan) +
-library(gclus) +
-library(ape) +
-source(file.choose()) #coldiss.R  +
-</code> +
-======Qu'​est-ce que l'​ordination?​======+
  
-L'ordination ​est un ensemble de méthodes pour décrire des échantillons dans de multiples dimensions (Clarke et Warwick 2011)Les méthodes d'​ordination sont donc très utiles ​pour simplifier et interpréter des données multivariées. Les écologistes parlent souvent de "faire une PCA" face à des données multidimensionnelles complexes et désordonnées. Programmer des méthodes ​d'ordination à partir de R est relativement simpleL'​interprétation ​des analyses d'​ordination peut par contre être plus difficile, surtout si vous n'​êtes pas sûr des questions biologiques que vous souhaitez explorer avec la méthode d'​ordination que vous utilisezUn examen attentif des objectifs de ces méthodes et de leur cadre d'​application est nécessaire pour obtenir de bons résultats !+**Mise à jour de mars 2021:** Ce wiki n'​est ​plus activement développé ou mis à jourLe matériel mis à jour pour la série ​d'ateliers ​du CSBQ est maintenant hébergé sur [[https://​github.com/​QCBSRworkshops/​workshop08|la page GitHub]] ​des ateliers R du CSBQ.
  
-Lorsque vous utilisez une méthode d'​ordination,​ un ensemble de variables est utilisé ​pour ordonner des échantillons (objets, sites, etc.) le long d'axes principaux représentant des combinaison de variables (Gotelli et Ellison 2004)L'​ordination permet donc de réduire ou de simplifier les données en créant de nouveaux axes intégrant la majeure partie de la variation présentes dans les donnéesA titre d'​exemple,​ un ensemble de données avec 24 variables peut être réduit à cinq composantes principales ​qui représentent les principaux gradients de variation entre les échantillonsLes méthodes d'​ordination sans contrainte ne sont pas adaptées au test d'​hypothèses biologiques,​ mais permettent l'​analyse exploratoire des données. Voir ([[http://ordination.okstate.edu/|the Ordination Website]]) pour un aperçu des différents méthodes (en Anglais)+Le matériel disponible inclut; 
 +  - La [[https://​qcbsrworkshops.github.io/​workshop08/​pres-fr/​workshop08-pres-fr.html|présentation Rmarkdown]] ​pour cet atelier; 
 +  - Un [[https://​qcbsrworkshops.github.io/​workshop08/​book-fr/​workshop08-script-fr.R|script R]] qui suit la présentation - //*en construction*//​. 
 +  - [[https://qcbsrworkshops.github.io/workshop08/​book-fr/​index.html|Le matériel écrit]] qui accompagne la présentation ​en format bookdown - //*en construction*//​.
  
-======Introduction aux données======+====== ​Atelier 8 : Modèles additifs généralisés ​======
  
-Nous allons utiliser deux principaux ensembles de données dans la première partie de cet atelier. "​DoubsSpe.csv"​ est une matrice de données d'​abondance d'​espèces de communautés de poissons dans laquelle la première colonne contient les noms des sites de 1 à 30 et les colonnes subséquentes correspondent aux différentes espèces de poissons. "​DoubsEnv.csv"​ est une matrice de données environnementales pour les mêmes sites. La première colonne contient donc les noms des sites de 1 à 30 et les colonnes suivantes les mesures de 11 variables abiotiques. Notez que les données utilisées pour les analyses d'​ordination sont généralement en  [[http://​en.wikipedia.org/​wiki/​Wide_and_narrow_data|format long (en anglais)]].+Développé par : Eric Pedersen ​et Zofia Taranu
  
-<code rsplus | Chargez DoubsSpe et DoubsEnv>​ +Révision par Cédric Frenette Dussault
-#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écutez cette ligne une seule fois. +
  
-#L'environnement:​ “DoubsEnv.csv” +**Résumé:​** ​L'objectif de l'​atelier d'​aujourd'​hui sera d'​examiner ce que nous entendons par un modèle non-linéaire et comment les GAMs (modèles additifs généralisésnous permettent de modéliser les relations non-linéairesNous examinerons également comment tracer et interpréter ces relations non-linéairescomment ajouter des interactions,​ comment prendre en compte ​la non-indépendance des données (//e.g.// erreurs autocorrélées) et comment inclure des effets aléatoires en se basant sur les ateliers précédents. Enfin, nous allons brièvement aborder la mécanique derrière le fonctionnement des GAMs.
-env<read.csv(file.choose(), row.names=1) +
-env<env[-8,] #Supprimer site 8 puisqu'​on l'a retiré de la matrice d'​abondances +
-#Exécutez cette ligne une seule fois  +
-</code>+
  
-======1Exploration des données======+**Lien vers la nouvelle [[https://​qcbsrworkshops.github.io/​workshop08/​workshop08-fr/​workshop08-fr.html|présentation Rmarkdown]]** ​
  
-=====1.1 Données sur les espèces=====+Lien vers la présentation Prezi associée : [[https://​prezi.com/​ddkjfsxlr-za/​| Prezi]]
  
-Nous pouvons utiliser les fonctions de résumé ​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).+Télécharger le script ​et les données ​pour cette leçon:  
 +  - [[http://​qcbs.ca/​wiki/​_media/​gam_e.r|Script]] 
 +  - [[http://​qcbs.ca/​wiki/​_media/​other_dist.csv|Données]]
  
-<code rsplus | Explorez DoubsSpe>​ +===== Objectifs de l'​atelier ===== 
-names(spe) # Les noms des colonnes ​ +  - Utiliser la librairie //mgcv// pour modéliser les relations non linéaires 
-dim(spe# Le nombre de lignes ​et de colonnes. ​ +  - Évaluer la sortie d'un GAM afin de mieux comprendre nos données 
-str(spe) # La structure ​interne de la matrice. ​ +  - Utiliser ​des tests pour déterminer si nos relations correspondent à des modèles non linéaires ou linéaires 
-head(spe# Les premières lignes.  +  - Ajouter des interactions non linéaires entre les variables explicatives 
-summary(spe# Les statistiques descriptives +  - Comprendre l'​idée d'une fonction de base (//basis function//) et la raison pour laquelle ça rend les GAMs si puissants ! 
-</​code>​+  - Comment modéliser la dépendance dans les données ​(autocorrélation, ​structure ​hiérarchique) en utilisant les GAMMs 
 +**Prérequis:​** Expérience du logiciel R (assez pour être en mesure d'​exécuter un script et d'​examiner les données et les objets dans Ret une connaissance de base de la régression simple ​(vous devez savoir ce qu'on entend par une régression linéaire et une ANOVA).
  
-Regardez la distribution des espèces+===== 0Le modèle linéaire ... et où il échoue =====
  
-<code rsplus | La distribution ​des espèces (DoubsSpe)>​ +Qu'​est-ce qu'un modèle linéaire ? La régression linéaire est ce que la plupart ​des gens apprennent avant tout en statistiques et est parmi les méthodes les plus performantes. Elle nous permet de modéliser une variable réponse en fonction de facteurs prédictifs et d'une erreur résiduelle. Tel que vu dans l'​atelier sur les [[http://​qcbs.ca/​wiki/​r_atelier4|modèles linéaires]],​ la régression fait cependant quatre suppositions importantes : 
-(ab<-table(unlist(spe))) #Les parenthèses signifient que la sortie s'​affiche immédiatement+  l'​erreur est distribuée normalement 
 +  - la variance des erreurs est constante 
 +  - chaque erreur est indépendante des autres ​(homoscédasticité) 
 +  - la réponse est linéaire: <m> {y} = {β_0} + {β_1}{x} </m>
  
-barplot(ab, las=1, xlab=”Abundance class”, ylab=”Frequency”,​ col=grey(5:0/5)) +Il n'y a qu'une façon pour qu'un modèle linéaire soit correctement appliqué ​:
-</​code>​+
  
-{{:spe_barplot.png?300|}} +{{::graphic0.1.png?350|}}
-Pouvez-vous voir qu'il y a une grande fréquence de zéros dans les données d'​abondance ?+
  
-Calculez le nombre d'absences. +et pourtant tant de façons pour qu'il ne le soit pas :
-<code rsplus | Absences>​ +
-sum(spe==0)  +
-</​code>​+
  
-Regardez la proportion de zéros dans les données de la communauté de poissons +{{::​graphic0.2.png|}}
-<code rsplus | Proportion de zéros>​ +
-sum(spe==0)/​(nrow(spe)*ncol(spe)) +
-</​code>​ +
-La proportion de zéros dans la matrice est de ~0.5+
  
-Calculez ​le nombre ​de sites où chaque espèce est présente. +Alors, comment résoudre ce problème ? Nous devons premièrement savoir ce que le modèle ​de régression cherche à faire : 
-<code rsplus | Nombre de sites où chaque espèce est présente>​ +   * ajuster une ligne qui passe au milieu des données, 
-spe.pres<​colSums(spe>​0) # Somme des sites où chaque ​espèce est présente.  +   * faire cela sans sur-ajuster les données, c'​est-à-dire en passant une ligne entre chaque ​point
-hist(spe.presmain=”Cooccurrence des espèces”las=1, xlab=”Fréquence”,​ breaks=seq(0,30, by=5), col=”grey”) +Les modèles linéaires le font en trouvant la meilleure ligne droite qui passe à travers les donnéesEn revancheles modèles additifs font cela en ajustant une courbe à travers les donnéesmais tout en contrôlant le degré de courbure de la ligne //(plus d'information sur cela plus bas)//.
-</​code>​ +
-Le plus grand nombre ​d'espèces se retrouvent dans un nombre intermédiaire de sites.+
  
-Calculez le nombre d'​espèces présentes à chaque site. Ici, nous utilisons une façon simpliste de calculer la richesse en espèces. Dans certains cas, le nombre d'​espèces présentes dans un site varie entre les sites car il peut y avoir une relation entre le nombre d'​individus comptés (l'​abondance) dans ce site et la richesse en espèces. [[http://​cc.oulu.fi/​~jarioksa/​softhelp/​vegan/​html/​diversity.html|richesse d'​espèces raréfiée (en anglais)]] est souvent une mesure plus appropriée que la richesse totale. La fonction rarefy () de vegan peut être utilisée pour calculer la richesse raréfiée des espèces.+===== 1Introduction aux GAMs =====
  
-<code rsplus | Richesse en espèces>​ +Examinons un exemplePremièrementnous allons générer ​des données et les représenter graphiquement.
-site.pres<- rowSums(spe>​0) # Nombre d'​espèces présentes dans chaque site +
-hist(site.presmain=”Richesse en espèces”,​ las=1, xlab=”Fréquence ​des sites”, ylab=”Nombre d'​espèces”,​ breaks=seq(0,​30,​ by=5), col=”grey”)  +
-</​code>​+
  
-=====1.Données sur l'​environnement=====+<file rsplus | Générer et tracer des données>​ 
 +library(ggplot2) 
 +set.seed(10) 
 +250 
 +runif(n,​0,​5) 
 +y_model ​3*x/(1+2*x) 
 +y_obs rnorm(n,​y_model,​0.1) 
 +data_plot ​qplot(x, y_obs) +  
 +            geom_line(aes(y=y_model)) +  
 +            theme_bw() 
 +print(data_plot) 
 +</​file>​
  
-Explorez les données environnementales pour détecter les colinéarités ​:+{{::​graphic1.1.jpg?​350|}}
  
-<code rsplus | Exploration de DoubsEnv>​ +Si nous modélisions cette relation par une régression linéaire, les résultats ne respecteraient pas les suppositions énumérées ci-dessus. Commençons par modéliser une régression en utilisant la méthode des moindres carrés en utilisant la fonction ''​gam()''​ de la librairie //mgcv// - donc en tant que modèle linéaire //(nous verrons plus bas comment utiliser cette fonction pour spécifier un terme non linéaire).//
-names(env) +
-dim(env) +
-str(env) +
-head(env) +
-summary(env) +
-pairs(env, main="​Données environnementales"​ )  +
-</code>+
  
-Dans ce casles 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 variablespeut être effectuée en utilisant la fonction decostand ​() de vegan. ​+<file rsplus | Modèle linéaire>​ 
 +library(mgcv) 
 +linear_model = gam(y_obs~x) 
 +model_summary=summary(linear_model) 
 +print(model_summary) 
 +data_plot = data_plot+ 
 +             ​geom_line(colour="​red"​, 
 +             aes(y=fitted(linear_model))) 
 +print(data_plot) 
 +</​file>​
  
-<code rsplus | Standardisation des données> +Nous pouvons constater à partir du sommaire que notre modèle linéaire explique une grande partie de la variance (R<sup>​2</​sup><sub>(adj)</​sub> ​ 0.639). Toutefoisles graphiques de diagnostic des résidus du modèle montrent que l'​écart type ne suit pas une distribution normale et que la variance n'est pas homoscédastique. De plusil reste un patron non-linéaire important. Essayons maintenant de résoudre ce problème en ajustant les données ​avec un terme non linéaire
-env.z<-decostand(env, method="​standardize"​) +
-apply(env.z2mean) # Les données ​sont maintenant centrées (moyennes~0)... +
-apply(env.z,​ 2, sd)   # et réduites (écart-type=1) +
-</​code>​+
  
-======2. Mesures ​d'association====== ​+Nous reviendrons sur ceci un peu plus tard, mais brièvement,​ les GAMs sont une forme non paramétrique de la régression où le <​m>​beta</​m>​*x<​sub>​i</​sub> ​d'une régression linéaire est remplacé par une fonction de lissage des variables explicatives,​ f (x<​sub>​i</​sub>​),​ et le modèle devient :
  
-L'​algèbre matricielle est à la base des méthodes d'​ordination. Une matrice est constituée de données ​(ex. valeurs mesuréesréparties en lignes et colonnes. Les analyses d'​ordination sont effectuées sur des matrices d'​association calculées à partir des matrices de données écologiques (telles que DoubsEnv ou DoubsSpe). La création d'une matrice d'​association permet le calcul de similarité et de distance entre les objets ou les descripteurs (Legendre et Legendre 2012). Avant de se lancer dans des analyses d'​ordination,​ il est important de passer du temps sur vos matrices de données. Explorer les mesures possibles d'​association qui peuvent être générées à partir de vos données avant de faire une ordination peut vous aider à mieux comprendre quelles mesures de distance sont appropriées pour les méthodes d'​ordination. Il peut être difficile de voir l'​objectif de chaque indice de dissimilarité,​ mais cette connaissance sera nécessaire pour mieux comprendre les méthodes d'​ordination canoniques présentées par la suite.+<​m>​{y_i} = f(x_{i}+ {ε_i}</​m>​
  
-EN RÉSUMÉ: Pour l'​ordination d'​objetsil faut calculer les distances entre eux. Ces distances peuvent être calculées de plusieurs façons, en prenant en compte l'​abondance ou les données de présence ​absence. Plus important encore, plusieurs propriétés sont de grande importance pour les mesures de distances, et elles seront explorées dans les exemples ci-dessous. Pour plus d'​informations sur les propriétés des mesures de distance et certains termes clés, cliquez sur la section cachée.+où y<​sub>​i</​sub>​ est la variable réponsex<​sub>​i<​/sub> est la covariable, et //f// est la fonction lissage.
  
-<hidden> +Étant donné que la fonction de lissage f(x<sub>i</​sub>​) est non linéaire et locale, l'​ampleur de l'​effet de la variable explicative peut varier en fonction de la relation entre la variable et la réponse. Autrement dit, contrairement à un coefficient fixe <​m>​beta</​m>,​ la fonction //f// peut changer tout au long du gradient x<​sub>​i</​sub>​. Le degré de lissage de //f// est contrôlée en utilisant une régression pénalisée qui est déterminée automatiquement à l'aide d'une méthode de validation croisée généralisée (GCV) de la librairie //mgcv// (Wood 2006).
-**Termes clés:​** ​+
  
-**Association -** «Terme général pour décrire toute mesure ou coefficient servant à quantifier la ressemblance ou différence entre les objets ou les descripteurs. Dans une analyse entre descripteurs,​ zéro signifie l'absence d'association.» ​(Legendre et Legendre 2012). +Avec ''​gam()'' ​les termes non linéaires sont spécifiés par des expressions de la forme: ​''​s(x)''​.
  
-**Similarité -** une mesure dont «le maximum ​(1est atteint lorsque deux objets sont identiques et le minimum lorsque deux objets sont complètement différents.» ​(Legendre et Legendre 2012)+<file rsplus | GAM> 
 +gam_model = gam(y_obs~s(x)) 
 +summary(gam_model) 
 +data_plot ​data_plot +   
 +     ​geom_line(colour="​blue",​aes(y=fitted(gam_model))) 
 +print(data_plot) 
 +</​file>​
  
-**Distance ​(aussi «dissimilarité»-** une mesure dont «le maximum (D=1est atteint lorsque deux objects sont complètement différent.» ​(Legendre et Legendre 2012). Distance ou dissimilarité ​(D= 1-S+La variance expliquée par notre modèle a augmenté de 20% (R<​sup>​2</​sup><​sub>​(adj)</​sub> ​0.859et quand on compare l'​ajustement du modèle linéaire ​(rougeau modèle non-linéaire ​(bleu), il est clair que l'​ajustement de ce dernier est relativement meilleur.
  
-Le choix d'une mesure d'​association dépend de vos données, mais aussi de ce que vous en savez d'un point de vue écologiquePar exemple, la distance euclidienne est une mesure de distance très commune, facile à utiliser et utile pour comprendre comment les différences entre deux échantillons sont basées sur la cooccurrence des espèces. Le calcul de la distance euclidienne prend en compte les zéros dans les données, ce qui signifie que deux échantillons ou sites ne contenant aucune espèce en commun (double-absence) peuvent sembler plus similaires que deux sites partageant quelques espèces. Dans ce cas, la distance euclidienne peut être trompeuse et il est souvent préférable de choisir une mesure de distance différente si beaucoup d'​espèces ont une abondance nulle dans votre matrice. Cette propriété est communément appelée le problème des «double zéros» en ordination.+{{::​graphic1.3.jpg?350|}}
  
-Quelques mesures de distance (d'après Gotelli et Ellison 2004):+La librairie //mgcv// comprend également une fonction ​''​plot''​ qui, par défaut, nous permet de visualiser la non-linéarité du modèle.
  
-^Mesure^Propriété^Description^ +<file rsplus | smooth plot> 
-^Euclidienne^Métrique^Distance entre deux points dans un espace en 2D.^ +plot(gam_model
-^Manhattan^Métrique^Distance entre deux points - la distance est la somme des différences entre coordonnées cartésiennes.^ +</file>
-^Corde^Métrique^Généralement utilisée pour déterminer les différences dues à la dérive génétique.^ +
-^Mahalanobis^Métrique^Distance entre un point et une distribution,​ où la distance est le nombre d'​écart-types du point correspondant à la moyenne de la distribution.^ +
-^Chi-carré^Métrique^Similaire à la distance euclidienne.^ +
-^Bray-Curtis^Semi-métrique^Dissimilarité entre deux échantillons ​(ou sitesoù la somme des valeurs minimales des espèces présentes dans les deux sites sont divisées par la somme des espèces répertoriées dans chaque site.^ +
-^Jaccard^Métrique^[[http://​en.wikipedia.org/​wiki/​Jaccard_index|Description]] ^ +
-^Sorensen'​s^Semi-métrique^Bray-Curtis correspond à 1 - Sorensen^ +
-</hidden>+
  
-=====2.1 Mesures ​de distance=====+Comment utilisons-nous les GAMs pour savoir si un modèle linéaire est suffisant pour modéliser nos données ? Nous pouvons utiliser les fonctions ''​gam()''​ et ''​anova()''​ pour tester formellement si une hypothèse de linéarité est justifiéeNous devons simplement le configurer ​de sorte que notre modèle non-linéaire est emboîté dans notre modèle linéaire; c'est à dire, nous devons créer un objet qui inclut à la fois ''​x''​ (linéaire) et ''​s(x)''​ (non-linéaire). En utilisant la fonction ''​anova()'',​ on vérifie si l'​ajout de ''​s(x)''​ au modèle avec seulement ''​x''​ comme covariable est justifié par les données.
  
-**Les données quantitatives des espèces** +<file rsplus | Test de linéarité>​ 
-Nous pouvons utiliser la fonction vegdist ​() pour calculer des indices de dissimilarité sur des données de composition de la communauté. Ceux-ci peuvent ensuite être visualisés sous forme de matrice si désiré.+linear_model = gam(y_obs~x) 
 +nested_gam_model = gam(y_obs~s(x)+x) 
 +print(anova(linear_model,​ nested_gam_model,​ test="​Chisq"​)) 
 +</​file>​
  
-<code rsplus | vegdist>​ +Le terme non linéaire est significatif:​
-spe.db<​-vegdist(spe,​ method="​bray"​) # distance de Bray (avec des données de présence-absence,​ correspond à Sorensen) ​  +
-spe.dj<​-vegdist(spe,​ method="​jac"​) # distance de Jaccard  +
-spe.dg<​-vegdist(spe,​ method="​gower"​) # distance de Gower   +
-spe.db<​-as.matrix(spe.db) # réarranger en format matrice (pour visualisation,​ ou pour exporter en .csv) +
-</​code>​+
  
-Une version condensée de la matrice spe.db représentant la distance entre les trois premières espèces de DoubsSpe ressemblerait à ceci: +   ​Analysis of Deviance Table 
-^ ^Espèce ​1^Espèce ​2^Espèce 3^ +   Model 1: y_obs ~ x 
-^Espèce 1^0.0^0.6^0.68^ +   ​Model ​2: y_obs ~ s(x) + x 
-^Espèce ​2^0.6^0.0^0.14^ +         ResidDf    ResidDev     ​Df ​       Deviance ​   Pr(>​Chi) ​    
-^Espèce 3^0.68^0.14^0.0+   ​1 ​   248.00 ​       ​6.5846                              ​ 
-Vous pouvez voir que lorsque l'on compare une espèce à elle-même (par exemple, Espèce ​à Espèce ​1), la distance est de 0 puisque les espèces sont identiques.+   ​   240.68        2.4988         7.3168    4.0858 ​     < 2.2e-16 *** 
 +   --- 
 +   ​Signif. codes: ​ 0 ‘***’ ​0.001 ‘**’ ​0.01 ‘*’ ​0.05 ‘.’ ​0.‘ ’ 1
  
-Ces mêmes mesures de distance peuvent être calculées à partir des données de présence-absence par l'​utilisation de l'​argument binary=TRUE dans la fonction vegdist(). Cela donnera des mesures de distance légèrement différentes.+----
  
-Vous pouvez également créer des représentations graphiques de ces matrices d'​association en utilisant la fonction coldiss. 
-<​hidden>​ 
-Cette fonction peut être sourcée à partir du script coldiss.R: 
-<code rsplus | coldiss> 
-windows() 
-coldiss(spe.db,​ byrank=FALSE,​ diag=TRUE) # Carte des points chauds Bray-Curtis ​ 
-windows() 
-coldiss(spe.dj,​ byrank=FALSE,​ diag=TRUE) # Carte des points chauds Jaccard ​ 
-windows() ​ 
-coldiss(spe.dg,​ byrank=FALSE,​ diag=TRUE) # Carte des points chauds Gower  
-</​code>​ 
-{{:​coldiss_Bray.png?​800|}} 
  
-La figure montre une matrice de dissimilarité dont les couleurs reflètent la mesure de distance. La couleur violet est associée aux zones de fortes dissimilarités. +**DÉFI 1**
-</​hidden>​+
  
-**Données environnementales quantitatives** +Nous allons maintenant essayer cela avec d'​autres données générées aléatoirement. Nous allons d'​abord générer ​les données. Ensuite, nous allons ajuster un modèle linéaire et un GAM à la relation ​entre ''​x_test''​ et ''​y_test_obs''​. 
-Regardons //les associations// ​entre les variables environnementales (aussi appelée mode Q): +Quels sont les degrés ​de libertés effectifs du terme non-linéaire ? Déterminez si l'​hypothèse ​de linéarité est justifiée pour ces données.
-<code rsplus | distances mesurées à partir des données environnementales>​ +
-env.de<-dist(env.z, method = "​euclidean"​) # matrice ​de distances euclidiennes des données ​envstandardisées +
-windows() # crée une nouvelle fenêtre graphique +
-coldiss(env.de,​ diag=TRUE) +
-</​code>​+
  
-Nous pouvons ensuite regarder //la dépendance//​ entre les variables environnementales (aussi appelée mode R): +<file rsplus| ​Générer des données>​ 
-<code rsplus | corrélations entre les données ​environnementales+<- 250 
-(env.pearson<-cor(env)) # coefficient r de corrélation de Pearson +x_test <- runif(n,-5,5
-round(env.pearson2 # arrondit les coefficients à deux décimales +y_test_fit ​<- 4*dnorm(x_test
-(env.ken<-cor(env, method="​kendall"​)) # coefficient tau de corrélation de rang de Kendall ​ +y_test_obs <- rnorm(n,y_test_fit, 0.2) 
-round(env.ken, 2)  +</file>
-</code>+
  
-La corrélation de Pearson mesure la corrélation linéaire entre deux variables. La corrélation de Kendall est une corrélation de rang qui quantifie la relation entre deux descripteurs ou deux variables lorsque les données sont ordonnées ​au sein de chaque variable.+++++ Réponse ​au défi 1 | 
  
-Dans certains cas, il peut y avoir des types mixtes de variables environnementales. Le mode Q peut alors être utilisé pour trouver des associations entre variables environnementales. C'est ce que nous allons faire avec l'​exemple fictif suivant:  +<file rsplus| ​Réponse au défi 1
-<code rsplus | Exemple+data_plot ​<- qplot(x_testy_test_obs 
-var.g1<-rnorm(300, 1+  ​geom_line(aes(y=y_test_fit))+ 
-var.g2<​-runif(30, 0, 5) +  ​theme_bw() 
-var.g3<​-gl(3, 10) +print(data_plot)
-var.g4<​-gl(2,​ 5, 30+
-(dat2<​-data.frame(var.g1,​ var.g2, var.g3, var.g4)+
-str(dat2) +
-summary(dat2) +
-</​code>​+
  
-Une matrice de dissimilarité peut être générée pour ces variables mixtes en utilisant la distance de Gower: +linear_model_test ​<- gam(y_test_obs~x_test) 
-<code rsplus | daisy> +nested_gam_model_test ​<- gam(y_test_obs~s(x_test)+x_test) 
-?daisy #Cette fonction peut gérer la présence de NA dans les données +print(anova(linear_model_test,​ nested_gam_model_testtest="Chisq"))
-(dat2.dg<-daisy(dat2metric="gower")) +
-coldiss(dat2.dg) +
-</​code>​+
  
-**Défi 1 - Niveau intermédiaire** +summary(nested_gam_model_test)$s.table 
-Discutez avec votre voisin: Comment pouvons-nous dire si des objets sont similaires avec un jeu de données multivariées?​ Faites une liste de toutes vos suggestions.+</​file>​
  
-**Défi ​- Solution** +  Analysis of Deviance Table 
-<hidden+   
-Discussion avec le groupe.  +   
-</hidden>+  Model 1: y_test_obs ~ x_test 
 +  Model 2: y_test_obs ~ s(x_test) + x_test 
 +       ​Resid. Df   ​Resid. Dev    Df      Deviance ​  Pr(>Chi)    ​ 
 +  ​1 ​    248.0      81.09                              ​ 
 +  ​2 ​    ​240.5 ​     7.46          7.5012 ​  ​73.629 ​   ​2.2e-16 *** 
 +  --- 
 +  Signif. codes: ​ 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
-**Défi 1 Niveau avancé** +                 ​edf ​  ​Ref.df ​       F p-value 
-Calculer à la mitaine // sans utiliser la fonction decostand ​() // les distances de Bray-Curtis et de Gower pour l'​abondance des espèces CHA, TRU et VAI dans les sites 1, 2 et 3.+  s(x_test7.602145 8.029057 294.0944       0
  
-**Défi 1 Solution** +Réponse: Oui la non-linéarité est justifiée. Les degrés de libertés effectifs (edf) sont >> 1.
-<hidden>+
  
-Formule pour calculer la distance de Bray-Curtisd[jk] = (sum abs(x[ij]-x[ik]))/​(sum (x[ij]+x[ik]))+{{::​challenge_1_soln.jpg?​350|}}
  
-Réduisez le jeu de données aux espèces CHA, TRU et VAI et aux sites 1, et 3 +++++ 
-<code rsplus | Subset>​ +===== 2. Plusieurs termes non linéaires =====
-spe.challenge<​-spe[1:​3,​1:​3] # les 3 premières lignes et 3 premières espèces (colonnes) +
-</​code>​+
  
-Déterminer l'​abondance totale des espèces pour chaque site d'intérêt (somme ​des trois lignes) qui correspondra au dénominateur de la distance de Bray-Curtis. +Avec les GAMs, il est facile ​d'ajouter ​des termes non linéaires et linéaires dans un seul modèleplusieurs termes non linéaires ou même des interactions non linéairesDans cette sectionnous allons utiliser un ensemble de données générées automatiquement par //mgcv//.
-<code rsplus | Abondance par site> +
->​(Abund.s1<​-sum(spe.challenge[1,])) +
-(Abund.s2<​-sum(spe.challenge[2,])) +
-(Abund.s3<​-sum(spe.challenge[3,​])) +
-</​code>​+
  
-Maintenant, calculez la différence de l'​abondance ​des espèces pour chaque paire de sites. Par exemple, quelle est la différence entre l'​abondance de CHA et TRU dans le site 1? Vous devez calculer les différences suivantes: +<file rsplus | Générer ​des données non-linéaires>​ 
-CHA et TRU site 1 +gam_data = gamSim(eg=5) 
-CHA et VAI site 1 +head(gam_data) ​ 
-TRU et VAI site 1 +</​file>​
-CHA et TRU site 2 +
-CHA et VAI site 2 +
-TRU et VAI site 2 +
-CHA et TRU site 3 +
-CHA et VAI site 3 +
-TRU et VAI site 3+
  
-<code rsplus | Différences d'​abondance>​ +Nous allons voir comment nous pouvons prédire la variable réponse y en fonction des autres variablesCommençons par un modèle de base comprenant un terme non linéaire ​(x1et un facteur qualitatif ​(X0 avec 4 niveaux).
-Spec.s1s2<​-0 +
-Spec.s1s3<​-0 +
-Spec.s2s3<​-0 +
-for (i in 1:3+
-  Spec.s1s2<​-Spec.s1s2+abs(sum(spe.challenge[1,​i]-spe.challenge[2,​i])+
-  Spec.s1s3<​-Spec.s1s3+abs(sum(spe.challenge[1,​i]-spe.challenge[3,​i])) +
-  Spec.s2s3<​-Spec.s2s3+abs(sum(spe.challenge[2,​i]-spe.challenge[3,​i])) } +
-</​code>​+
  
-Maintenant, utilisez les différences calculées comme numérateur et l'​abondance totale de l'​espèce comme dénominateur pour retrouver l'​équation de la distance de Bray-Curtis. +<file rsplus | Générer des données non linéaires ​
-<code rsplus | Distance de Bray-Curtis+basic_model = gam(y~x0+s(x1), data= gam_data
-(db.s1s2<​-Spec.s1s2/​(Abund.s1+Abund.s2)) #1 comparé à 2 +basic_summary = summary(basic_model) 
-(db.s1s3<​-Spec.s1s3/​(Abund.s1+Abund.s3)) #1 comparé à 3 +print(basic_summary$p.table
-(db.s2s3<​-Spec.s2s3/​(Abund.s2+Abund.s3)#2 comparé à 3  +print(basic_summary$s.table) 
-</code>+plot(basic_model
 +</file>
  
-Vérifiez vos résultats ​en utilisant ​la fonction vegdist ​() +Ici, la sortie de ''​p.table''​ fournit le tableau de résultats ​pour chaque terme paramétrique et le tableau ''​s.table''​ nous donne les résultats du terme non linéaire. Notez que pour le second tableau, ​la courbure du terme non linéaire ''​s(X1)''​ est indiquée par le paramètre edf (degrés de libertés effectifs); plus la valeur de l'edf est élevée, plus la non-linéarité est forte. Une valeur élevée ​(8-10 ou plus) signifie que la courbe est fortement non linéaire, alors qu'une courbe avec un edf égal à 1 est une ligne droiteEn revanchedans la régression linéaire, les degrés de libertés du //modèle// sont équivalents au nombre de paramètres libres non redondants p dans le modèle (et les degrés de libertés //​résiduels//​ sont égaux à n-p). Nous reviendrons plus tard sur le concept d'edf.
-<code rsplus | Vérification des résultats>​ +
-(spe.db.challenge<​-vegdist(spe.challengemethod="​bray"​)) +
-</code>+
  
-Une matrice comme celle-ci est calculée et devrait être correspondre à vos calculs manuels: +   ​print(basic_summary$p.table) 
-^ ^Site 1^Site 2^ +               Estimate Std. Error   t value     ​Pr(>​|t|) 
-^Site 2^0.5^--^ +   ​(Intercept) 8.550030 ​ 0.3655849 23.387258 ​1.717989e-76 
-^Site 3^0.538^0.0526^+   x02         2.418682  ​0.5165515 ​ 4.682364 3.908046e-06 
 +   x03         ​4.486193  ​0.5156501 ​ 8.700072 9.124666e-17 
 +   ​x04 ​        ​6.528518  ​0.5204234 12.544629 1.322632e-30 
 +   > print(basic_summary$s.table) 
 +              edf   ​Ref.df ​       F      p-value 
 +   s(x1) 1.923913 2.406719 42.84268 1.076396e-19
  
-Pour la distance de Gower, procédez de la même façon, mais utiliser l'​équation appropriée:​ 
-Distance de Gower: d[jk] = (1/M) sum(abs(x[ij]-x[ik])/​(max(x[i])-min(x[i]))) 
-<code rsplus | Distance de Gower> 
-# Calculer le nombre de colonnes 
-M<​-ncol(spe.challenge) 
  
-# Calculer les différences d'​abondance ​de chaque espèce entre paires de sites +Dans notre modèle ​de basel'edf du terme non linéaire s(x1est ~ 2, ce qui indique une courbe non linéaireLe graphique du modèle illustre bien la forme de ce terme non linéaire :
-Spe1.s1s2<​-abs(spe.challenge[1,1]-spe.challenge[2,​1]) +
-Spe2.s1s2<​-abs(spe.challenge[1,​2]-spe.challenge[2,​2]) +
-Spe3.s1s2<​-abs(spe.challenge[1,​3]-spe.challenge[2,3]) +
-Spe1.s1s3<​-abs(spe.challenge[1,​1]-spe.challenge[3,​1]) +
-Spe2.s1s3<​-abs(spe.challenge[1,​2]-spe.challenge[3,​2]) +
-Spe3.s1s3<​-abs(spe.challenge[1,​3]-spe.challenge[3,​3]) +
-Spe1.s2s3<​-abs(spe.challenge[2,​1]-spe.challenge[3,​1]) +
-Spe2.s2s3<​-abs(spe.challenge[2,​2]-spe.challenge[3,​2]) +
-Spe3.s2s3<​-abs(spe.challenge[2,​3]-spe.challenge[3,​3])+
  
-# Calculer l'​étendue d'​abondance de chaque espèces parmi les sites +{{::​graphic2.2.jpg?425|}}
-Range.spe1<​-max(spe.challenge[,​1]) - min (spe.challenge[,​1]) +
-Range.spe2<​-max(spe.challenge[,​2]) - min (spe.challenge[,​2]) +
-Range.spe3<​-max(spe.challenge[,​3]) - min (spe.challenge[,​3])+
  
-# Calculer la distance de Gower +Nous pouvons ajouter un second terme x2, mais spécifier une relation linéaire avec Y (//i.e.// les GAMs peuvent inclure à la fois des termes linéaires et non linéaires dans le même modèle). Ce nouveau terme linéaire x2 sera présenté dans le tableau ''​p.table'',​ pour lequel une estimation du coefficient de régression sera indiquéeDans le tableau ''​s.table'',​ nous retrouvons encore une fois le terme non linéaire s(x1et son paramètre de courbure.
-(dg.s1s2<​-(1/M)*((Spe2.s1s2/Range.spe2)+(Spe3.s1s2/Range.spe3))) +
-(dg.s1s3<​-(1/M)*((Spe2.s1s3/Range.spe2)+(Spe3.s1s3/Range.spe3))) +
-(dg.s2s3<​-(1/​M)*((Spe2.s2s3/​Range.spe2)+(Spe3.s2s3/​Range.spe3)))+
  
-# Vérifier vos résultats +<file rsplus | deux termes> 
-(spe.db.challenge<-vegdist(spe.challengemethod="​gower"​)+two_term_model ​<- gam(y~x0+s(x1)+x2data=gam_data
-</code> +two_term_summary ​<- summary(two_term_model) 
-</hidden>+print(two_term_summary$p.table) 
 +print(two_term_summary$s.table) 
 +</file>
  
 +Pour évaluer si la relation entre y et x2 est non linéaire, on peut modéliser x2 avec une fonction non linéaire. Tel que vu auparavant, nous pouvons utiliser une ANOVA pour tester si le terme non linéaire est nécessaire.
  
 +<file rsplus | Deux termes non linéaires >
 +two_smooth_model <- gam(y~x0+s(x1)+s(x2),​ data=gam_data)
 +two_smooth_summary <- summary(two_smooth_model)
 +print(two_smooth_summary$p.table)
 +print(two_smooth_summary$s.table)
 +plot(two_smooth_model,​page=1)
 +</​file>​
  
-=====2.2 Transformations des données de composition des communautés=====+{{::​graphic2.4.jpg?​600|}}
  
-Les données de composition des communautés peuvent également ​être standardisées ou transforméesLa fonction decostand ​() de vegan fournit des options ​de standardisation et de transformation ​de ce type de données+Lorsqu'​il y a plus d'une variable d'​incluse dans le modèle, comme ci-dessus, la réponse ajustée peut-être partitionnée entre les contributions de chaque variableIci, nous pouvons évaluer l'​effet de chaque variable où l'axe des ordonnées représente la contribution ​(effet) de chaque covariable à la réponse ajustée centrée sur 0. Si l'​intervalle ​de confiance chevauche zéro pour certaines valeurs ​de x, cela indique que l'​effet est non significatif. Lorsque la contribution varie selon l'axe x, un changement ​de cette variable cause un changement ​de la variable réponse.
  
-Transformer les abondances en données de présence-absence:​ +<file rsplus | ANOVA
-<code rsplus | decostand, présence-absence+anova(basic_model,​two_term_model,​two_smooth_modeltest="Chisq") 
-spe.pa<​-decostand(spemethod="pa")  +</file>
-</code>+
  
-D'​autres transformations peuvent être utilisées pour corriger l'​influence d'​espèces rares, par exemple, la transformation de Hellinger+    Analysis of Deviance Table 
-<code rsplus | Hellinger et Chi-carré+    Model 1y ~ x0 + s(x1) 
-#La transformation Hellinger +    Model 2: y ~ x0 + s(x1) + x2 
-spe.hel<-decostand(spe,​ method="​hellinger"​) # vous pouvez aussi simplement écrire "​hel"​+    Model 3: y ~ x0 + s(x1) + s(x2) 
 +      Resid. Df Resid. Dev      Df Deviance ​ Pr(>Chi)    ​ 
 +    ​1 ​   394.08 ​    ​5231.6 ​                               
 +    ​2 ​   393.10     ​4051.3 0.97695 ​  ​1180.2 ​2.2e-16 *** 
 +    3    385.73 ​    ​1839.5 7.37288 ​  ​2211.8 < 2.2e-16 *** 
 +    --- 
 +    Signif. codes: ​ 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
-#​Transformation de chi-carré +Le meilleur modèle est le modèle avec deux fonctions non linéaires pour x1 et pour x2.
-spe.chi<​-decostand(spe,​ method="​chi.square"​) +
-</​code>​+
  
-**Défi option 2 Niveau avancé** +----
-Calculez les distances de Hellinger et de Chi-carré sur les données "​spe"​ sans utiliser decostand ().+
  
-**Défi option 2 - Solution** 
-<​hidden>​ 
-La transformation de Hellinger est une transformation qui diminue l'​importance accordée aux espèces rares. 
  
-<code rsplus | Solution>​ +**DÉFI 2**
-# Hellinger  +
-# Calculer l'​abondance des espèces par site  +
-(site.totals=apply(spe,​ 1, sum))+
  
-# Réduire ​les abondances d'espèces en les divisant par les totaux par sites +Créez deux nouveaux modèles avec la variable x3 : un modèle avec x3 comme paramètre linéaire et un autre modèle avec x3 avec un paramètre non linéaire. Utilisez des graphiques, ​les tables des coefficients et la fonction ​''anova()''​ afin de déterminer s'il est nécessaire d'​inclure x3 dans le modèle.
-(scale.spe<​-spe/​site.totals)+
  
-# Calculer la racine carrée des abondances d'​espèces réduites +++++ Réponse au défi 2| 
-(sqrt.scale.spe<​-sqrt(scale.spe))+
  
-# Comparer les résultats +<file rsplus | Réponse au défi 2> 
-sqrt.scale.spe +three_term_model <- gam(y~x0+s(x1)+s(x2)+x3,​ data=gam_data) 
-spe.hel +three_smooth_model <- gam(y~x0+s(x1)+s(x2)+s(x3),​ data=gam_data) 
-sqrt.scale.spe-spe.hel # ou: sqrt.scale.spe/​spe.hel+three_smooth_summary <summary(three_smooth_model)
  
-# Chi-carré +print(three_smooth_summary$p.table) 
-# Premièrement calculer le total des abondances d'​espèces par site +print(three_smooth_summary$s.table) 
-(site.totals<​-apply(spe, 1, sum))+plot(three_smooth_model,page=1) 
 +# edf = 1 -> le terme est donc linéaire.
  
-Ensuite calculer la racine carrée du total des abondances d'espèces +anova(two_smooth_model,​three_term_model,​test="​Chisq"​) 
-(sqrt.spe.totals<-sqrt(apply(spe,​ 2, sum)))+le terme x3 n'est pas significatif 
 +</file>
  
-# Réduire les abondances d'​espèces en les divisant par les totaux par sites et les totaux par espèces +++++ 
-scale.spe2<​-spe +===== 3Interactions =====
-for (i in 1:​nrow(spe)) { +
-  for (j in 1:​ncol(spe)) { +
-   ​(scale.spe2[i,​j]=scale.spe2[i,​j]/​(site.totals[i]*sqrt.spe.totals[j])) ​  }}+
  
-#Ajuster les abondances en les multipliant par la racine carrée du total de la matrice des espèces +Il y a deux façons ​de modéliser une interaction entre deux variables: 
-(adjust.scale.spe2<-scale.spe2*sqrt(sum(rowSums(spe))))+  * si une variable est quantitative et l'​autre est qualitative,​ on utilise l'​argument ''​by''​ -> s(x, by=facteur),​ 
 +  * si les deux variables sont quantitatives,​ on inclut les deux termes sous une même fonction non linéaire -> s(x1, x2). 
 +L'​argument ''​by''​ permet de faire varier un terme non linéaire selon les différents niveaux d'un facteurNous allons examiner ceci en utilisant notre variable qualitative ''​x0''​ et examiner si la non-linérité de ''​s(x2)''​ varie selon les différents niveaux de ''​x0''​. Pour déterminer si les courbes diffèrent significativement entre les niveaux du facteur, nous allons utiliser une ANOVA sur l'​interaction.
  
-#Vérifier les résultats +<file rsplus | Interaction qualitative > 
-adjust.scale.spe2 +categorical_interact <- gam(y~x0+s(x1)+s(x2,​by=x0),​data=gam_data) 
-spe.chi +categorical_interact_summary <- summary(categorical_interact) 
-adjust.scale.spe2-spe.chi # or: adjust.scale.spe2/​spe.chi +print(categorical_interact_summary$s.table) 
-</​code>​ +plot(categorical_interact,​page=1) 
-</hidden>+# ou nous pouvons utiliser la fonction vis.gam où theta représente la rotation du plan x-
 +vis.gam(categorical_interact,​view=c("​x2","​x0"​),​theta=40,​n.grid=500,​border=NA) ​ 
 +anova(two_smooth_model,​ categorical_interact,​test="​Chisq"​) 
 +</file> 
 +{{::​graphic3.1b.png?​350|}}
  
 +Nous pouvons constater à partir du graphique que les formes des termes non linéaires sont comparables entre les quatre niveaux de ''​x0''​. L'​ANOVA le confirme également (déviance = 98,6, p = 0,2347).
  
 +Ensuite, nous allons examiner l'​interaction non linéaire entre deux termes quantitatifs,​ ''​x1''​ et ''​x2''​. Cette fois-ci, l'​argument ''​by''​ est supprimé.
  
-=====2.3 Groupement=====+<file rsplus | Interaction quantitative>​ 
 +smooth_interact <- gam(y~x0+s(x1,​x2),​data=gam_data) 
 +smooth_interact_summary <- summary(smooth_interact) 
 +print(smooth_interact_summary$s.table) 
 +plot(smooth_interact,​page=1,scheme=3) 
 +# plot(smooth_interact,​page=1,scheme=1) donne un graphique comparable à vis.gam() 
 +vis.gam(smooth_interact,​view=c("​x1","​x2"​),​theta=40,n.grid=500,border=NA)  
 +anova(two_smooth_model,​smooth_interact,​test="​Chisq"​) 
 +</​file>​
  
-Les matrices d’association nécessaires pour utiliser les méthodes de groupementLe groupement n’est pas une méthode statistique en tant que telle puisqu’elle ne teste pas d’hypothèse,​ mais permet de déceler des structures dans les données en partitionnant soit les objets, soit les descripteurs. Les objets similaires sont agrégés en sous-groupes ce qui permet de mettre en relief des cassures (contrastes) entre les données. En tant que biologiste, il peut être intéressant de tenter de séparer une série de sites en groupes en fonction de leurs caractéristiques environnementales ou de leur composition en espèces.+{{::​graphic3.2b.png?350|}}
  
-Les résultats d’un groupement sont généralement représentés sous forme de dendrogramme ​(arbre), dans lequel les objets sont agrégés ​en groupesIl a plusieurs familles ​de méthodes ​de groupements,​ mais nous présenterons uniquement un aperçu de trois méthodes : le groupement agglomératif hiérarchique à liens simples, le groupement agglomératif hiérarchique à liens complets et la méthode ​de WardPour plus de détails sur les différentes familles de méthodes de groupementconsulter Legendre ​et Legendre 2012 (chapitre 8)+L'​interaction entre ''​s(x1)''​ et ''​s(x2)''​ est significative et le graphique ​en deux dimensions illustre très bien cette interaction non linéaireLa relation entre et x1 change en fonction ​de la valeur ​de x2. Vous pouvez changez ​la valeur ​de l'​argument ''​theta''​ pour tourner l'axe du graphiqueSi vous prévoyez exécuter un grand nombre ​de graphiquessupprimez l'​argument ''​n.grid = 500'',​ car ceci fait appel à des calculs intensifs ​et ralentit R.
  
-Dans les méthodes hiérarchiques,​ les éléments des petits ensembles se regroupent en groupes plus vastes de rang supérieur, et ainsi de suite (par exemple : espèces, genres, familles, ordre)Avant de faire le groupement, il faut créer une matrice d’association entre les objets. Une matrice de distances est le choix par défaut des fonctions de groupement dans R. La matrice d’association est premièrement classée en ordre croissant de distances. Ensuite, les groupes sont formés de manière hiérarchique selon les critères spécifiques à chaque méthode.+===== 4Changer la fonction ​de base =====
  
-Prenons un exemple tout simple ​d'une matrice ​de distances euclidiennes entre 5 objets dont on a ordonné ​les distances en ordre croissant+Sans entrer dans le détail, sachez qu'il est possible de modifier le modèle de base que nous avons vu avec : 
 +  * des fonctions plus complexes en modifiant la fonction de base (par exemple, cyclique),​ 
 +  * d'​autres distributions : tout ce que vous pouvez faire avec un GLM (tel que spécifier l'​argument ''​family''​) est possible avec les GAMs, 
 +  * des modèles à effets mixtes en utilisant la fonction //gamm// ou la fonction //gamm4// de la librairie //​gamm4//​. 
 +Nous allons ​d'abord jeter un coup d’œil au changement ​de la fonction de base puis une introduction rapide aux autres distributions et les GAMMs (modèles additifs généralisés à effets mixtes) suivra. Commençons par regarder un cas où modifier la fonction de base peut être utile, soit avec des données cycliques.
  
-{{ :groupement.exemple.png |}}+Imaginez que vous avez une série temporelle de données climatiques,​ divisées en mesures mensuelles, et que vous voulez déterminer s'il y a une tendance de température annuelle. Nous allons utiliser la série temporelle de température de Nottingham pour cette section ​:
  
-Pour le groupement agglomératif à liens simples, les deux objets les plus proches se regroupent en premier. Ensuite, un deuxième groupe est formé à partir des deux objets les plus proches suivant ​(il se peut que ce soit deux objets différents,​ ou bien un objet et le groupe formé précédemment), et ainsi de suite. Cette méthode forme généralement de longues chaînes de groupes ​(dans l'​exemple ci-hautles objets ​à 5 se regroupent successivement). À l’inversepour le groupement agglomératif à liens completun objet se regroupe à un autre objet/​groupe seulement lorsqu’il est aussi lié à l’élément le plus éloigné de ce groupe. Ainsi, quand deux groupes fusionnenttous les éléments des deux groupes sont liés à la distance considérée ​(ci-hautle groupe 3-4 ne se lie au groupe 1-2 qu'à la distance à laquelle tous les autres éléments sont déjà liés). C’est pour cette raison que le groupement à liens complets forme généralement plusieurs petits groupes séparés et qu’elle peut être plus appropriée pour relever des contrastes ou des discontinuités dans les données.+<file rsplus | Fonction de base cyclique>​ 
 +data(nottem) 
 +n_years <- length(nottem)/​12 
 +nottem_month <rep(1:12times=n_years) 
 +nottem_year <- rep(1920:​(1920+n_years-1),each=12) 
 +nottem_plot <- qplot(nottem_month,nottem 
 +                    colour=factor(nottem_year) 
 +                    geom="​line"​+ theme_bw() 
 +print(nottem_plot) 
 +</​file>​
  
-Comparons ces deux méthodes en utilisant les données d’abondances ​de poissons de la rivière Doubs. +En utilisant les données ​//nottem//, nous avons créé trois nouveaux vecteurs ; ''​n_years''​ correspond au nombre ​d'​années ​de données ​(20 ans), ''​nottem_month''​ est un codage qualitatif pour les 12 mois de l'​année,​ pour chaque année échantillonnée (série ​de 1 à 12répétée 20 fois) et ''​nottem_year''​ est une variable où l'​année correspondant à chaque mois est fournie.
-Les données ​d’abondances ont été au préalable transformées par la méthode Hellinger. Puisque ​les fonctions ​de groupement requièrent une matrice ​de distancesla première étape sera de générer ​une matrice de distances Hellinger.+
  
-<code rsplus | vegdist, distance de Hellinger>​ +Données mensuelles ​des années 1920 à 1940:
-spe.dhel<​-vegdist(spe.hel,​method="​euclidean"​) #crée une matrice de distances Hellinger à partir ​des données d’abondance transformées+
  
-#Pour voir la différence entre les deux types d’objets +{{::​graphic4.1.jpg?​500|}} ​
-head(spe.hel)# données d’abondances transformées Hellingerhead(spe.dhel)# matrice de distances de Hellinger entre les sites +
-</​code>​+
  
-La plupart des méthodes ​de groupement sont disponible dans la fonction hclust() de la librairie stats +Pour modéliser cela, nous devons utiliser ce qu'on appelle un "​spline cubique cyclique",​ ou ''​cc'',​ pour modéliser les effets du mois et de l'​année. 
 +<file rsplus | gam cyclique > 
 +year_gam <- gam(nottem~s(nottem_year)+s(nottem_month,​ bs="​cc"​)) 
 +summary(year_gam)$s.table 
 +plot(year_gam,​page=1,​ scale=0) 
 +</​file>​
  
-<code rsplus | hclust, comparaison du groupement à liens simples et à liens complets>​ +{{::​graphic4.2.jpg?700|}}
-#Faire le groupement à liens simples +
-spe.dhel.single<​-hclust(spe.dhel,​ method="​single"​) +
-plot(spe.dhel.single)+
  
-#Faire le groupement à liens complet +Il y a une hausse d'​environ 1 - 1,5ºC au cours de la série, mais au cours d'une année, il y a une variation d'​environ 20ºC. ​ 
-spe.dhel.complete<​-hclust(spe.dhel, method="​complete"​) +Les données réelles varient autour de ces valeurs prédites et ceci représente donc la variance inexpliquéeIci, nous pouvons voir l'un des avantages très intéressants de l'​utilisation des GAMsNous pouvons soit tracer la surface réponse ​(valeurs préditesou les termes ​(contribution de chaque covariable) tel qu'​indiqué ci-hautVous pouvez imaginer ce dernier en tant qu'une illustration de la variation des coefficients de régression et comment leur contribution (ou taille de leur effetvarie au fil du temps. Dans le premier graphique, nous voyons que les contributions positives de la température sont survenues après 1930.
-plot(spe.dhel.complete) +
-</​code>​+
  
-{{ :clust_single.png |}}{{ :​clust_complete.png |}}+Sur des échelles de temps plus longues, en utilisant par exemple des données paléolimnologiques,​ d'​autres ([[http://www.aslo.info/​lo/​toc/​vol_54/​issue_6_part_2/​2529.pdf|Simpson & Anderson 2009;. Fig 3c]]) ont utilisé des GAMs pour tracer la contribution (effet) de la température sur la composition d'​algues dans les lacs afin d'​illustrer comment les contributions significatives ont seulement eu lieu au cours de deux périodes extrêmement froides (c'​est-à-dire,​ la contribution est importante lorsque les intervalles de confiance ne recoupent pas la valeur de zéro à environ 300 et 100 ans AVJC). Cela a permis aux auteurs de non seulement déterminer combien de variance est expliquée par la température au cours des derniers siècles, mais aussi de repérer dans le temps cet effet significatif. Si cela vous intéresse, le code pour tracer soit la surface de réponse (''​type = "​response"''​) ou les termes (''​type = "​terms"''​) est disponible ci-dessous. Lorsque les termes sont sélectionnés,​ vous obtiendrez la même figure que celle ci-dessus.
  
-Est-ce que les deux dendrogrammes sont très différents?​+++++ Graphique de contribution vs réponse ajustée|  
 +<file rsplus | contribution plot> 
 +pred<-predict(year_gam,​ type = "​terms",​ se = TRUE) 
 +I<​-order(nottem_year) 
 +plusCI<​-I(pred$fit[,​1] + 1.96*pred$se[,​1]) 
 +minusCI<​-I(pred$fit[,​1] - 1.96*pred$se[,​1]) 
 +xx <- c(nottem_year[I],​rev(nottem_year[I])) 
 +yy <- c(plusCI[I],​rev(minusCI[I])) 
 +plot(xx,​yy,​type="​n",​cex.axis=1.2) 
 +polygon(xx,​yy,​col="​light grey",​border="​light grey"​) 
 +lines(nottem_year[I],​ pred$fit[,​1][I],​lty=1) 
 +abline(h=0) 
 +</​file>​ 
 +++++
  
-On remarque que pour le groupement à liens simple, plusieurs objets s’enchaînent (par exemple les sites 19, 29, 30, 20, 26, etc.) alors que des groupes plus distincts peuvent être observés dans le groupement à liens complets+===== 5Intro rapide aux GAMMs =====
  
-La méthode de Ward diffère légèrement ​des deux méthodes précédentes. Le critère utilisé est la méthode des moindres carrés (comme dans les modèles linéaires). Ainsi, des objets/​groupes fusionnent de façon à minimise la variance intragroupes. Pour débuter, chaque objet est considéré comme un groupe. À chaque étape, la paire de groupes à fusionner est celle qui résulte à la plus petite augmentation de la somme des carrés des écarts intra-groupes.+==== La non-indépendance ​des données ====
  
-La méthode ​de Ward est également disponible sous la fonction hclust(). Par contre, le dendrogramme produit par défaut montre ​les distances ​au carréAfin de comparer ce dendrogramme à celui du groupement à liens simples et à liens complets, il faut calculer ​la racine carrée des distances+Lorsque les observations ne sont pas indépendantes,​ les GAMs peuvent être utilisés soit pour incorporer:​ 
 +  * une structure ​de corrélation pour modéliser les résidus autocorrélés (autorégressif ​(AR), moyenne mobile (MA), ou une combinaison des deux (ARMA))  
 +  * des effets aléatoires qui modélisent l'​indépendance entre les observations d'un même site. 
 +En plus de changer la fonction de basenous pouvons aussi complexifier ​le modèle en intégrant une structure d'​auto-corrélation (ou même des effets mixtes) en utilisant ​les fonctions //gamm// ou //​gamm4//​. 
 +Pour commencer, nous allons jeter un coup d’œil ​au premier cas ; un modèle avec autocorrélation temporelle dans les résidusRé-examinons le modèle ​de la température de Nottingham; nous allons vérifier si les résidus sont corrélés en faisant appel à la fonction (partielle) d'​autocorrélation.
  
-<code rsplus | hclust, méthode de Ward+<file rsplus| ​Erreurs corrélées
-#Faire le groupement de Ward  +par(mfrow=c(1,​2)) 
-spe.dhel.ward<​-hclust(spe.dhelmethod="ward.D2") +acf(resid(year_gam),​ lag.max = 36main = "ACF") 
-plot(spe.dhel.ward)+pacf(resid(year_gam),​ lag.max = 36, main = "​pACF"​) 
 +</​file>​
  
-#Refaire le dendrogramme en utilisant la racine carrée des distances +{{::​graphic5.1.jpg?700|}}
-spe.dhel.ward$height<​-sqrt(spe.dhel.ward$height) +
-plot(spe.dhel.ward) +
-plot(spe.dhel.ward,​ hang=-1) # hang=-1 permet d’afficher les objets sur la même ligne +
-</​code>​+
  
-{{ :clust_ward.png |}} +Les graphiques des fonctions d'​autocorrélation suggèrent qu'un modèle AR de faible ordre est nécessaire (avec un ou deux intervalles de temps décalés), donc nous pouvons évaluer deux modèles; ajouter un AR(1) ou un AR(2) au modèle de la température de Nottingham et évaluer le meilleur avec une ANOVA.
-{{ :​clust_wardfinal.png |}}+
  
-Les groupements générés par la méthode de Ward ont tendance à être plus sphériques et à contenir des quantités plus similaires d’objets.+<file rsplus| modèles AR> 
 +year_gam <- gamm(nottem~s(nottem_year)+s(nottem_month,​ bs="​cc"​)) 
 +year_gam_AR1 <- gamm(nottem~s(nottem_year)+s(nottem_month,​ bs="​cc"​),​ 
 +                     ​correlation = corARMA(form = ~ 1|nottem_year,​ p = 1)) 
 +year_gam_AR2 <- gamm(nottem~s(nottem_year)+s(nottem_month,​ bs="​cc"​),​ 
 +                     ​correlation = corARMA(form = ~ 1|nottem_year,​ p = 2)) 
 +anova(year_gam$lme,​year_gam_AR1$lme,​year_gam_AR2$lme) 
 +</​file>​
  
-Quelle méthode choisir?+                   Model df      AIC      BIC    logLik ​  ​Test ​  ​L.Ratio p-value 
 +  year_gam$lme ​        ​1 ​ 5 1109.908 1127.311 -549.9538 ​                         
 +  year_gam_AR1$lme ​    ​2 ​ 6 1101.218 1122.102 -544.6092 1 vs 2 10.689206 ​ 0.0011 
 +  year_gam_AR2$lme ​    ​3 ​ 7 1101.598 1125.962 -543.7988 2 vs 3  1.620821 ​ 0.2030 
 +Le modèle avec la structure AR(1) prévoit une augmentation significative comparativement au premier modèle (LRT = 10,69, p = 0,0011), mais il y a très peu d'​intérêt à considérer le modèle AR(2) (LRT = 1,62, b = 0,203). 
 +==== Modélisation avec effets mixtes ====
  
-Le choix de la bonne mesure d’association et de la bonne méthode ​de groupement dépend de l’objectif. Qu’est-ce qu’il est plus intéressant de démontrer : des gradients? des contrastes? Il est également important de tenir en compte les propriétés de la méthode utilisée dans l’interprétation des résultatsSi plus d’une méthode semble adéquate pour répondre à une question biologique, comparer ​les dendrogrammes serait une bonne option. Encore une foisle groupement n’est pas une analyse statistique,​ mais il est possible de tester les résultats ​et d’identifier des partitions ayant un sens biologique. Il est également possible de déterminer le nombre de groupes optimal et de performer des tests statistiques sur les résultats. Les méthodes de groupement peuvent aussi être combinées à une ordination pour distinguer des groupes de sites. Ces avenues ne seront pas explorées dans cet atelier. Pour aller plus loinconsulter Borcard et al. 2011.+Comme nous l'​avons vu dans la section précédente,​ ''​bs''​ spécifie ​la fonction ​de base sous-jacentePour les facteurs aléatoires (origine et pente linéaire)nous utilisons ''​bs = "​re"'' ​et pour les pentes aléatoires non linéairesnous utilisons ''​bs = "​fs"''​.
  
 +Trois types d'​effets aléatoires différents sont possibles lors de l'​utilisation des GAMMs (où //fac// représente une variable qualitative utilisée pou l'​effet aléatoire et //x0// est un effet quantitatif fixe) :
 +  * //​interceptes aléatoires//​ ajustent la hauteur des termes du modèle avec une valeur constante de pente : s(fac, bs="​re"​)
 +  * //pentes aléatoires//​ ajustent la pente d'une variable explicative numérique: s(fac, x0, bs="​re"​)
 +  * //surfaces lisses aléatoires//​ ajustent la tendance d'une prédiction numérique de façon non linéaire: s(x0, fac, bs="​fs",​ m=1) où l'​argument m = 1 met une plus grande pénalité au lissage qui s'​éloigne de 0, ce qui entraîne un retrait vers la moyenne.
  
 +Nous examinerons d'​abord un GAMM avec un interception aléatoire. Tel que vu précédemment,​ nous allons utiliser ''​gamSim()''​ pour générer un ensemble de données, cette fois-ci avec une composante d'​effet aléatoire. Ensuite, nous construirons un modèle avec un intercepte aléatoire en utilisant //fac// comme facteur aléatoire.
  
-======3. Ordination sans contrainte======+<file rsplus| Intercepte aléatoire > 
 +# Générez des données 
 +gam_data2 <- gamSim(eg=6) 
 +head(gam_data2)
  
-Les analyses d'​ordination sans contrainte permettent d'​organiser des échantillons,​ des sites ou des espèces le long de gradients continus ​(ex. écologiques ou environnementaux). Les ordinations sans contrainte se différencient des analyses canoniques ​(voir plus loin dans cet atelier ​"analyses canoniques"​) ​par le fait que ces techniques ne tentent pas de définir une relation entre des variables dépendantes et indépendantes.+# Faites rouler un modèle avec intercepte aléatoire 
 +gamm_intercept <- gam(y ~ s(x0+ s(fac, bs="re"), data=gam_data2) 
 +summary(gamm_intercept) 
 +</​file>​
  
-L'​ordination sans contrainte peut être utilisée pour: +Notez le terme aléatoire ​dans le tableauVous pouvez le visualiser:
-- Évaluer les relations //au sein// d'un ensemble de variables (et non pas entre séries de variables). +
-- Trouver les éléments clés de variation au sein d'​échantillons,​ de sites, d'​espèces,​ etc. +
-- Réduire les dimensions d'un jeu de données multivariées sans perte importante d'​informations. +
-- Créer de nouvelles variables pour une utilisation ​dans des analyses ultérieures (comme la régression)Les composantes principales des axes d'​ordination sont en effet des combinaisons linéaires des variables d'​origine. +
-[[http://​www.umass.edu/​landeco/​teaching/​multivariate/​schedule/​ordination1.pdf|Source]]+
  
-=====3.1 Analyses en Composantes Principales ​(Principal Component AnalysisPCA)=====+<file rsplus| graphique intercepte aléatoire>​ 
 +plot(gamm_interceptselect=2 
 +# select=2 parce que le terme aléatoire se trouve sur la 2e ligne du tableau sommaire. 
 +</​file>​
  
-L'Analyse en Composantes Principales (ou PCA) fur originellement décrite ​par Pearson (1901) bien qu'elle soit le plus souvent attribué à Hotelling (1933) qui l'a proposé indépendammentCette méthode ainsi que nombre de ces implications sont pour l'analyse de données sont présentées dans l'article fondateur ​de Rao (1964). La PCA est utilisée pour générer, à partir d'un large jeu de données, un nombre restreint de variables clefs qui permettent de représenter au maximum la variation présente ​dans le jeu de donnéesEn d'​autres termesla PCA est utilisée pour générer des combinaisons de variables à partir d'un ensemble plus grand de variables tout en conservant la majorité de variation ​de l'ensemble des données. La PCA est une technique d'​analyse puissante ​pour l'​analyse dans descripteurs quantitatifs (tels que les abondances d'​espèces),​ mais ne peut pas être appliquée aux données binaires (telles que l'​absence/présence des espèces).+Une fonction de traçage vraiment intéressante que nous allons maintenant utiliser est le ''​plot_smooth'' de la librairie //​itsadug//​. Contrairement au graphique ​par défaut ​''​plot.gam'',​ cette fonction présente ​l'effet additionné du GAMM avec l'option ​de ne pas inclure les courbes aléatoires ​dans le graphiqueIcinous allons premièrement tracer l'effet combiné ​de x0 (sans les niveaux ​de l'effet aléatoire) et ensuite ​une courbe ​pour les quatre niveaux de //fac//:
  
-À partir d'un jeu de données contenant ​des variables à distribution normalele premier axe de PCA (ou axe de composante principalecorrespond à la droite qui traverse la plus grande dimension de l’ellipsoïde décrivant la distribution multi-normale des donnéesLes axes suivants traversent de façon similaire cet ellipsoïde selon l'​ordre décroissant de ces dimensionsAinsiil est possible d'obtenir un maximum de p axes principaux à partir d'un jeu de données contenant p variables.+<file rsplus | Graphique ​des 4 niveaux du facteur fac> 
 +par(mfrow=c(1,2), cex=1.1) 
 +plot_smooth(gamm_intercept,​ view="​x0",​ rm.ranef=TRUE,​ main="​intercept + s(x1)", rug=FALSE) 
 +plot_smooth(gamm_intercept,​ view="​x0",​ cond=list(fac="​1"​),​  
 +            main="... + s(fac)"​col='orange', ylim=c(8,​21),​ rug=FALSE) 
 +plot_smooth(gamm_intercept,​ view="​x0",​ cond=list(fac="​2"​),​ add=TRUE, col='​red'​) 
 +plot_smooth(gamm_intercept,​ view="​x0",​ cond=list(fac="​3"​),​ add=TRUE, col='​purple'​) 
 +plot_smooth(gamm_intercept,​ view="​x0",​ cond=list(fac="​4"​),​ add=TRUE, col='​turquoise'​) 
 +</​file>​
  
-Pour cela, la PCA effectue un rotation du système d'axes originels défini par les variables de façon à ce que les axes d'​ordination successifs soient orthogonaux entre eux et correspondent aux dimensions successives du maximum de variance observée dans le nuage de points (voir ci-dessus)Les nouvelles variables produites par la PCA sont non-corrélées entre elles (les axes d'​ordination étant orthogonaux) et peuvent alors être utilisés dans d'​autres types d'​analyse telles que des régressions multiples (Gotelli et Ellison, 2004). Les composantes principales situent la position des objets dans le nouveau système de coordonnées calculé par la PCA. La PCA s'​effectue sur une matrice d'​association entre variables, et a pour caractéristique de préserver les distances euclidiennes et de détecter des relations linéaires, uniquement. En conséquence,​ les abondances brutes des espèces doivent être soumises à une pré-transformation (comme une transformation d'​Hellinger) avant de réaliser une PCA.+{{::​graphic5.3.jpg?650|}}
  
-//Pour faire une PCA, vous avez besoin ​:// +Ensuite, nous allons générer et tracer un modèle avec une pente aléatoire ​:
-- Un ensemble de variables (sans distinction entre les variables indépendantes ou dépendantes,​ c-est-à-dire un ensemble d'​espèces OU un ensemble de variables environnementales). +
-- De sites (objets, échantillons) dans lesquels sont mesurés les mêmes variables. +
-- En général, il est préférable d'​avoir un plus grand nombre de sites que de variables dans le jeu de données (plus de lignes que de colonnes).+
  
-La PCA est particulièrement utile pour des matrices de données contenant plus de deux variables, mais il est plus facile de décrire son fonctionnement avec un exemple bidimensionnel. Dans l'​exemple suivant ​(d'​après Clarke et Warwick 2001), la matrice contient les données d'​abondance de deux espèces dans neuf sites:+<file rsplus | pente aléatoire > 
 +gamm_slope <- gam(y ~ s(x0+ s(x0fac, bs="​re"​),​ data=gam_data2) 
 +summary(gamm_slope)
  
-^Site^Espèce ​1^Espèce 2^ +plot_smooth(gamm_slope,​ view="​x0",​ rm.ranef=TRUE,​ main="​intercept + s(x0)",​ rug=FALSE) 
-^A^6^2^ +plot_smooth(gamm_slope,​ view="​x0",​ cond=list(fac="​1"​), ​ 
-^B^0^0^ +            ​main="​... + s(fac)",​ col='​orange',​ylim=c(7,22), rug=FALSE) 
-^C^5^8^ +plot_smooth(gamm_slope,​ view="​x0",​ cond=list(fac="​2"​),​ add=TRUE, col='​red'​) 
-^D^7^6^ +plot_smooth(gamm_slope,​ view="​x0",​ cond=list(fac="​3"​),​ add=TRUE, col='​purple'​) 
-^E^11^6^ +plot_smooth(gamm_slope,​ view="​x0",​ cond=list(fac="​4"​),​ add=TRUE, col='​turquoise'​) 
-^F^10^10^ +</​file> ​      
-^G^15^8^ +
-^H^18^14^ +
-^I^14^14^+
  
-La représentation des sites en deux dimensions devrait ressembler à ceci: +{{::graphic5.4.jpg?650|}}    
-{{:pcaex_1.png?500|}}+
  
-Ce nuage de points est une ordination. Il présente ​la distribution des espèces entre les sites, mais vous pouvez imaginer qu'il est plus difficile de visualiser ​un tel graphique en présence de plus de deux espècesDans ce cas, l'​objectif est de réduire le nombre de variables en composantes principales. Pour réduire les données bidimensionnelles précédentes à une dimension, une PCA peut être effectuée :+Nous allons maintenant inclure à la fois un intercepte et une pente aléatoires.
  
-{{:​pcaex_2.png?​500|}}+<file rsplusIntercepte et pente aléatoires > 
 +gamm_int_slope <- gam(y ~ s(x0) + s(fac, bs="​re"​)  
 +                      + s(fac, x0, bs="​re"​),​ data=gam_data2) 
 +summary(gamm_int_slope)
  
-Dans ce casla première composante principale est orientée dans le sens de la plus grande variation dans les pointsces points étant perpendiculaires à la ligne.+plot_smooth(gamm_int_slopeview="​x0"​rm.ranef=TRUE, main="​intercept + s(x0)",​ rug=FALSE) 
 +plot_smooth(gamm_int_slope,​ view="​x0",​ cond=list(fac="​1"​),​  
 +            main="​... + s(fac) + s(fac, x0)", col='​orange',​ ylim=c(7,​22),​ rug=FALSE) 
 +plot_smooth(gamm_int_slope,​ view="​x0",​ cond=list(fac="​2"​),​ add=TRUE, col='​red',​ xpd=TRUE) 
 +plot_smooth(gamm_int_slope,​ view="​x0",​ cond=list(fac="​3"​),​ add=TRUE, col='​purple',​ xpd=TRUE) 
 +plot_smooth(gamm_int_slope,​ view="​x0",​ cond=list(fac="​4"​),​ add=TRUE, col='​turquoise',​ xpd=TRUE) 
 +</​file>​
  
-Une seconde composante principale est alors ajoutée perpendiculairement à la première:+{{::​graphic5.5.jpg?​650|}} ​
  
-{{:pcaex_3.png?​500|}}+Notez que les pentes aléatoires sont statique ​:
  
-Dans le diagramme final, les deux axes de PCA sont pivotés et les axes sont maintenant les composantes principales ​(et non plus les espèces):+<file rsplus| graphique des pentes aléatoires>​ 
 +plot(gamm_int_slope,​ select=3 
 +# select=3 parce que la pente aléatoire se trouve sur la 3e ligne du tableau sommaire. 
 +</​file>​
  
-{{:pcaex_4.png?500|}}+Enfin, nous allons examiner un modèle avec une surface lisse aléatoire.
  
-Pour les PCA avec plus de deux variablesles composantes principales sont ajoutées de la façon suivante ​(Clarke et Warwick 2001):+<file rsplus| Surface lisse aléatoire > 
 +gamm_smooth <- gam(y ~ s(x0) + s(x0fac, bs="​fs",​ m=1), data=gam_data2) 
 +summary(gamm_smooth) 
 +</​file>​
  
-PC1 = axe qui maximise la variance des points qui sont projetées perpendiculairement à l'​axe. +Icisi les pentes aléatoires variaient selon x0, nous auront vue des courbe variable pour chaque niveau ​:
-PC2 = axe perpendiculaire à PC1mais dont la direction est à nouveau celle maximisant la variance lorsque ​les points y sont projetés perpendiculairement. +
-PC3 et ainsi de suiteperpendiculaire aux deux premiers axes dont la direction est à nouveau celle maximisant la variance lorsque les points y sont projetés perpendiculairement.+
  
-Lorsqu'​il y a plus de deux dimensionsla PCA produit un nouvel espace dans lequel tous les axes de PCA sont orthogonaux (ce qui signifie ​que la corrélation entre chaque combinaison de deux axes est nulle), et où les axes de PCA sont ordonnés selon la proportion de variance des données d'​origine qu'ils expliquent.+<file rsplus| graphiques des fonctions aléatoires>​ 
 +plot(gamm_smoothselect=1)  
 +# select=1 parce que le terme se trouve sur la 1e ligne du tableau sommaire. 
 +</​file>​
  
-Les données "​spe"​ comprennent 27 espèces de poissons. Pour simplifier cette diversité à un petit nombre de variables ou pour identifier différents groupes de sites associés à des espèces particulièresune PCA peut être effectuée.+Finalementtous ces modèles mixes peuvent ​être compare en utilisant la fonction anova() pour trouver le meilleur modèle.
  
-Exécuter une PCA sur les données d'​abondance d'​espèces soumises à la transformation d'​Hellinger : +===== 6Autres distributions =====
-<code rsplus | PCA avec rda() de (vegan)>​ +
-#Exécuter la PCA avec la fonction rda()- cette fonction calcule à la fois des PCA et des RDA +
-spe.h.pca<​-rda(spe.hel)+
  
-#​Extraire ​les résultats +Pour vous donner un bref aperçu de l'​utilisation des GAMs lorsque la variable réponse ne suit pas une distribution normale ou que les données sont des abondances ou proportions ​(par exemple, distribution Gamma, binomiale, Poisson, binomiale négative), l'​exemple qui suit utilise un ensemble de données où une répartition binomiale sera nécessaire,​ y compris une modélisation d'une relation non linéaireLa variable réponse représente le nombre de succès (l'​événement a eu lieuen fonction des défaillances au cours d'une expérience.
-summary(spe.h.pca +
-</​code>​+
  
-Résultats: ​+<file rsplus| Loading the data> 
 +gam_data3 <- read.csv("​other_dist.csv"​) 
 +summary(gam_data3) 
 +str(gam_data3) 
 +</​file>​
  
-{{:pca_outputfr_1.png?800|}} +  '​data.frame'​: 514 obsof  4 variables: 
-{{:pca_outputfr_2.png?800|}} +   $ prop num  1 1 1 1 0 1 1 1 1 1 ..
-{{:pca_outputfr_3.png?800|}}+   $ totalint  4 20 20 18 18 18 20 20 20 20 ... 
 +   $ x1   : int  550 650 750 850 950 650 750 850 950 550 ... 
 +   $ fac  : Factor w/ 4 levels "​f1","​f2","​f3",​..:​ 1 1 1 1 1 1 1 1 1 1 ...
  
-**Interprétation des résultats de PCA** +''​prop'' ​est la variable réponseégal à la proportion ​de //succès / (succès + échecs)//. Notez qu'il existe ​de nombreux cas où la proportion est égal à 1 ou 0 qui indique que les résultats ​ont toujours été des succès ou des échecs, respectivement,​ à ce moment mesuré durant l'​expérience.\\ 
-Cette sortie R contient les valeurs propres ou «eigenvalue» de la PCA. La valeur propre ​est la valeur de la variation ramenée à la longueur d'un vecteuret correspond ​à la quantité ​de variation expliquée par chaque axe d'ordination ​de la PCA. Comme vous pouvez le voir, la fonction summary fournit de nombreuses informations. Parmi les résultats, la proportion de variance ​des données expliquée par les variables sans contraintes est une information importanteDans cet exemple, la variance totale des sites expliquée par les espèces ​est de 0,5 (50%). Le résumé vous indique également quelle proportion de la variance totale expliquée est répartie entre chaque composantes principales de la PCA: le premier axe de PCA explique 51.33% ​de la variation tandis que le second axe explique 12.78%Vous pouvez également extraire certaines parties des résultats: +''​x1'' ​est le temps écoulé depuis le début ​de l'​expérience ​(variable explicative).\\ 
-<code rsplus | Résultats ​de PCA> +''​total''​ représente ​le nombre ​de //succès + échecs// observé au moment x1i de l'​expérience.\\ 
-summary(spe.h.pca, display=NULL# seulement les valeurs propres +''​fac''​ est un facteur qui code pour l'​essai 1 à 4 de l'​expérience ​(nous n'​utiliserons pas cette variable dans cette section).
-eigen(cov(spe.hel)) # vous pouvez aussi trouver les valeurs propres par cette ligne de code +
-</​code>​+
  
-Les scores (c'​est-à-dire les coordonnées) des sites ou des espèces peuvent également être extraits d'une PCA. Ces scores permettent, ​par exemple, d'​utiliser une composante principale comme une variable dans une autre analyse, ou de faire des graphiques supplémentairesPar exemple, vous pouvez obtenir ​par PCA une variable unique issue du jeu de données "​spe"​ puis l'​utiliser pour la corréler ​par régression ​à une autre variable, ou déterminer un gradient spatialPour extraire les scores d'une PCA, utiliser ​la fonction ​ scores ():  +Commençons ​par la visualisation ​des donnéesNous sommes intéressés ​par le nombre ​de succès ​par rapport aux échecs ​à mesure que x1 augmenteÉtant donné qu'il y a des mesures répétées pour la valeur de x1 (essais ​à 4avec nombreuses observations par essai), nous pouvons d'abord présenter la proportion ​de succès en moyenne par boîte de temps (x1):
-<code rsplus | scores()>​ +
-spe.scores<​-scores(spe.h.pca,​ display="​species",​ choices=c(1,2)) # scores des espèces selon les premier et deuxième axes  +
-site.scores<​-scores(spe.h.pcadisplay="​sites",​ choices=c(1,​2)) # scores des sites selon les premier et deuxième axes  +
-#Remarque: si vous ne spécifiez pas le nombre de composantes principales à l'aide de choices = c (1,2 +
-#(ou choices = c (12)), les scores selon toutes les composantes principales seront extraits. +
-</​code>​+
  
-La PCA des données d'​abondances de poissons produit autant de composantes principales qu'il y d'​espèces (i.e. de colonnes dans le jeu de données), soit 27 composantes principales. Le nombre de variables à traiter n'est donc pas directement réduit par la PCA. Pour réduire le nombre de variables, il est alors nécessaire de déterminer quelles composantes principales sont significatives et doivent être conservées,​ par exemple à l'aide du critère de Kaiser-Guttman. Ce critère compare la variance expliquée par chaque composante principale à la moyenne de la variance expliquée par l'​ensemble des composantes principales. Un histogramme illustrant la significativité des différentes composantes principale peut ensuite être tracé à l'aide du code ci-dessous :  +<file rsplus| ​Visualisation ​des données
-<code rsplus | Les axes> +emptyPlot(range(gam_data3$x1), c(0,1), h=.5, 
-# Identification ​des axes significatifs de la PCA à l'aide du critère de Kaiser-Guttman +          main="Probability of successes", ​ylab="Probability",xlab="x1")
-ev<​-spe.h.pca$CA$eig +
-ev[ev>mean(ev)] +
-n<​-length(ev) +
-bsm<​-data.frame(j=seq(1:n), p=0+
-bsm$p[1]=1/n +
-for (i in 2:n+
-  bsm$p[i]=bsm$p[i-1]+(1/​(n=1-i))} +
-bsm$p=100*bsm$p/​n +
-bsm +
-barplot(ev, ​main="valeurs propres", ​col="grey", ​las=2) +
-abline(h=mean(ev),​ col="red"​) ​ +
-legend("​topright",​ "​moyenne des valeurs propres",​ lwd=1, col=2, bty="​n"​) +
-</​code>​+
  
-{{:​pca_sigaxes_sp.png?500|}}+avg <- aggregate(prop ~ x1, data=gam_data3,​ mean, na.rm=TRUE) 
 +lines(avg$x1,​ avg$prop, col="​orange",​lwd=2) 
 +</​file>​
  
-Cet histogramme montre que la proportion de la variance expliquée par chaque composante chute en-dessous de la proportion moyenne expliquée par l'​ensemble des composantes après le sixième axe d"​ordination PC6. En consultant le résumé de nouveau, on peut constater que la proportion cumulée de la variance expliquée par les cinq premières composantes principales est de 85%.+{{::​graphic1_otherdist.jpg?450|}}
  
-La PCA n'est pas seulement appropriée pour les données ​de composition d'espècesmais peut également être exécutée sur des variables environnementales standardisées:​ +Notez comment la probabilité ​de succès augmente avec x1. D'après vousest-ce que cette tendance est linéaire ​ou non linéaire? Nous allons tester cela en utilisant un GAM logistique ​(nous utilisons une distribution ''​binomiale''​ puisque la variable réponse représente des proportions).
-<code rsplus | PCA pour les variables environnementales>​ +
-#Exécuter la PCA +
-env.pca<-rda(env.z) # ou rda(env, scale=TRUE)+
  
-#Extraction des résultats +<file rsplus| GAM logistique>​ 
-summary(env.pca+prop_model <- gam(prop~ s(x1), data=gam_data3,​ weights=total,​ family="​binomial"​
-summary(env.pca, scaling=2)  +prop_summary <- summary(prop_model) 
-</​code>​+print(prop_summary$p.table
 +print(prop_summary$s.table)
  
-«Scaling» réfère à quelle portion de la PCA est redimensionnée aux valeurs propres. Scaling = 2 signifie que les scores des espèces sont mises à l'​échelle des valeurs propres, alors que scaling = 1 signifie que les scores des sites sont mises à l'​échelle des valeurs propres. Scaling = 3 signifie qu'à la fois les scores des espèces et des sites sont mis symétriquement à l'​échelle de la racine carrée des valeurs propres. En scaling = 1,les distances euclidiennes entre les sites (lignes de la matrice de donnéessont conservées tandis qu'en scaling = 2 les corrélations entre espèces (les colonnes de la matrice de données) sont conservées. Cela implique que lorsque vous regardez un biplot de PCA en Scaling = 2, l'​angle entre les descripteurs représente leur corrélation.+plot(prop_model) 
 +</​file>​
  
-<code rsplus | Axes significatifs>​ +               ​Estimate ​ StdError  z value   Pr(>|z|
-ev<-env.pca$CA$eig +  (Intercept 1.173978  ​0.02709613 ​ 43.32641  0
-ev[ev>mean(ev)] +
-n<​-length(ev) +
-bsm<​-data.frame(j=seq(1:​n),​ p=0+
-bsm$p[1]=1/​n +
-for (i in 2:n) { +
-  bsm$p[i]=bsm$p[i-1]+(1/​(n=1-i))} +
-bsm$p=100*bsm$p/​n +
-bsm +
-barplot(ev, main="​valeurs propres",​ col="​grey",​ las=2) +
-abline(h=mean(ev),​ col="​red"​)  +
-legend("​topright",​ "​moyenne des valeurs propres",​ lwd=1, col=2, bty="​n"​) +
-</​code>​ +
-Comparer cet histogramme de valeurs propres avec celui que vous avez créé pour la PCA sur les abondances d'​espèces.+
  
-Bien que beaucoup d'​informations puissent être extraites d'une PCA par la fonction summary de la PCA, l'​interpétation et la communication des résultats est souvent facilitée en traçant un biplotSur un biplot de PCA, l'axe des x correspond à la première composante principale et l'axe des y à la deuxième composante principale.+          edf       Ref.df    Chi.sq     ​p-value 
 +  s(x1)   ​4.591542 ​ 5.615235 ​ 798.9407 ​  ​2.027751e-164
  
-La fonction plot () permet de tracer des biplot sur lequels les sites figurent en chiffres noirs et les espèces sont représentées en rouge. Le biplot de la  //PCA des abondances d'​espèces//​ peut être appelé comme suit: 
-<code rsplus | Graphique PCA des abondances d'​espèces > 
-plot(spe.h.pca) 
-</​code> ​ 
  
-{{:​spe_PCA1.png?800|}}+Qu'​est-ce que l'​ordonnée représente dans ce modèle? 
 +  * Rappel : le modèle utilise le nombre de succès vs échecs pour calculer le //logit//, qui est le logarithme du rapport entre les succès et échecs:
  
-La construction de biplots de PCA s'​articule en trois étapes:  +<mlogit log({N_{success}/​N_{failures}}) </m>
-<code rsplus | Biplot de PCA> +
-plot(spe.h.pca,​ type=”n”) # Produit une figure vierge +
-points(spe.h.pca, dis=”sp”,​ col=”blue”) # ajoute les points correspondant aux espèces +
-#utilizer text() plutôt que points() si vous préférez que les codes des espèces s'​affichent (nom des colonnes) +
-points(spe.h.pca,​ dis=”sites”,​ col=”red”# ajoute les points correspondant aux sites  +
-</code>+
  
-{{:spe_PCA2.png?800|}}+- Si succès = échecs, le rapport est de 1 et le logit est 0 (log (1) = 0).\\ 
 +- Si les succès ont un nombre plus grand que les échecs, le ratio est supérieur à 1 et le logit a une valeur positive (par exemple, log(2) = 0,69).\\ 
 +- Si les succès ont un nombre plus petit que les échecs, le ratio est inférieur à 1 et le logit a une valeur négative (par exemple, log(0,5) = -0.69).\\
  
-Pour créer de plus beaux biplotsessayez ce code: +Doncl'ordonnée ​est le //logit//, et indique s'il y a en moyenne plus de succès que d'​échecsIcil'​estimé est positif ce qui signifiequ'en moyenneil n'y a plus de succès que d'​échecs.
-<code rsplus | Graphique PCA des abondances d'espèces v.2> +
-#Scaling 1 +
-windows() +
-plot(spe.h.pca) +
-windows() +
-biplot(spe.h.pca) +
-windows() +
-# scaling 1 = distance biplot :  +
-# distances entre les objets ​est une approximation de leur distance euclidienne +
-# les angles entre les descripteurs ne réflètent PAS leur corrélation +
-plot(spe.h.pca,​ scaling=1, type="​none",​  +
-     ​xlab<​-c("​PC1 (%)", round((spe.h.pca$CA$eig[1]/sum(spe.h.pca$CA$eig))*100,​2)),​  +
-     ​ylab<​-c("​PC2 (%)", round((spe.h.pca$CA$eig[2]/sum(spe.h.pca$CA$eig))*100,2))) +
-points(scores(spe.h.pcadisplay="​sites"​choices=c(1,2), scaling=1),​ +
-       ​pch=21,​ col="​black",​ bg="​steelblue",​ cex=1.2) +
-text(scores(spe.h.pca,​ display="​species",​ choices=c(1),​ scaling=1),​ +
-     ​scores(spe.h.pca,​ display="​species",​ choices=c(2),​ scaling=1),​ +
-     ​labels=rownames(scores(spe.h.pca,​ display="​species",​ scaling=1)),​ +
-     ​col="​red",​ cex=0.8) ​  +
-</​code>​+
  
-Le code ci-dessus a produit trois biplots mais le dernier est le plus attrayant : +Qu'est-ce que le terme de lissage indique? 
-{{:spe_PCA3.png?800|}}+  * Ceci représente la façon dont les chances de succès vs échecs changent sur l'​échelle de x1 (l'​échelle du temps dans cet exemple). Donc, puisque l'edf > 1, la proportion de succès augmente plus rapidement au fil du temps (si par exemple, la réponse représente ​le nombre d'​individus de l'​espèce A vs l'​espèce B et que nous augmentons la concentration des nutriments au fil du temps, ces résultats indiqueront que l'​espèce A est de plus en plus observée alors que les concentrations de nutriments approchent de l'​optimum de cette espèce au cours de l'​expérience).
  
-Sur ce graphique, les sites scores sont indiqués par des points bleus et les noms d'​espèces sont en rouge. Il est également possible de représenter les sites par leurs noms.+=== Visualiser la tendance au fil du temps ===
  
-Comment interpréter ce type de graphique ?  +Enfin, nous allons voir les différentes façons ​de représenter ces relations graphiquement.
-Ce biplot permet d'​observer qu'il existe certains groupes de sites homogènes du point de vue de la composition de leur communautés de poissons. On y voit également que l'​espèce «ABL» n'a pas la même prévalence dans la majorité des sites que les autres espèces plus proches du centre du graphique.+
  
-Les biplots ne doivent pas seulement être interprétés en termes de proximitémais également d'​angles. Deux variables séparées d'un angle de 90 degrés ne sont pas corrélées. Deux variables très rapprochés sont fortement corrélées. Deux variables aux directions opposées sont corrélées négativement. ​+<file rsplus| Tracer la tendance>​ 
 +par(mfrow=c(1,2)) 
 +plot(prop_model,​ select=1, scale=0, shade=TRUE) 
 +abline(h=0)
  
-Maintenant regardons le biplot de la //PCA environnement//:​ +plot_smooth(prop_modelview="​x1"​,main=""​) 
-<code rsplus | Graphique PCA des variables environnementales>​ +(diff <- find_difference(out$fv$fitout$fv$CI, xVals=out$fv$x1)) 
-#Scaling 2 +addInterval(0lowVals=diff$start, highVals = diff$endcol='​red'​lwd=2
-windows(+abline(v=c(diff$startdiff$end), lty=3, col='​red'​
-plot(env.pca) +text(mean(c(diff$startdiff$end)), 2.1, "sign. more \n success", col='red'font=3
-windows() +</file>
-# scaling 2 = graphique de corrélations : +
-# les distances entre les objets ne sont PAS des approximations de leur distance euclidienne +
-# les angles entres les descripteurs reflètent leur corrélation +
-plot(env.pcascaling=2type="none" +
-     xlab<-c("PC1 (%)"round((env.pca$CA$eig[1]/​sum(env.pca$CA$eig))*100,2)), +
-     ylab<-c("PC2 (%)"round((env.pca$CA$eig[2]/​sum(env.pca$CA$eig))*100,2)), +
-     ​xlim<​-c(-1,​1),​ ylim=c(-1,1)+
-points(scores(env.pca,​ display="​sites",​ choices=c(1,2), scaling=2),  +
-       ​pch=21, col="​black",​ bg="​darkgreen",​ cex=1.2)  +
-text(scores(env.pca, display="​species",​ choices=c(1)scaling=2)+
-     ​scores(env.pca,​ display="​species",​ choices=c(2), scaling=2), +
-     ​labels<​-rownames(scores(env.pcadisplay="species", ​scaling=2)),​ +
-     col="red"cex=0.8+
-</code+
  
-{{:env_PCA1.png?800|}}+{{::​graphic2_otherdist.jpg?800|}}
  
-Rappelez-vous qu'un biplot de PCA est en fait un nuage de points dans lequel ​les axes sont des combinaisons linéaires des variables ​d'origineIl existe donc beaucoup ​de façons différentes ​de tracer un biplotPar exemplevous pouvez utiliser la fonction ggplot () et les compétences acquises ​de l'atelier 4 pour tracer votre graphique ​d'​ordination dans ggplot.+Quels renseignements ces graphiques nous apportent-ils vis à vis les succès et échecs ? 
 +  * **graphique de gauche**: contribution (ou effet partiel si nous avions plus qu'une variable explicative) au fil du temps. La valeur logit augmente, donc les succès augmentent et les échecs diminuent. 
 +  * **graphique de droite**: valeurs ajustées, ordonnée incluse (somme ​des effets si nous avions plus d'une variable explicative dans le modèle)Nous voyons ici que la valeur logit est estimée près de zéro au début ​de l'​expérience ; cela signifie qu'il y a des quantités égales de succès et d'​échecsPeu à peu, les succès augmentent et à environ x1=400 il y a beaucoup plus de succès que d'​échecs (l'effet est significativement différent de zéro). Nous avons également montré comment nous pouvons utiliser le graphique ​pour déterminer à quelle valeur de x1 cela se produit.
  
 +Enfin, pour nous aider à interpréter les résultats, nous pouvons re-transformer l'​effet sur une échelle de proportions avec la fonction ''​plot_smooth''​ de la librairie //​itsadug//:​
  
-**Utilisation des axes de PCA comme variable explicative composite**+<file rsplus| Graphique transformé>​ 
 +par(mfrow=c(1,​1)) 
 +plot_smooth(prop_model,​ view="​x1",​ main="",​ 
 +            transform=plogis,​ ylim=c(0,​1)) 
 +abline(h=.5,​ v=diff$start,​ col='​red',​ lty=2) 
 +</​file>​
  
-Dans certains cas, l'​utilisateur cherche à réduire un grand nombre de variables environnementales en un plus faible nombre de variables composites. Lorsque les axes de PCA représentent des gradients écologiques (i.e. lorsue les variables environnementales sont corrélées de façon cohérente avec les axes de PCA), l'​utlisateur peut utiliser les scores des sites le long axes de PCA dans de nouvelles analyses (au lieu d'​utiliser les variables environnementales brutes). En d'​autres termes, étant donné que les scores des sites le long des axes de PCA représentent des combinaisons linéaires des descripteurs,​ ils peuvent être utilisés comme proxy des conditions écologiques dans de nouvelles analyses.+{{::​graphic3_otherdist.jpg?450|}}
  
-Dans l'exemple ci-dessus, ​le premier axe de PCA peut être identifié comme un gradient écologique allant ​des sites oligotrophes riches en oxygène aux sites eutrophes pauvres en oxygène: de gauche à droitele premier groupe de sites montrent les plus hautes altitudes (alt) et pente (slope), et les plus faibles débit (deb) et distance ​à la source (das). Le second groupe ​de sites possède les plus hautes valeurs ​de concentration en oxygène (oxy) et les plus faibles concentrations en nitrates (nit). Un troisième groupe ​de sites montrent des valeurs intermédiaire pour l'​ensemble de ces variables.+Comme nous l'avons déjà vu avec le graphique précédent ​des valeurs logitsnous voyons qu'à approximativement x1=400 ​la proportion ​de succès augmente ​de façon significative au-dessus ​de 0,5.
  
-Dans ce cas, si l'​objectif est d'​identifier si une espèce particulière est associée au gradient oligotrophe-eutrophe,​ il est possible de corréler l'​abondance de cette espèce aux scores des sites le long du premier axe de PCA. Par exemple, si l'​utilisateur veut identifier si l'​espèce TRU est associée à des eaux oligotrophes ou eutrophes, il lui est possible d’utiliser le modèle linéaire suivant: ​ 
  
-<code rsplus | Graphique PCA des variables environnementales>​ +===== 7Les GAMs en coulisse ======
-Sites_scores_Env_Axis1<​- scores(env.pca, display="​sites",​ choices=c(1), scaling=2) +
-spe$ANG +
-plot( Sites_scores_Env_Axis1,​ spe$TRU) +
-summary(lm(spe$TRU~Sites_scores_Env_Axis1)) +
-abline(lm(spe$TRU~Sites_scores_Env_Axis1)) +
-</​code> ​+
  
-Ce modèle simple montre que l'​abondance de l'​espèce TRU est significativement liée aux scores des sites le long du premier axe de PCA (t = -5.30, p = 1.35e-05, adj-R2 = 49.22%), c’est-à-dire qu'​elle dépend ​d'​un ​gradient oligotrophe-eutrophe. dans ce cas l'espèce TRU préfère donc les eaux oligotrophes.+Nous allons maintenant prendre quelques minutes pour regarder comment fonctionnent les GAMsCommençons en considérant ​d'abord un modèle qui contient une fonction lisse d'une covariable, x<​sub>​i</​sub>​ :
  
-**Défi 3** +<m>{y_i} = f(x_{i}+ {ε_i}</m>
-Exécuter une PCA sur l'​abondance des espèces de mites. Quels sont les axes significatifs ? Quels sont groupes de sites pouvez-vous identifier? Quelles espèces sont liées à chaque groupe de sites?  +
-<code rsplus | Données mites> +
-mite.spe<​-data(mite# données disponibles dans vegan +
-</code+
  
-**Défi ​Solution** +Pour estimer la fonction //f//, nous avons besoin de représenter l'​équation ci-dessus de manière à ce qu'​elle devienne un modèle linéaire. Cela peut être fait en choisissant une base, b<sub>i</sub>(x), définissant l'​espace des fonctions dont //f// est un élément:
-<hidden> +
-Votre code ressemble certainement à celui-ci: +
-<code rsplus | Mite PCA>+
  
-# Transformation de Hellinger  +<m>f(x) = sum{i=1}{q}{b_{i}(x)β_{i}}</m>
-mite.spe.hel<-decostand(mite.spe,​ method="​hellinger"​) +
-mite.spe.h.pca<​-rda(mite.spe.hel) +
-  +
-# Quels sont les axes significatifs?​  +
-ev<​-mite.spe.h.pca$CA$eig  +
-ev[ev>mean(ev)+
-n<​-length(ev) +
-bsm<​-data.frame(j=seq(1:n), p=0) +
-bsm$p[1]=1/n +
-for (i in 2:n) { +
-  bsm$p[i]=bsm$p[i-1]+(1/​(n=1-i))} +
-bsm$p=100*bsm$p/+
-bsm +
-barplot(ev, main="​Valeurs propres",​ col="​grey",​ las=2) +
-abline(h=mean(ev),​ col="​red"​)  +
-legend("​topright",​ "​Moyenne des valeurs propres",​ lwd=1, col=2, bty="​n"​) ​+
  
-# Résultats ​    +**Un exemple simple : une base polynomiale**
-summary(mite.spe.h.pca,​ display=NULL)  +
-windows()+
  
-# Représentation graphique de la PCA +Supposons que //f// est considéré comme un polynôme d'​ordre 4de sorte que l'​espace des polynômes d'​ordre 4 et moins contiennent //f//. Une base de cet espace serait alors :
-plot(mite.spe.h.pca,​ scaling=1, type="​none",​  +
-     ​xlab=c("​PC1 (%)", round((mite.spe.h.pca$CA$eig[1]/sum(mite.spe.h.pca$CA$eig))*100,2)), +
-     ​ylab=c("​PC2 (%)", round((mite.spe.h.pca$CA$eig[2]/sum(mite.spe.h.pca$CA$eig))*100,​2)))  +
-points(scores(mite.spe.h.pca,​ display="​sites",​ choices=c(1,​2),​ scaling=1),​ +
-       ​pch=21,​ col="​black",​ bg="​steelblue",​ cex=1.2) +
-text(scores(mite.spe.h.pca,​ display="​species",​ choices=c(1),​ scaling=1),​ +
-     ​scores(mite.spe.h.pca,​ display="​species",​ choices=c(2),​ scaling=1),​ +
-     ​labels=rownames(scores(mite.spe.h.pca,​ display="​species",​ scaling=1)),​ +
-     ​col="​red",​ cex=0.8) ​  +
-</​code> ​+
  
-Votre graphique ressemblera à ceci. +<m>b_{1}(x)=1, ​  b_{2}(x)=x, ​  ​b_{3}(x)=x^{2}, ​  ​b_{4}(x)=x^{3}, ​  ​b_{5}(x)=x^{4}</m>
-{{:​mite_pca.png?​800|}}+
  
-Bien que les sites soient tous de composition semblable ​(aucun groupe distinct de sites n'​apparait sur le biplot), certaines espèces semblent souvent être présentes ensemble, par exemple Spec 01, Spec 10, Spec 14 et Spec 15. +et f(xdevient:
-</​hidden>​+
  
-=====3.2. Analyse des Correpondances ​(Correspondence Analysis, CA)=====+<​m>​f(x) ​β_{1} + x_{i}β_{2} +  x^{2}_{i}β_{3} + x^{3}_{i}β_{4}(x+ x^{4}_{i}β_{5}</​m>​
  
-L'une des hypothèses clefs de la PCA postule que les espèces sont liées les unes aux autres de façon linéaire, ​et qu'​elles répondent de façon linéaire aux gradients écologiques. Ce n'est cependant pas nécessairement ​le cas dans les données écologiques (e.g. beaucoup d'​espèces montrent en effet un distribution unimodale le long des gradients environnementaux). Utiliser une PCA avec des données contenant des espèces à distribution unimodale, ou un grand nombre de zéros (absence des espèces), peut conduire à un phénomène statistique appelé "​horseshoe effect"​ (ou effet fer à cheval) se produisant le long de gradients écologiques. Dans de tels cas, l'​Analyse des Correspondances (CA) permet de mieux représenter les données (voir Legendre et Legendre pour plus d'​informations). Comme la CA préserve les distances de Chi2 entre objets (tandis que la PCA préserve les distances euclidiennes),​ cette technique est, en effet, plus appropriée pour ordonner les jeux de données contenant des espèces à distribution unimodale, et a, pendant longtemps, était l'une des techniques les plus employées pour analyser les données d'​absence-présence ou d'​abondances d'​espèces. Lors d'une CA, les données brutes sont d'​abord transformées en une matrice Q des contributions cellule-par-cellule à la statistique Chi2 de Pearson, puis la matrice résultante est soumise à une décomposition en valeurs singulières afin de calculer les valeurs propres et vecteurs propres de l'​ordination.+et le modèle complet devient:
  
-Le résultat d'une CA représente donc une ordination dans laquelle les distances de Chi2 entre objets sont préservées (au lieu de la distance euclidienne dans une PCA), le distance de Chi2 n'​étant pas influencée par la présence de double-zéros. Ainsi, la CA constitue une méthode d’ordination puissante pour l'​analyse des abondances brutes d'​espèces (i.e. sans pré-transformation). Contrairement à la PCA, la CA peut être appliquée sur des données quantitatives ou binaires ​(telles que des abondance ou absence-présence d'​espèces). Comme dans une PCA, le critère de Kaiser-Guttman peut être utilisé pour identifier les axes significatifs d'une CA, et les scores des objets le long des axes d'​ordination peuvent être extraits pour être utlisés dans des régressions multiples par exemple.+<​m>​y_{i} = β_{1} + x_{i}β_{2} +  x^{2}_{i}β_{3} + x^{3}_{i}β_{4}(x+ x^{4}_{i}β_{5} + ε_{i}</​m>​
  
-Exécuter une CA sur les données d'​abondance d'​espèces:​+Chaque fonction de base est multipliée par un paramètre à valeur réelle, β<​sub>​i</​sub>,​ et est ensuite additionnée pour donner la courbe finale **f(x)**.
  
-<code rsplus | CA par cca() de vegan> +{{::​polynomial_basis_example.png?600|}}
-#Effectuer une CA à l'aide de la fonction cca() (NBcca() est utilisée à la fois pour les CA et CCA) +
-spe.ca <- cca(spe) +
-  +
-# Identifier les axes significatifs +
-ev<​-spe.ca$CA$eig +
-ev[ev>​mean(ev)] +
-n=length(ev) +
-barplot(ev, main="​Eigenvalues",​ col="​grey",​ las=2) +
-abline(h=mean(ev),​ col="​red"​) +
-legend("​topright",​ "​Average eigenvalue",​ lwd=1, col=2, bty="​n"​) +
-</​code> ​+
  
-{{ :​ca_guttmankaiser.png |}}+En faisant varier le coefficient β<​sub>​i</​sub>,​ on peut faire varier la forme de f(x) pour produire une fonction polynomiale d'​ordre 4 ou moins.
  
-D’après cet histogramme,​ à partir du sixième axe d’ordination CA6, la proportion de variance expliquée diminue sous la proportion moyenne expliquée par l'​ensemble des axes. La sortie R de la CA ci-dessous montre également que les cinq premiers axes d'​ordination explique ​une proportion cumulée de variance expliquée ​de 84.63%.+**Un autre exemple : une base de spline cubique**
  
-<code rsplus | Extraire les résultats ​d'une CA> +Un spline cubique est une courbe construite à partir de sections ​d'un polynôme cubique reliées entre elles de sorte qu'​elles sont continues en valeurChaque section du spline a des coefficients différents.
-summary(spe.h.pca)  +
-summary(spe.h.pca, diplay=NULL) +
-</​code>​+
  
-{{ :ca_summary.png |}}+{{::​cubic_spline_fr.png?550|}}
  
-Les résultats d'une CA sont présentés sous R de la même façon que ceux d'​une ​PCA. On y observe que le premier axe CA1 explique 51.50% ​de la variation ​de l'​abondance ​des espèces tandis que le second axe CA2 explique 12.37% ​de la variation.+Voici une représentation ​d'​une ​fonction lisse utilisant une base de régression spline cubique ​de rang 5 avec des nœuds situés à incréments ​de 0,2:
  
-<code rsplus | Construction des biplots>​ +{{::​graphic6.1.jpg?300|}}
-par(mfrow=c(1,​2)) +
-#### scaling 1 +
-plot(spe.ca, scaling=1, type="​none",​ main='​CA - biplot scaling 1', xlab=c("​CA1 (%)", round((spe.ca$CA$eig[1]/​sum(spe.ca$CA$eig))*100,​2)),​ +
-ylab=c("​CA2 (%)", round((spe.ca$CA$eig[2]/​sum(spe.ca$CA$eig))*100,​2)))+
  
-points(scores(spe.ca, display="​sites",​ choices=c(1,​2),​ scaling=1), pch=21, col="​black",​ bg="​steelblue"​cex=1.2)+Dans cet exemple, les nœuds sont espacés uniformément à travers la gamme des valeurs observées de x. 
 +Le choix du degré de finesse du modèle est pré-déterminé par le nombre de noeudsqui était arbitraireY a-t-il une meilleure façon de sélectionner les emplacements des nœuds?
  
-text(scores(spe.ca,​ display="​species",​ choices=c(1),​ scaling=1),​ +**Contrôler le degré de lissage avec des splines de régression pénalisés**
-     ​scores(spe.ca,​ display="​species",​ choices=c(2),​ scaling=1),​ +
-     ​labels=rownames(scores(spe.ca,​ display="​species",​ scaling=1)),​col="​red",​ cex=0.8)+
  
-#### scaling 2 +Au lieu de contrôler le lissage ​(non linéarité) en modifiant le nombre de nœudsnous gardons celle-ci fixée à une taille un peu plus grande que raisonnablement nécessaire et on contrôle le lissage du modèle en ajoutant une pénalité sur le niveau de courbure.
-plot(spe.cascaling=1, type="​none",​ main='​CA ​biplot scaling 2', xlab=c("​CA1 (%)", round((spe.ca$CA$eig[1]/​sum(spe.ca$CA$eig))*100,​2)),​ +
-     ​ylab=c("​CA2 (%)", round((spe.ca$CA$eig[2]/​sum(spe.ca$CA$eig))*100,​2)),​ ylim=c(-2,​3))+
  
-points(scores(spe.cadisplay="​sites",​ choices=c(1,2), scaling=2), pch=21, col="​black",​ bg="​steelblue",​ cex=1.2) +Doncplutôt que d'​ajuster le modèle en minimisant ​(comme avec la méthode des moindres carrés:
-text(scores(spe.ca,​ display="​species",​ choices=c(1),​ scaling=2),​ +
-     ​scores(spe.ca,​ display="​species",​ choices=c(2),​ scaling=2),​ +
-     ​labels=rownames(scores(spe.ca,​ display="​species",​ scaling=2)),​col="​red",​ cex=0.8) +
-</​code>​+
  
-{{ :​ca_biplot.png ​|}}+<​m> ​||y - XB||^{2</m>
  
-Ces biplots montrent qu'un groupe de sites (à gauche) possède des communautés similaires de poissons caractérisées par de nombreuses espèces dont GAR, TAN, PER, ROT, PSO et CAR; dans le coin supérieur droit, un second groupe de sites se caractérisent par les espèces LOC, VAI et TRU; le dernier groupe de sites dans le coin inférieur droit montrent des communautés abondantes ​en BLA, CHA et OMB.+il peut être modélisé ​en minimisant :
  
-**Défi 4** +<m> ||y - XB||^{2} + {lambda}int{0}{1}{[f^{''​}(x)]^{2}dx} </m>
-Exécuter une CA sur les données d'abondance des //espèces d'acariens// ​(données mite). Quels sont les axes importants? Quels groupes de sites pouvez-vous identifier? Quelles espèces sont liées à chaque groupe de sites? ​+
  
-**Défi 4 - Solution** +Quant <mlambda </m> tend vers <m> infty </m>, le modèle devient linéaire. Pour la sélection du meilleur paramètre de lissage, <​m>​lambda</​m>,​ on utilise une approche de validation croisée. Si <​m>​lambda</​m>​ est trop élevé, les données seront trop lissées et si elle est trop faible, les données ne seront pas assez lissées. Idéalement,​ il serait bon de choisir une valeur <​m>​lambda</​m>​ de sorte que le //f prédit// est aussi proche que possible du //f observé//. Un critère approprié pourrait être de choisir <​m>​lambda</​m>​ pour minimiser ​:
-<hidden> +
-Votre code devrait s'​apparenter à celui-ci:+
  
-<code rsplus | CA sur les données mite> +<mM = 1/n sum{i=1}{n}{(hat{f_{i}} - f_{i})^{2}} </m>
-# CA  +
-mite.spe.ca<​-cca(mite.spe)+
  
-# Quels sont les axes importants?  +Étant donné que //f// est inconnuM est estimé en utilisant une technique de validation croisée généralisée qui laisse de côtéà chaque tourune donnée et estime la capacité moyenne des modèlesconstruits sur les données restantes, de prédire la donnée qui a été mise de côté ​(pour plus de détailsconsultez Wood 2006).
-ev<​-mite.spe.ca$CA$eig +
-ev[ev>​mean(ev)] +
-n=length(ev) +
-barplot(evmain="​Eigenvalues"​col="​grey"​las=2) +
-abline(h=mean(ev)col="​red"​) +
-legend("​topright"​"​Average eigenvalue",​ lwd=1, col=2, bty="​n"​)+
  
-# Résultats ​  +**Illustration du principe de validation croisée**
-summary(mite.spe.ca,​ display=NULL)+
  
-# Biplot +{{::​illustration_of_smooth_sel.png?600|}}
-windows() +
-plot(mite.spe.ca, scaling=1, type="​none",​ +
-     ​xlab=c("​PC1 (%)", round((mite.spe.ca$CA$eig[1]/​sum(mite.spe.ca$CA$eig))*100,​2)),​ +
-     ​ylab=c("​PC2 (%)", round((mite.spe.ca$CA$eig[2]/​sum(mite.spe.ca$CA$eig))*100,​2))) +
-points(scores(mite.spe.ca,​ display="​sites",​ choices=c(1,​2),​ scaling=1),​ +
-       ​pch=21,​ col="​black",​ bg="​steelblue",​ cex=1.2) +
-text(scores(mite.spe.ca,​ display="​species",​ choices=c(1),​ scaling=1),​ +
-     ​scores(mite.spe.ca,​ display="​species",​ choices=c(2),​ scaling=1),​ +
-     ​labels=rownames(scores(mite.spe.ca,​ display="​species",​ scaling=1)),​ +
-     ​col="​red",​ cex=0.8) +
-</​code>​+
  
-Et votre biplot devrait ressembler ​à celui-ci:+Dans le premier panneau, la courbe correspond ​à un ajustement faible par rapport aux données et ne fait pas mieux avec le point manquant. Dans le troisième panneau, la courbe ajuste le bruit aussi bien que le signal et la variabilité supplémentaire induite l'​amène à prédire la donnée manquante plutôt mal. Dans le deuxième panneau, cependant, nous voyons que l'​ajustement de la courbe du signal sous-jacent est très bien, le lissage passe à travers le bruit et la donnée manquante est raisonnablement bien prédite.
  
-{{ :ca_mite_biplot.png |}} +{{::gcv.png?600|}}
  
-</​hidden>​+**Note supplémentaire sur les degrés de liberté effectifs (edf)**
  
-=====3.3 Analyse en coordonnées principales (Principal Coordinates Analysis, PCoA)=====+Combien de degrés de liberté a un GAM ?
  
-La PCA, comme la CA, impose une préservation des distances entre objets: ​la distance euclidienne dans le cas de la PCA, et la distance ​de Chi2 dans la CA. Si l'objectif est d'​ordonner les objets sur la base d'une autre mesure de distance plus appropriée au problème, la PCoA constitue une technique de choix. Dans une PCA, les données sont pivotées de façon à ce que la première composante principale ​(correspondant à une combinaison linéaire des descripteursexplique la plus forte proportion de variation possible; la contribution de chaque descripteur (espèces ou variables environnementales) à chaque composante principale peut alors être évaluée d'après son score. +Au lieu de fournir ​la sortie de la validation croisée en termes ​de <​m>​lambda</​m>​ (un paramètre qui détermine ​la complexité ​de l'ajustement), ​la fonction ​''gam()'' ​utilise un terme appelé ​les degrés ​de liberté effectifs ​(edf), de manière cohérente ​à quantifier ​la complexité ​du modèle ​(pour justifier notre intuition que les degrés ​de liberté offrent une manière cohérente de quantifier la complexité du modèle).  
-La PCoA est une seconde méthode d'ordination sans contrainte dans laquelle ​les points sont ajoutés les uns après les autres à l'​espace d'​ordination en utilisant la distance euclidienne //ou n'​importe quelle mesure ​de distance ​(dissimilaritémétrique vous choisissez//​. Un premier point est ainsi placé dans l'​espèce d'​ordinationpuis un second point placé ​à la valeur de distance ​du premier, puis un troisième et ainsi de suite en ajoutant autant d'​axes ​(de dimensionsque nécessaire+Parce que le nombre ​de paramètres libres ​des splines ​de lissage (tel que les GAMs) est souvent difficile à définir, ​les edf sont liés à <​m>​lambda</​m>,​ où l'​effet ​de la pénalité est de réduire les degrés ​de libertésDonc, si le nombre de noeuds est arbitrairement réglé à k = 10, k-1 définit la limite supérieure ​des degrés de libertés associés à un terme de lissage. Ce nombre diminue alors que la pénalité lambda augmente jusqu'​à ce que le meilleur modèle soit trouvé par validation croisée.
-Il est parfois difficile ​de choisir entre effectuer une PCA ou une PCoA. La PCA permet toutefois de réduire ​des données multivariables en un faible nombre ​de dimensions tandis ​que la PCoA est utile pour visualiser ​les distances entre sites (ou objets). La PCoA est aussi particulièrement adaptées pour les jeux de données présentant plus de colonnes que de lignesPar exemple, si des centaines d'​espèces ont été observées dans un petit nombree ​de quadrats, une approche basée sur une PCoA utilisant ​la distance de Bray-Curtis (voir ci-dessous) peut être plus adaptée.+
  
-PCoA avec DoubsSpe (transformé Hellinger):  +===== Références =====
-<code rsplus | PCoA> +
-# En utilisant la fonction cmdscale() +
-?cmdscale +
-cmdscale(dist(spe.hel),​ k=(nrow(spe)-1),​ eig=TRUE)+
  
-# En utilisant la fonction pcoa() +Il existe beaucoup d'​information sur les GAMs...
-?pcoa +
-spe.h.pcoa<​-pcoa(dist(spe.hel))+
  
-# Extraction des résultats +Simon Woodl'auteur ​de la librairie ​//mgcv//, a un site très utile avec des conférences ​et des notes introductives sur la façon ​d'utiliser ​les GAMs : 
-spe.h.pcoa  +  * http://people.bath.ac.uk/​sw283/​mgcv/​ 
- +Il a aussi écrit ​un livre//​Generalized Additive Models: An Introduction with R//que nous avons utilisé comme référence ​pour cet atelier.
-# Représentation graphique +
-biplot.pcoa(spe .h.pcoaspe.hel, dir.axis2=-1) +
-</​code>​ +
- +
-Les résultats de cette PCoA sont:  +
- +
-{{:​pcoa_outputfr_1.png?​800|}} +
- +
-{{:​pcoa_outputfr_2.png?​800|}} +
- +
-Et le graphique:  +
- +
-{{:​pcoa_spe.png?​500|}} +
- +
-Vous pouvez aussi exécuter cette PCoA avec une autre mesure de distance (ex. Bray-Curtis):​  +
- +
-<code rsplus | PCoA avec Bray-Curtis>​ +
-spe.bray.pcoa<​-pcoa(spe.db) # il s'agit de la matrice de distances de Bray-Curtis qu'on a générée plus tôt  +
-spe.bray.pcoa +
-biplot.pcoa(spe.bray.pcoa,​ spe.hel, dir.axis2=-1) +
-# Le choix d'une mesure de distance est très important car ça influence les résultats!  +
-</code> +
- +
-**Défi 5** +
-Exécuter une PCoA sur les données d'​abondance des //espèces d'​acariens// transformées Hellinger (données mite). Quels sont les axes importants? Quels groupes de sites pouvez-vous identifier? Quelles espèces sont liées à chaque groupe de sites? Comment les résultats de cette PCoA se comparent-ils ​avec ceux de la PCA? +
- +
-**Défi 5 - Solution** +
-<​hidden>​ +
-<code rsplus | PCoA avec les données d'​acariens>​ +
-mite.spe.h.pcoa<​-pcoa(dist(mite.spe.hel)) +
-mite.spe.h.pcoa +
-windows() +
-biplot.pcoa(mite.spe.h.pcoa,​ mite.spe.hel,​ dir.axis2=-1) +
-</​code>​ +
- +
-Représentation graphique:  +
- +
-{{:​pcoa_mite.png?​500|}} +
- +
-Les espèces 16 et 31 sont plus éloignées ​des autres espèces en termes de distance, ​et donc leur distribution entre les sites est très différente de celle des autres espèces ​d'acariens. Les sites dont les étiquettes se chevauchent sont de bons exemples de sites à forte similarité en termes de communautés d'​acariens. +
-</hidden>​ +
- +
-=====3.4Positionnement multidimensionnel non-métrique (Nonmetric Multidimensional Scaling, NMDS)===== +
- +
-Les méthodes d'​ordination non contrainte présentées ci-dessus permettent d'​organiser les objets (ex. les sites) caractérisés par des descripteurs (ex. les espèces) dans un espace comprenant l'​ensemble des dimensions décrites par l’ellipsoïde représentant le nuage des points de données. En d'​autres termes, la PCA, la CA et la PCoA calculent un grand nombre d'axes d'​ordination (nombre proportionnel au nombre de descripteurs) représentant la variation des descripteurs entre sites et préservant les distances entre objets (distance euclidienne dans une PCA, distance de Chi2 dans une CA et distance définie par l'​utilisateur dans une PCoA). L'​utilisateur peut ensuite sélectionner les axes d'​intérêt (généralement les deux premiers axes d'​ordination) pour représenter les objets dans un biplot. Le biplot produit représente ainsi correctement les distances entre objets (ex. la similarité des sites), mais ne permet pas de représenter l'​ensemble des dimensions de la variation dans l'​espace d'​ordinations (étant donnée que l'Axe 3, l'Axe 4,..., l'Axe n'​apparaissent pas sur le biplot, mais contribuent tout de même à expliquer la variation entre objets)+
- +
-Dans certains cas, la priorité n'est pas de préserver la distance exacte entre les objets, mais au contraire de représenter ​aussi fidèlement que possible les relations entre objets selon un petit nombre d'axes (généralement deux ou trois) spécifiés par l'​utilisateur. Dans de tels casle positionnement multidimensionnel non-métrique (NMDS) est la solution. Si l'​utilisateur définit un nombre d'axe égal à deuxle biplot produit par le NMDS correspond à la meilleure solution graphique ​pour représenter en deux dimensions la similarité entre objets (les objets dissimilaires étant les plus éloignées,​ et les objets similaires étant les plus proches)De plus, le NMDS permet à l'​utilisateur de choisir la mesure de distance qu'il souhaite pour ordonner les objets. +
- +
-Afin de trouver la meilleure représentation des objets, le NMDS applique une procédure itérative qui vise à positionner les objets dans le nombre spécifié de dimensions de façon à minimiser une fonction de stress (variant de 0 à 1) qui mesure la qualité de l'​ajustement de la distance entre objets dans l'​espace d'​ordination. Ainsi, plus la valeur du stress sera faible, plus la représentation des objets dans l'​espace d'​ordination sera exacte. Un second moyen d'​évaluer l'​exactitude d'un NMDS consiste à construire un diagramme de Shepard qui représente les distances entre objets sur le biplot d'​ordination en fonction de leurs distances réelles. Le R2 obtenu à partir de la régression entre ces deux types de distance mesure la qualité de l'​ajustement du NMDS. +
- +
- +
-<code rsplus | NMDS sur les données d’abondance spe avec une distance de Bray-Curtis et k=2 axes d'​ordination>​ +
-# NMDS +
-spe.nmds<​-metaMDS(spe,​ distance='​bray',​ k=2)+
    
-### Extraction des résultats +Le matériel ​de cet atelier a également été obtenu à partir des blogs et des tutoriels suivants ​
-spe.nmds +  http://www.fromthebottomoftheheap.net/blog
- +  http://www.sfs.uni-tuebingen.de/~jvanrij/Tutorial/​GAMM.html 
-### Évaluation de la qualité de l'​ajustement et construction du diagramme de Shepard +  * http://www.sfs.uni-tuebingen.de/~jvanrij/​LSA2015/​AnswersLab2.html
-spe.nmds$stress +
-stressplot(spe.nmds,​ main='​Shepard plot'​) +
- +
-# Construction du biplot +
-windows() +
-plot(spe.nmds,​ type="​none",​ main=paste('​NMDS/​Bray - Stress=',​ round(spe.nmds$stress,​ 3)), +
-     ​xlab=c("​NMDS1"​),​ +
-     ​ylab=c("​NMDS2"​)) +
-points(scores(spe.nmds,​ display="​sites",​ choices=c(1,​2)),​ +
-       ​pch=21,​ col="​black",​ bg="​steelblue",​ cex=1.2) +
-text(scores(spe.nmds,​ display="​species",​ choices=c(1)),​ +
-     ​scores(spe.nmds,​ display="​species",​ choices=c(2)),​ +
-     ​labels=rownames(scores(spe.nmds,​ display="​species"​)),​ +
-     ​col="​red",​ cex=0.8) +
-</​code>​ +
- +
-{{ :​shepard_plot.png |}} +
- +
-Le diagramme ​de Shepard identifie une forte corrélation entre les distances observées ​et les distances de l'​ordination (R2 > 0.95), et donc une bonne qualité de l'​ajustement du NMDS. +
- +
-{{ :nmds_biplot.png |}} +
- +
-Le biplot du NMDS identifie un groupe de sites caractérisés par les espèces BLA, TRU, VAI, LOC, CHA et OMB, tandis que les autres espèces caractérisent un groupe de sites situés dans le coin supérieur droit du biplot. Quatre sites situés dans le coin inférieur droit sont fortement différents des autres. +
- +
- +
-**Défi 6** +
-Exécuter un NMDS sur les données d'​abondance des //espèces d'​acariens// (données mite) en deux dimensions à partir de distances de Bray-Curtis. Évaluer la qualité de l'​ajustement et interpréter le biplot. +
- +
-**Défi 5 - Solution** +
-<​hidden>​ +
-<code rsplus>​ +
-### NMDS +
-mite.spe.nmds<-metaMDS(mite.spe, distance='​bray',​ k=2) +
- +
-### Extraction des résultats +
-mite.spe.nmds +
- +
-### Évaluation ​de la qualité de l'​ajustement +
-mite.spe.nmds$stress +
-stressplot(mite.spe.nmds,​ main='​Shepard plot'​) +
- +
-### Construction du biplot +
-windows() +
-plot(mite.spe.nmds,​ type="​none",​ main=paste('​NMDS/Bray - Stress=',​ round(mite.spe.nmds$stress,​ 3)), +
-     ​xlab=c("​NMDS1"​),​ +
-     ​ylab=c("​NMDS2"​)) +
-points(scores(mite.spe.nmds,​ display="​sites",​ choices=c(1,​2)),​ +
-       ​pch=21,​ col="​black",​ bg="​steelblue",​ cex=1.2) +
-text(scores(mite.spe.nmds,​ display="​species",​ choices=c(1)),​ +
-     ​scores(mite.spe.nmds,​ display="​species",​ choices=c(2)),​ +
-     ​labels=rownames(scores(mite.spe.nmds,​ display="​species"​)),​ +
-     ​col="​red",​ cex=0.8) +
-</code> +
- +
-{{ :​nmds_mite_shepard.png |}} +
- +
-La corrélation entre distance observée et distance d'​ordination (R2 > 0.91) et la valeur de stress relativement faible identifient une bonne qualité de l'​ajustement du NMDS. +
- +
-{{ :nmds_mite_biplot.png |}} +
- +
-Aucun groupe de sites ne peut être précisément identifié à partir du biplot, ce qui montre que la plupart des espèces sont présentes dans la plupart des sites, i.epeu de sites présentent des communautés distinctes.  +
-</hidden>​ +
- +
-=====3.5. En résumé===== +
- +
-L’ordination constitue une puissante méthode d'​analyse pour étudier les relations entre objets caractérisés par différents descripteurs (ex. des sites décrits par leurs communautés biologiques,​ ou leurs variables environnementales),​ mais de nombreuse méthodes d'​ordination existent. Ces méthodes diffèrent principalement par le type de distance qu'​elles préservent,​ le type de variables qu'​elles autorisent, et le nombre de dimensions de l'​espace d'​ordination. Pour mieux guider votre choix de la méthode d'​ordination à utiliser, le tableau ci-dessous identifie les caractéristiques de chacune des quatre méthodes d'​ordination présentées lors de cet atelier.+
  
-{{ :​resume_ordination.jpg |}}+Enfin, les pages d'​aide,​ disponibles via ''?​gam''​ dans R sont une excellente ressource.
  
-Lors du prochaine atelier, vous verrez comment identifier les relations entre variables environnementales et communautés biologiques décrivant un même ensemble de sites, à l'aide des méthodes d'​analyses canoniques.