3.3. Analyse en Coordonnées Principales
Sommaire
0. Introduction
Suivez le lien
Librairies nécessaires
Intro
Recommendations
Ordination en espace réduit
4, 5, 6 + Dimensions
Exemple: 1 Dimension
2 Dimensions
3 Dimensions
> install.packages("vegan")
> install.packages("gclus")
> install.packages("ape")
> require(vegan)
> require(gclus)
> require(ape)
> source(file.choose()) #coldiss.R
Ordination?
...
Tout le monde en fait, mais qu'est-ce que ça signifie?
- Avez-vous déjà "ordonné"?
- Quel type de données?
- Communautés?
- Génétiques?
- Quel logiciel?
- Canoco?
- R? vegan? ecodist? ape?
- créez votre propre script
- référez-vous au code fourni uniquement en cas de besoin
- éviter de copier/coller le code à partir du script
Données et code R disponible à:
qcbs.ca/wiki/r/r_atelier9
Naviguez dans la page wiki et téléchargez:
- Script R
- DoubsEnv.csv
- DoubsSpe.csv
- coldiss.R
# N'oubliez pas le chemin!
- L'algèbre matricielle est complexe et difficile à comprendre
- Une compréhension générale est suffisante pour utiliser efficacement les méthodes d'ordination
Mais que se passe-t-il si nous voulons nous intéresser à la réponse de différentes espèces d'algues?
1. Exploration des données
Données environnementales
Fréquences des espèces
Données environnementales
Exploration des données
decostand()
Richesse en espèces
Fréquences des espèces
Données de poissons de la rivière Doubs
Observer le nombre d'espèces présentes dans chaque site:
Explorer le contenu des données environnementales:
> ?decostand
Observer la distribution de fréquence des espèces:
Chargement des données espèces (DoubsSpe.csv):
Explorer le contenu des données espèces:
> sum(spe==0)
> site.pres <- rowSums(spe>0)
> barplot(site.pres, main="Species richness",
las=1, xlab="Sites",
ylab="Number of species", col="grey")
> (ab=table(unlist(spe)))
> barplot(ab, las=1, xlab=”Abundance class”,
ylab=”Frequency”, col=grey(5:0/5))
Standardiser les variables environnementales est indispensable car il est impossible de comparer des variables d'unités différentes:
> spe <- read.csv(file.choose(), row.names=1)
> spe <- spe[-8,] #supprimer le site vide
> names(env) #noms des objets
> dim(env) #dimensions
> str(env) #structure des objets
> head(env) #5 premières lignes
> summary(env) #résumé statistique
STANDARDISATION NÉCESSAIRE !
> env.z <- decostand(env, method="standardize")
> names(spe) #noms des objets
> dim(spe) #dimensions
> str(spe) #structure des objets
> head(spe) #5 premières lignes
> summary(spe) #résumé statistique
Chargement des données environnementales (DoubsEnv.csv):
Quelle proportion de zéros?
Explorer la colinéarité en vérifiant les corrélations entre variables:
Données de Verneaux (1973) :
- caractérisation des communautés de poissons
- 27 espèces
- 30 sites
- 11 variables environnementales
Cette fonction centre-réduit les données pour permettre la fiabilité des analyses suivantes:
> sum(spe==0)/(nrow(spe)*ncol(spe))
> env <- read.csv(file.choose(), row.names=1)
> env <- env[-8,] #supprimer le site vide
> pairs(env, main="Bivariate Plots of the Environmental Data" )
> apply(env.z, 2, mean)
> apply(env.z, 2, sd)
ATTENTION: N'EXÉCUTER QU'UNE SEULE FOIS!
2. Similarité/Dissimilarité/Groupement
Au-délà de la 1ère dimension
Similarité/Dissimilarité
Défi #1
Solution #1
2.1 Mesure d'association
Comparaison des sites
ATTENTION!
2.2 Transformation des données de communautés
Transformation des données de communautés
Mesures de distance
des communautés
Transformation des données de communautés
Comparaison des sites de la rivière Doubs
Visualisation d'une matrice de distances
- L'algébre matricielle est au coeur de plusieurs méthodes d'analyses multivariées
- Utile pour comprendre vos données
- Certains types d'ordination ou de groupement nécessitent des mesures appropriées
- Transformer des abondances en présence-absence
- Euclidienne
- Manhattan
- Corde
- Hellinger
- Chi-carré
- Bray-Curtis
La fonction vegdist() comprend les mesures de distances communes:
Discuter avec votre voisin: Comment savoir si deux objets caractérisés par des données multidimensionnelles sont similaires?
Faites une liste de vos suggestions.
> spe.pa <- decostand(spe, method="pa")
Il faut comprendre les forces et limites de vos données de communautés avant d'effectuer une ordination, et standardiser/transformer vos données en conséquence:
Le choix de la mesure de distance
nécessite une bonne compréhension de vos données!
Similarité: S = 1 - D
Distance: D = 1 - S
> ?vegdist
- Réduire le poids des espèces rares:
# Chaque mesure est utile dans différentes situations
- Les jeux de données écologiques correspondent souvent à de grandes matrices
- L'ordination calcule les relations entre espèces, ou entre objets
- Ces relations peuvent être simplifiées par des mesure de dissimilarité
> spe.hel <- decostand(spe, method="hellinger")
> spe.chi <- decostand(spe, method="chi.square")
Comment la composition des communautés diffère-t-elle entre les 30 sites de la rivière Doubs?
- abondances/comptages/présence-absence?
- distributions asymétriques?
- beaucoup d'espèces rares?
- sur-abondance d'espèces dominantes?
- problème de double-zéros?
- Réduire le poids des espèces abondantes:
- Explorer différentes mesures de distance entre objets permet de mieux comprendre le fonctionnement de l'ordination
> spe.db.pa <- vegdist(spe,method="bray")
Les mesures de distance sont très utiles pour comprendre vos données, mais l'ordination préserve plus d'information.
> spe.pa <- decostand(spe, method="log")
Centre de la science de la biodiversité du Québec
Série d'ateliers R
Comparaison
Groupement de Ward
Groupement à liens complet
Groupement
Groupement à liens simples
2.3 Groupement
Groupement hiérarchique
Groupement de Ward
Quelle méthode choisir?
Comparaison
Aperçu de 3 méthodes hiérarchiques
- Faire le groupement de Ward:
- À partir d'une matrice de distances, on classe les objets en ordre croissant
- Créer une matrice de distances à partir des données de la rivière Doubs transformées Hellinger
> spe.dhel.ward <- hclust(spe.dhel, method="ward.D2")
- les deux objets les plus proches se regroupent
- ensuite les groupes se lient à la distance à laquelle les objets qu'ils contiennent sont tous liés
- les deux objets les plus proches se regroupent
- ensuite les deux objets les plus proches suivants
- et ainsi de suite.
- Les éléments de petits ensembles se regroupent en groups plus vastes de rang supérieur
- ex. espèces, genres, familles, ordres...
> spe.dhel <- vegdist(spe.hel, method="euclidean")
- Groupement agglomératif à liens simples
- Groupement agglomératif à liens complets
- Groupement de Ward
- Permet de mettre en lumière des structures dans les données en partitionnant les objets
- Représenté sous forme de dendrogramme (arbre)
- Pas une méthode statistique!
- Dessiner le dendrogramme en utilisant la racine carrée des distances:
- Dépend de votre objectif
- démontrer des gradients? des contrastes?
- Si plus d'une méthode semble adéquate --> comparer les dendrogrammes
- Encore une fois: ceci n'est pas une analyse statistique
- Mais! il est possible:
- de déterminer le nombre de groupes optimal
- de faire des tests statistiques sur les résultats
- de combiner le groupement à l'ordination pour distinguer des groupes de sites
- Faire le groupement à liens simples:
- Le critère pour lier les groupes: la méthode des moindres carrés
- les groupes fusionnent de façon à minimiser la variance intragroupe
- à 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
> spe.dhel.single <- hclust(spe.dhel, method="single")
> plot(spe.dhel.single)
> spe.dhel.ward$height <- sqrt(spe.dhel.ward$height)
> plot(spe.dhel.ward, hang=-1)
Liens complets
Les groupes sont plus distincts
- Faire le groupement à liens complets:
Liens simples
Les objets ont tendance à s'enchaîner (ex. 19, 29, 30, 20, 26, etc.)
hang=-1 permet d'afficher les objets sur la même ligne (le bout des groupes est à distance 0)
Les objets ont tendance à former des groupes plus sphériques et homogènes
> spe.dhel.complete <- hclust(spe.dhel, method="complete")
> plot(spe.dhel.complete)
Atelier 9:
Analyses multivariées
Ordination non-contraintes
3. Ordination non contraintes
- Évaluer les relations au sein d'un jeu de données (et non entre jeu de données, i.e. analyses canoniques).
- Trouver les composantes clefs de la variation entre échantillons, sites, espèces...
- Réduire le nombre de dimensions des données multivariables sans perte importante d'information.
- Créer de nouvelles variables utilisables dans d'autres analyses (e.g. régression).
Site web: http://qcbs.ca/wiki/r_atelier9
3.1. Analyse en Composantes Principales
PCA - Ce qu'il vous faut
PCA - Principes
PCA - Interprétation des sorties R
PCA sur les données Poissons!
PCA - Interprétation des sorties R
PCA - Cas multidimensionnel
PCA - Principes
PCA - Interprétation des sorties R
Analyse en Composantes Principales (ACP ou PCA)
PCA - Illustration graphique
- Un jeu de données correspondant à des variables réponses (ex. composition de communautés) OU à des variables explicatives (ex. variables environnementales)
- La PCA (tout comme la RDA) est implémentée par la fonction rda() de la librairie vegan:
L'abondante information produite par une PCA est plus facile à comprendre et à interpréter à l'aide de biplots permettant de visualiser les patrons présents dans les données.
> ?rda
- PC1 --> axe qui maximise la variance des points projetés perpendiculairement sur les axes.
- PC2 --> doit être orthogonal à PC1, mais sa direction maximise la variance des points projetés.
- PC3 --> et ainsi de suite: orthogonale à PC1 et PC2...
- Liste des valeurs propres associées à chaque Composantes Principales (ici 27 PCs sont identifiées, soit le nombre de dimensions des données)
- Effectuer une PCA sur les abondances de poissons transformées Hellinger:
- Un biplot peut être rapidement créé via la fonction "plot()":
> spe.h.pca <- rda(spe.hel)
- Liste de la proportion de variance expliquée par chaque Composante Principale (et la proportion cumulée)
- Échantillons correspondant à des mesures du même jeu de variables
- Généralement, un jeu de données plus long que large est préféré
- Préserve, en 2 dimensions, le maximum de variation des données
- Il en résulte des variables synthétiques orthogonales entre elles (et donc non corrélées)
*La valeur propre est la valeur du changement dans la longueur d'un vecteur, et ici représente la quantité de variation expliquée par chaque Composante Principale.
> plot(spe.h.pca)
- Total de variance expliquée par les descripteurs (ici les espèces de poissons)
- Dans une PCA, la proportion "Total" et "Unconstrained" de variance expliquée est identique
Quand il y a plus de deux dimensions, la PCA produit
un nouvel espace dans lequel tous les axes sont orthogonaux (i.e. la corrélation entre les axes =0) et où les axes sont ordonnés selon le pourcentage de variation des données brutes qu'ils expliquent (valeur propre).
Le graphique final subit une rotation afin que les deux axes correspondent aux composantes principales
(et non plus aux espèces)
51.3% de 0.5025 égal 0.258
En 2D, les sites seraient disposés de cette façon
0.258 + 0.064 + .... = 0.5025 Variance totale expliquée
> summary(spe.h.pca)
la première composante principale est celle qui maximise la variation observée, i.e. la meilleure droite entre les sites
La seconde composante principale est ajouté perpendiculairement au premier axe
QUIZ !!!
Sélection des PCs significatives
Accéder à une partie de la sortie R
Critère de Kaiser-Guttman
PCA - Variables environnementales
PCA - Interprétation des sorties R
PCA - Variables environnementales
PCA - Interprétation des sorties R
- Extraire les valeurs propres de chaque PC:
Sélectionner les PCs qui capturent plus de variance que la moyenne de tous les PCs
Un PCA peut aussi être effectué sur les variables environnementales standardisées pour comparer les sites ou évaluer les corrélations entre variables,...
Vous pouvez extraire les scores (des sites ou des espèces) pour réaliser des graphiques ou les utiliser dans de nouvelles analyses:
La sortie R est très dense, mais vous pouvez accéder au besoin à des informations spécifiques:
> ev <- env.pca$CA$eig
La force de la PCA est de condenser la variance contenue dans un grand jeu de données en jeu de données synthétiques moins nombreuses.
- critère de Kaiser-Guttman
- modèle Broken Stick
- etc...
- Extraire les valeurs propres de chaque PCs:
- Sélectionner les valeurs propres supérieures à la moyenne:
- Extraire le score des espèces pour PC1 et PC2:
- Extraire les valeurs propres et leur contribution à la variance expliquée:
- Effectuer une PCA sur les variables environnementales standardisées
> ev <- spe.h.pca$CA$eig
> ev[ev>mean(ev)]
> spe.scores <- scores(spe.h.pca,
display="species",
choices=c(1,2))
> summary(spe.h.pca, display=NULL)
> env.pca <- rda(env.z)
Dans notre cas, 27 PCs sont identifiées, mais seuls les premières contribuent de façon importante à la variance expliquée tandis que les autres représentent le bruit des données et peuvent être écartés...
- Sélectionner les valeurs propres supérieures à la moyenne:
- Species refère aux colonnes de votre jeu de données, ici différentes espèces de poissons
- Les scores correspondent aux coordonnées de chaque espèce le long de chaque PC.
- Il existe deux façons principales de représenter une ordination en 2D...ici la sortie R nous informe que le cadrage utilisé par défaut est de type 2...
- Extraire le score des sites pour PC1 et PC2::
- Calculer les valeurs propres:
> ev[ev>mean(ev)]
- Site réfère aux lignes de votre jeu de données, ici différentes stations d'échantillonage le long de la rivière
- Les scores correspondent aux coordonnées de chaque site le long de chaque PC
> 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")
> eigen(cov(spe.hel))
> summary(env.pca)
> summary(env.pca, scaling=2) #default
Comment choisir les PCs significatives?
> site.scores <- scores(spe.h.pca,
display="sites",
choices=c(1,2))
Biplot de PCA avec biplot()
Biplot de PCA avec plot()
Type de scaling
Cadrage 2 avec biplot()
Biplot avancés
Cadrage 1 avec biplot()
ggvegan
Solution #3
rgl and vegan3d
Solution #3
Défi #3
- Biplot 3D interactif avec rgl
- Un ensemble d'outils pour créer des biplots avec gglot2:
- Transformation Hellinger des données:
- En extrayant certains parties de la sorties R, il est possible de créer des biplots plus détaillés et clairs:
- Avec la fonction biplot() de la librairie base, des flèches sont représentées pour montrer les directions et les angles entre descripteurs:
> mite.spe.hel <- decostand(mite.spe, method="hellinger")
- Cadrage 1: préserve au maximum la distance euclidienne (dans l'espace d'ordination) entre objets (ex. sites); les angles entre descripteurs (ex. espèces) ne sont pas informatifs. Meilleur cadrage pour interpréter les relations entre objets!
> require(rgl)
> require(vegan3d)
> ordirgl(spe.h.pca)
Exécuter une PCA sur les données d'abondance d'acariens :
> biplot(spe.h.pca)
> data(mite)
> mite.spe.h.pca <- rda(mite.spe.hel)
> install.packages("devtools")
> require("devtools")
> install_github("ggvegan", "gavinsimpson")
> require(ggvegan)
> autoplot()
C'est l'heure du quiz !
> plot(spe.h.pca, scaling=1, type="none",
xlab = c("PC1 (%)", round((spe.h.pca$CA$eig[1]/sum(spe.h.pca$CAeig))*100,2)),
ylab = c("PC2 (%)", round((spe.h.pca$CA$eig[2]/sum(spe.h.pca$CAeig))*100,2)))
> points(scores(spe.h.pca, display="sites", choices=c(1,2), scaling=1),
pch=21, col="black", bg="steelblue", cex=1.2)
> text(scores(spe.h.pca, display="species", choices=c(1), scaling=1),
scores(spe.h.pca, display="species", choices=c(2), scaling=1),
labels=rownames(scores(spe.h.pca, display="species", scaling=1)),
col="red", cex=0.8)
- Recherche des axes significatifs par critère de Guttman-Kaiser:
- Les descripteurs distants de 180 degrés sont négativement corrélés
- Les descripteurs distants de 90 degrés ne sont pas corrélés
- Les descripteurs distants de 0 degrés sont positivement corrélés
- Par défaut, le cadrage est de type 2
- La fonction plot() est rapide, mais le graphique qui en résulte permet difficilement d'interpréter les angles entre descripteurs
1. Que signifie PCA?
a) Principal Coordinate Analysis
b) Principle Compound Analysis
c) Principal Component Analysis
d) Principle Component Analysis
- Cadrage 2 (option par défaut): les distances entre objets ne sont pas des approximations de leurs distances euclidiennes, mais les angles entre descripteurs reflètent leurs corrélations. Meilleur cadrage pour interpréter les relations entre descripteurs !
- La distance entre sites reflètent leurs différences de composition en poissons
# voir la fonction arrows() pour ajouter des vecteurs
> ev <- mite.spe.h.pca$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")
> spe.sc<-scores(spe.h.pca,choices=1:2,scaling=,display="sp")
> arrows(0,0,spe.sc[,1],spe.sc[,2],length=0)
- Quels sont les axes significatifs?
- Quels groupes de sites pouvez-vous identifier?
- Quelles espèces caractérisent ces groupes de sites?
ATTENTION !
- La PCA est une méthode linéaire basée sur quelques hypothèse clefs:
- distribution multinormale des données
- nombre limité de zéros
- le gradient d'intérêt doit causer la majorité de la variance dans le jeu de données
- Certains de ces problèmes peuvent être réglés en utilisant des transformations appropriées des données avant d'effectuer une PCA
- Dans certains cas, tels que les études couvrant de longs gradients environnementaux, il est préférable d'utiliser d'autres méthodes d'ordination non-contraintes (ex. CA)
C'est l'heure du quiz !
# Le non-respect de ces hypothèses peut causer une forme de fer à cheval sur les biplots (horseshoe effect), sur lesquels les extrémités du fer à cheval sont proches mais représentent en réalité des conditions opposées du gradient
2. Laquelle de ces méthodes est la meilleure pour visualiser les distances entre composition des communautés de différents sites?
a) Principal Component Analysis
b) Principal Coordinate Analysis
c) Principal Response Curves
d) Redundancy Analysis
3.2. Analyse des Correspondances
C'est l'heure du quiz !
3. Que représente une valeur propre dans une PCA ?
a) Une mesure de significativité
b) La proportion de variance expliquée par une composante principale
c) Les coordonnées d'un objet dans l'espace d'ordination
d) Le secret de la vie, de l'univers, et de pas mal tout le reste
C'est l'heure du quiz !
Données non-centrées, beurk !
Les 2 premiers axes
expliquent 100% de la variation
Défi #4
Solution #4
CA: axes significatifs
CA: sorties R
Comment effectuer une CA ?
CA: biplots
Principes de la CA
CA: interprétation des résultats
CA: interprétation des biplots
Solution #4
Distances euclidiennes vs distances de Chi2
- Comme pour une PCA, le critère de Kaiser-Guttman identifie les axes significatifs:
- Les résultats d'une CA sont présentés de la même façon que pour une PCA et peuvent être appelés en utilisant:
Exécuter une CA sur les données d'abondance des espèces d'acariens.
> mite.spe.ca <- rda(mite.spe)
- La CA est implémentée dans la librairie vegan par la fonction cca:
- Dans de tels cas, la CA devrait être préférée à la PCA car elle préserve les distances de Chi2 entre objets... et représente donc mieux les réponses unimodales
> data(mite)
- Recherche des axes significatifs par critère de Guttman-Kaiser:
- Le critère de Kaiser-Guttman montre que les 5 premiers axes contribuent majoritairement à la variance expliquée.
> spe.ca <- cca(spe)
> summary(spe.ca)
- La PCA préserves les distances euclidiennes entre objets et postule une relation linéaire entre descripteurs
- ...mais dans certains cas (ex. gradients longs), les espèces présentent une réponse unimodale aux gradients environnementaux
- 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érise par les espèces LOC, VAI et TRU; le dernier groupe de sites dans le coin inférieur droit montre des communautés abondantes en BLA, CHA et OMB.
Live Long and Ordinate
> 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")
> ev <- mite.spe.ca$CA$eig
> ev[ev>mean(ev)]
> n <- length(ev)
> barplot(ev, main="Eigenvalues", col="grey", las=2)
> abline(h=mean(ev), col="red")
> legend("topright", "Average eigenvalue",
lwd=1,col=2, bty="n")
- Ici, une CA est effectuée sur les abondances de poissons
- Quels sont les axes importants?
- Quels groupes de sites se distinguent?
- Quelles espèces caractérisent à chaque groupe de sites?
PCoA et distances non-métriques
Biplot de PCoA avec biplot.pcoa()
PCoA et distances non-métriques
Solution #5
Solution #4
Défi #5
PCoA - Interprétation des sorties R
PCoA
PCoA - Interprétation des sorties R
PCoA - Données Poissons
Solution #5
- Construisons maintenant le biplot (sans espèces):
- La PCoA peut aussi être utilisée pour capturer de l'information à partir de distance non-métriques, telles que la distance de Bray-Curtis:
- La fonction biplot.pcoa() permet de visualiser en 2D les distances entre sites, et les espèces associées à chaque site:
- Les fonctions cmdscale() et pcoa() des librairies {stats} et {ape} permettent d'effectuer une PCoA:
> biplot.pcoa(mite.spe.h.pcoa, mite.spe.hel)
- Transformation Hellinger:
> biplot.pcoa(spe.bray.pcoa)
> spe.bray.pcoa <- pcoa(spe.db)
> spe.bray.pcoa
> biplot.pcoa(spe.h.pcoa, spe.hel)
Exécuter une PCoA sur les données d'abondances des espèces d'acariens transformées Hellinger.
> mite.spe.hel <- decostand(mite.spe, method="hellinger")
> ?cmdscale
> ?pcoa
- Effectuer une PCoA sur les abondances de poissons transformées Hellinger:
> mite.spe.h.pcoa<-pcoa(dist(mite.spe.hel))
- Vecteurs: vecteurs propres associés à chaque valeur propre contenant les coordonnées de chaque site dans l'espace euclidien.
- Observer la sortie R et noter la présence de valeurs propres négatives. Elles sont liées à l'impossibilité de représenter des distances non-métriques dans un espace euclidien sans corrections (cf. Legendre & Legendre 2012):
> spe.h.pcoa <- pcoa(dist(spe.hel))
- En PCA, le maximum de la variation des données est préservée
- En PCoA, la distance entre objets est préservée autant que possible dans un espace multidimensionnel
- 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?
- Valeurs propres
- Valeurs propres relatives
- Modèle broken stick: évalue les axes significatifs
- Valeurs propres cumulées: cumul des valeurs propres relatives
- Broken stick cumulés: cumul des valeurs du modèle broken stick
# Ce sont les résultats les plus utiles pour des analyses subséquentes puisqu'ils capturent fidèlement la distance entre objets.
> spe.bray.pcoa <- pcoa(spe.db,
correction="cailliez")
> spe.bray.pcoa
> spe.h.pcoa
# LA PCoA EST PARTICULIÈREMENT CONSEILLÉE POUR DES JEUX DE DONNÉES PLUS LARGES QUE LONGS
3.4. Positionnement multidimensionnel non-métrique
Principes du nMDS
nMDS
Solution #6
Défi #6
Solution #6
Qualité de l'ajustement
NMDS - Données Poissons
- La fonction metaMDS de la librairie vegan permet de réaliser un nMDS:
> spe.nmds<-metaMDS(spe, distance='bray', k=2)
> mite.spe.nmds<-metaMDS(mite.spe, distance='bray', k=2)
- nMDS
- équivalent non-métrique de la PCoA
- basé sur un algorithme itératif d'optimisation pour identifier la meilleure représentation possible des données dans l'espace d'ordination
- de plus en plus populaire
- Dans un nMDS, l'utilisateur peut ainsi spécifier:
- le nombre de dimensions
- la mesure de distance
- En PCA, CA et PCoA, les objets sont ordonnés dans un petit nombre de dimensions (i.e. axes) généralement > 2
- En conséquence, les biplots 2D ne représentent pas toute la variation présente dans les données.
- Parfois, l'objectif est cependant de représenter les données dans un nombre défini de dimensions.
- Comment effectuer une ordination pour illustrer l'ensemble de la variation des données ?
> spe.nmds
- 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.
> spe.nmds <- (spe, distance='bray', k=2)
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.
- Le biplot du nMDS identifie un groupe de sites caractérisées 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.
- Le diagramme de Shepard 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.
> 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)
- Évaluation de la qualité de l'ajustement
Exécuter un nMDS sur les données d'abondance des espèces d'acariens en deux dimensions à partir de distances de Bray-Curtis.
Évaluer la qualité de l'ajustement et interpréter le biplot.
- Évaluer la qualité de l'ajustement
> mite.spe.nmds$stress
> stressplot(mite.spe.nmds, main='Shepard plot')>
- distance spécifie la mesure de distance choisie
- k spécifie le nombre de dimensions
> spe.nmds$stress
> stressplot(spe.nmds, main='Shepard plot')
Solution #6
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.
Conclusion
- Beaucoup de méthodes d'ordination existent, mais leurs spécificités doivent guider le choix de la méthode à utiliser: