Prezi logo
Log in
Loading…

CSBQ Atelier R 10

CQ

CSBQ QCBS

Transcript

0. Intro

1. Explorer et préparer les données

2. Analyses canoniques

2.2. RDA partielle

2.3. Partitionnement de la variation

3. Arbre de régression multivarié

4. Analyse discriminante linéaire

Sommaire

Pour débuter

Librairies requises

Intro

> install.packages("vegan")

> install.packages("labdsv")

> install.packages("plyr")

> install.packages("MASS")

> install.packages("mvpart")

> install.packages("MVPARTwrap")

> install.packages("rdaTest")

  • Cet atelier est une suite logique à l'atelier 9 mais avec plus:
  • de données environnementales!
  • de test d'hypothèses!
  • de matrices!
  • de plaisir!

Vous savez comment faire des ordinations, ce qui est vraiment utile pour l'exploration des données...

Mais que faire pour tester des hypothèses ou faire des prédictions?

Aujourd'hui...

nous explorerons des méthodes sous contrainte qui nous permettront d'intégrer les variables environnementales pour expliquer et prédire la variation dans les donnés

Recommendations:

  • créez votre propre script
  • référez-vous au code fourni si nécessaire
  • évitez de copier-coller le code fourni
  • souvenez-vous de changer votre répertoire de travail pour le dossier dans lequel vous avez enregistré les fichiers requis

Les données et le script sont disponible sur:

qcbs.ca/wiki/r/r_atelier10

Téléchargez depuis la page wiki:

  • script
  • DoubsEnv.csv
  • DoubsSpe.csv
  • données fictives
  • librairie mvpart
  • librairie MVPARTwrap
  • librairie rdaTest

Après l'installation, charger les librairies

> library(packagename)

Données environnementales

Données rivière Doubs

Données rivière Doubs

Explorer les données

decostand()

Fréquences des espèces

Transformation des données

Explorer le contenu des données environnementales:

> ?decostand

Chargez les données de poissons (DoubsSpe.csv):

Combien de zéros?

Explorer les données d'abondances de poissons:

> sum(spe==0)

> spe <- read.csv(file.choose(), row.names=1)

> spe <- spe[-8,] #enlever le site sans poisson

> (ab=table(unlist(spe)))

> barplot(ab, las=1, xlab=”Abundance class”,

ylab=”Frequency”, col=grey(5:0/5))

Il faut standardiser les variables environnementales mesurées en différentes unités pour pouvoir en comparer les effets:

> names(env) #noms des objets

> dim(env) #dimensions

> str(env) #structure des objets

> head(env) #premières lignes

> summary(env) #statistiques descriptives

Nous sommes en présence d'une grande abondance de zéros

Nous ne voulons pas considérer les absences d'espèces comme similarité entre les sites dans nos analyses

Nous appliquerons la transformation Hellinger, appropriée pour les données d'abondances d'espèces. Cette transformation peut être exécutée avec la fonction decostand() de vegan{}

408!

CENTRER-RÉDUIRE !

> env.z <- decostand(env, method="standardize")

Chargez les données de l'environnement (DoubsEnv.csv):

> names(spe) #noms des objets

> dim(spe) #dimensions

> str(spe) #structure des objects

> head(spe) #premières lignes

> summary(spe) #statistiques descriptives

Grande proportion de zéros

Quelle proportion de zéros?

données de Verneaux (1973):

  • communautés de poissons
  • 27 espèces
  • 30 sites
  • 11 variables environnementales

> env <- read.csv(file.choose(), row.names=1)

> env <- env[-8,] #enlever le site sans données

Explorer la colinéarité entre les variables en visualisant les corrélations:

La fonction centre et réduit les variables. Ces variables centrées-réduites seront utilisées pour la plupart des analyses présentées aujourd'hui.

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

> spe.hel <- decostand(spe, method="hellinger")

# N'EXÉCUTEZ QU'UNE FOIS

> pairs(env, main="Bivariate Plots of the Environmental Data" )

>50% !

> apply(env.z, 2, mean)

> apply(env.z, 2, sd)

2.1. Analyse canonique de redondances

Esthétique des graphiques

Mais qu'est-ce que ça signifie?

Analyses canoniques

triplots RDA

RDA avec R

Test statistique de la RDA

Sélection des variables

Cadrage de type 1

Sommaire de la RDA

Cadrage de type 2

Analyse canonique de redondances (RDA)

On utilise anova.cca pour tester la significativité du modèle global et des axes

Dans le type 1, les sites (points bleus) les plus similaires sont plus près les uns des autres

Il est possible de faire des triplot de RDA rapidement à l'aide de plot() mais sans ajouter d'argument supplémentaire, le résultat peut être très moyen...

Il est possible de changer plusieurs éléments et de mettre votre touche personnelle aux graphiques tant pour les points, le texte, les vecteurs...

Pour ne retenir que les variables explicatives significatives, nous utilisons la sélection progressive de variables avec la fonction ordiR2step()

  • Préparer les données
  • Faire la RDA
  • Extraire les résultats

Il y a deux options de cadrage:

  • Type 1: Les distances entre objets sont des approximations de leur similarité
  • Type 2: Les angles entre les variables reflètent leur corrélation

> anova.cca(spe.rda.signif, step=1000)

> anova.cca(spe.rda.signif, step=1000, by="axis")

  • Il s'agit d'un terme générique pour des analyses qui
  • identifient les relations entre les assemblages d'espèces de différents sites et les variables environnementales mesurées à ces sites
  • sont utilisées pour tester des hypothèses
  • ... et faire des prédictions

>points(scores(spe.rda.signif, display="sites", choices=c(1,2), scaling=1),

pch=21, col="black", bg="steelblue", cex=1.2)

> ordiR2step(rda(spe.hel~1, data=env.z),scope=formula(spe.rda),direction="forward,R2scope=TRUE,pstep=1000)

> env.signif<-subset(env.z, select=c("alt","oxy","dbo"))

> env.z<-subset(env.z, select=-das)

> spe.rda<-rda(spe.hel~.,data=env.z)

> summary(spe.rda, display=NULL)

  • Une extension de la régression multiple appliquée aux données multivariées
  • Modélise les effets d'une matrice de variables explicatives sur une matrice de variables réponses
  • Les variables peuvent être quantiatives, qualitatives ou binaires (1 ou 0)
  • Les variables explicatives devraient être centrées-réduites

Dans le type 2, les flèches et leurs angles sont très informatifs. Les variables ayant des flèches plus longues sont plus importantes pour expliquer la variation.

Un des aspects les plus puissants de la RDA est la visualisation à la fois des variables environnementales, des espèces et des sites

Faisons maintenant le triplot de la RDA...

>arrows(0,0,scores(spe.rda.signif, display="species", choices=c(1), scaling=1),scores(spe.rda.signif, display="species", choices=c(2), scaling=1), col="black",length=0)

Combien de variables sont conservées?

Quel est le R2 de la RDA effectuée avec les variables significatives seulement?

"Constrained proportion":

variation de Y expliquée par X (R2: 73.41%)

proportion non expliquée de Y: 26.59%

> spe.rda.signif<-rda(spe.hel~., data=env.signif)

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

  • les variables dbo et oxy sont très négativement corrélées
  • l'espèce TRU est très corrélé à la variable alt

Centre de la science de la biodiversité du Québec

Série d'ateliers R

Défi 1 (Solution)

Défi 1

Faire une RDA sur les données d'acariens disponibles sur vegan en fonction des variables environnementales

> data(mite)

> data(mite.env)

> mite.spe<-mite

  • Quelle est la proportion de variance expliquée par le modèle?
  • Quelles sont les variables explicatives significatives?

R2aj: 19.2%

RDA partielle sur les données Doubs

Faire le triplot de la RDA partielle

Pourquoi c'est intéressant?

RDA partielle

Défi #2

Solution #2

Comment ça fonctionne?

RDA partielle

Test statistique sur la RDA partielle

  • Faire la RDA partielle:

> mite.spe.subs <- rda(mite.spe.hel ~ Shrub + Topo

+ Condition(SubsDens + WaterCont + Substrate), data=mite.env)

#Faire la RDA partielle

> spechem.physio<-rda(spe.hel, envchem, envtopo)

> summary(spechem.physio, display=NULL)

#ou

> spechem.physio2<-rda(spe.hel ~ pH + dur + pho + nit + amm

+ oxy + dbo + Condition(alt + pen + deb), data=env.z)

> summary(spechem.physio2, display=NULL)

En gros...

  • Une RDA est effectuée pour déterminer les effets de W sur Y
  • On extrait les résidus (nouvelle matrice Yres|W)
  • contenant les variables réponses Y sans les effets de W
  • Une RDA est effectuée pour déterminer les effets de X sur Yres|W
  • Déterminer les effets de la chimie de l'eau sur les abondances d'espèces entre les sites en tenant compte des effets de la topographie
  • Un cas spécial de la RDA
  • Implique un tableau de variables explicatives additionnel
  • Permet d'isoler les effets des variables X sur les variables réponse Y en contrôlant - en ajustant - pour les effets des covariables W

Faire une RDA partielle sur les données d'acariens:

  • contrôler pour les variables de substat SubsDens, Watercont and Substrate

> R2adj<- RsquareAdj(spechem.physio)$adj.r.squared

> R2adj

> anova.cca(spechem.physio, step=1000)

> anova.cca(spechem.physio, step=1000, by="axis")

> plot(spechem.physio, scaling=2, main="Triplot partial RDA - scaling 2",

type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))

> points(scores(spechem.physio, display="sites", choices=c(1,2),

scaling=2),pch=21, col="black", bg="steelblue", cex=1.2)

> text(scores(spechem.physio, display="species", choices=c(1),scaling=2),

scores(spechem.physio, display="species", choices=c(2), scaling=2),

labels=rownames(scores(spechem.physio, display="species", scaling=2)),

col="grey", cex=0.8)

> arrows(0,0,scores(spechem.physio, display="bp", choices=c(1),scaling=2),

scores(spechem.physio, display="bp", choices=c(2), scaling=2),col="red")

> text(scores(spechem.physio, display="bp", choices=c(1), scaling=2)+0.05,

scores(spechem.physio, display="bp", choices=c(2), scaling=2)+0.05,

labels=rownames(scores(spechem.physio, display="bp", choices=c(2),

scaling=2)),col="red", cex=1)

  • Pour déterminer les effets de variables environnementales sur les compositions en espèces des sites en prenant en compte la variation due à des variables qui ne sont pas d'intérêt pour notre étude
  • Contrôle pour l'effet linéaire de variables connues
  • Isoler les effets d'un ou de plusieurs descripteur(s)
  • Analyser des échantillons appariés

> summary(mite.spe.subs, display="NULL")

> (R2adj <- RsquareAdj(mite.spe.subs)$adj.r.squared)

> anova.cca(mite.spe.subs, step=1000)

> anova.cca(mite.spe.subs, step=1000, by="axis")

> envtopo<-env.z[,c(1:3)] # topographie

> envchem<-env.z[,c(4:10)] # chimie de l'eau

  • Est-ce que le modèle est significatif?
  • Quels sont les axes significatifs?
  • Quelle est la variance expliquée par les variables de substrat?

La RDA est significative, mais seulement le premier axe canonique est significatif. Shrub et Topo expiquent 9.81% de la variation dans la communauté d'acariens et les covariables de substrat ont expliqué 42.84% de la variation.

Le R2adj de la RDA partielle est 8.33%

X

W

~

Y

Atelier 10: Analyses multivariées avancées

Partitionnement de la variation

Pour les pros...

Solution #3

Défi #3

Solution #3

Partitionnement de la variation

Diagramme de Venn

Solution #3

  • Il est possible de ne conserver que les variables significatives de chaque sous-groupe de variables...
  • ... et de tester les fractions du partitionnement qui sont significatives

> plot(spe.part.all, digits=2)

Pour séparer la variation expliquée d'une matrice réponse entre 2, 3 ou 4 matrices de variables explicatives

Faire le partitionnement de la variation des données d'acariens du aux variables significatives de substrat et aux (autres) variables significatives:

Utile pour lier les abondances d'espèces entre sites en fonction de types différents de variables environnementales

  • grande vs petite échelle
  • biotique vs abiotique
  • chimique vs physique

> #Varpart

> spe.part <- varpart(spe.hel, envchem.pars, envtopo.pars)

> plot(spe.part, digits=2)

> #Tests of significance

> anova.cca(rda(spe.hel, envchem.pars), step=1000) # [a+b]

> anova.cca(rda(spe.hel, envtopo.pars), step=1000) # [b+c]

> env.pars <- cbind(envchem.pars, envtopo.pars)

> anova.cca(rda(spe.hel, env.pars), step=1000) # [a+b+c]

> anova.cca(rda(spe.hel, envchem.pars, envtopo.pars),

step=1000) # [a]

> anova.cca(rda(spe.hel, envtopo.pars, envchem.pars),

step=1000) # [c]

envchem

envtopo

> #Varpart

> mite.spe.part <- varpart(mite.spe.hel, mite.subs.pars, mite.other.pars)

> plot(mite.spe.part, digits=2)

> #Tests of significance

> anova.cca(rda(mite.spe.hel, mite.subs.pars), step=1000) # [a+b]

> anova.cca(rda(mite.spe.hel, mite.other.pars), step=1000) # [b+c]

> mite.env.pars <- cbind(mite.subs.pars, mite.other.pars)

> anova.cca(rda(mite.spe.hel, mite.env.pars), step=1000) # [a+b+c]

> anova.cca(rda(mite.spe.hel, mite.subs.pars, mite.other.pars),

step=1000) # [a]

> anova.cca(rda(spe.hel, mite.other.pars, mite.subs.pars),

step=1000) # [c]

> mite.subs=mite.env[,c(1,2,3)] # first set of variables

> mite.other=mite.env[,c(4,5)]

> #RDA of mite.subs

> rda.mite.subs <- rda(mite.spe.hel~., data=mite.subs)

> R2a.all.chem <- RsquareAdj(rda.mite.subs)$adj.r.squared

> ordiR2step(rda(mite.spe.hel~1, data=mite.subs),

scope= formula(rda.mite.subs),direction= "forward",

R2scope=TRUE, pstep=1000)

> names(mite.subs)

> mite.subs.pars <- mite.subs[, c( 2,3 )]

> #RDA for mite.other

> rda.mite.other <- rda(mite.spe.hel~., data=mite.other)

> R2a.all.other <- RsquareAdj(rda.mite.other)$adj.r.squared

> ordiR2step(rda(mite.spe.hel~1, data=mite.other),

scope=formula(rda.mite.other), direction= "forward",

R2scope=TRUE, pstep=1000)

> names(mite.other)

> mite.other.pars <- mite.other[, c(1,2)]

> #RDA of chemistry variables

> spe.chem <- rda(spe.hel~., data=envchem)

> #Select significant chemistry variables

> R2a.all.chem <- RsquareAdj(spe.chem)$adj.r.squared

> ordiR2step(rda(spe.hel~1, data=envchem), scope= formula(spe.chem),

direction= "forward", R2scope=TRUE, pstep=1000)

> names(envchem)

> envchem.pars <- envchem[, c( 4, 6, 7 )]

> #RDA with other environmental variables

> spe.topo <- rda(spe.hel~., data=envtopo)

> R2a.all.topo <- RsquareAdj(spe.topo)$adj.r.squared

> ordiR2step(rda(spe.hel~1, data=envtopo), scope= formula(spe.topo),

direction= "forward", R2scope=TRUE, pstep=1000)

> names(envtopo)

> envtopo.pars <- envtopo[, c(1,2)]

  • Quelle proportion de la variation est expliquée par chaque sous-ensemble?
  • Quelles sont les fractions significatives?
  • Toutes les fractions sont significatives (p<0.001)

> library(vegan)

> ?varpart

> spe.part.all <- varpart(spe.hel, envchem,envtopo)

> spe.part.all

Legendre and Legendre 2012

Toutes les fractions sont significatives

Cool diagramme de Venn

> install.packages("venneuler")

> library(venneuler)

> v <- venneuler(c(A=0.241,"A&B"= 0.233, B=0.112))

> v$labels <- c("Chem\n24%","Topo \n11%")

> v$colors <- c(.78,.6)

> plot(v)

Site: http://qcbs.ca/wiki/r_atelier10

MRT sur les données Doubs

Un peu de vocabulaire

L'arbre de régression

Espèces discriminantes

Comparer les arbres

Sommaire

Validation croisée

MRT vs RDA

Arbre de régression multivarié (MRT)

valeur seuil

> summary(doubs.mrt)

Nous utiliserons la fonction mvpart() de la librairie... mvpart{}

> # En utilisant le critère CVRE (10 banches)

> doubs.mrt.cvre<- mvpart((data.matrix(spe.hel)~., env, legend=FALSE,

margin=0.01, cp=0, xv="pick", xval=nrow(spe.hel),xvmult=100,which=4)

> # En choisissant de faire un arbre à 4 branches

> doubs.mrt.4<- mvpart((data.matrix(spe.hel)~., env, legend=FALSE,

margin=0.01, cp=0, xv="pick", xval=nrow(spe.hel),xvmult=100,which=4)

Est-ce que certaines espèces contribuent beaucoup à la variation entre les sites à chaque division?

Quels sont les sites inclus dans chaque feuille?

Nous pouvons trouver les réponses à ces questions en générant un sommaire plus détaillé à l'aide de MRT().

Cette fonction est disponible sous MVPARTwrap{} et fait également appel à rdaTest{}

  • Feuille: groupe terminal de sites
  • Noeud: point où les donées se divisent
  • caractérisé par une valeur seuil d'une variable explicative
  • Branche: chaque lignée formée par un noeud

> ?mvpart

> # Enlever "distance from source" des données environnementales brutes

> env <- subset(env, select = -das)

> # Créer l'arbre

> doubs.mrt<- mvpart((data.matrix(spe.hel)~., env, legend=FALSE,

margin=0.01, cp=0, xv="pick", xval=nrow(spe.hel),xvmult=100,which=4)

  • Une technique de groupement sous contrainte qui:
  • fait le partitionnement d'une matrice réponse quantitative sous la contrainte d'une matrice de variables explicatives
  • implique deux procédures s'effectuant en parallèle:
  • le partitionnement sous contrainte (ou en d'autre mots, la construction de l'arbre)
  • la sélection du partitionnement optimal par validation croisée
  • résulte en un arbre
  • diagrammes à bandes des abondances d'espèces aux sites inclus dans cette feuille
  • erreur relative
  • nombre de sites

Avantages du MRT:

  • n'assume pas de relation linéaire entre les matrices Y et X
  • est robuste en présence de valeurs manquantes
  • est robuste en présence de colinéarité entre descripteurs
  • ne nécessite pas la transformation des variables explicatives
  • l'arbre est facile à interpréter

Le code ne change absolument pas...

La sélection de la taille de l'arbre s'effectue sur le graphique de validation croisée

> doubs.mrt.wrap<-MRT(doubs.mrt, percent=10,

species=colnames(spe.hel))

> summary(doubs.mvpart.wrap)

xv= "pick" # pour choisir où "élaguer" (couper) l'arbre

xval=now(spe.hel) # nombre de validations croisées

xvmult=100 # nombre de permutations par validation croisée

  • Les espèces discriminantes principales sont TRU, VAI et ABL.
  • ABL est l'espèce la plus indicatrice des sites à basse altitude

CP veut dire "complexity parameter". En d'autres mots, la variance expliquée à chaque noeud. La première ligne, à zéro noeuds, c'est le R2 de l'arbre global.

la réciproque de l'erreur (le R2): 43.7%

Solution #4

Défi #4

Créer un arbre de régression multivarié avec les données sur les acariens.

> mite.mrt<-mvpart(data.matrix(mite.spe.hel)~.,

mite.env,legend=FALSE, margin=0.01, cp=0,

xv="pick", xval=nrow(mite.spe.hel), xvmult=100,which=4)

> summary(mite.mrt)

> mite.mrt.wrap<-MRT(mite.mrt, percent=10,

species=colnames(mite.spe.hel))

> summary(mite.mrt.wrap)

  • Selectionner l'arbre de taille minimale à une erreur type du CVRE minimal
  • Quelle est la proportion de variation expliquée par l'arbre?
  • Combien de branches contient-il?
  • Quelles sont les espèces discriminantes?
  • 25.6% de la variation dans la composition de la communauté d'acariens entre les sites est expliquée par le partitionnement des sites à une valeur seuil de contenu en eau du substrat de 385.1 mg/l.
  • LCIL est une espèce discriminante des sites à haute teneur en eau

Groupements prédits

Défi 5: solution

Défi 5

Défi 5 solution

Faire des prédictions

Groupement a priori

Faire la LDA

sommaire LDA

LDA

Le sommaire de la LDA peut être utilisé pour déterminer dans quel groupe la LDA a classifié chaque objet (site) et indique la précision du groupement

Nous utiliserons le facteur de groupement de la colonne 4 du nouveau tableau spa.group

  • Pour notre exemple, nous ferons un groupement a priori pour les données de poissons de la rivière Doubs
  • Nous choisirons la latitude puisque nous avons les coordonnées x et y de chaque site dans l'espace (DoubsSpa.csv)
  • Nous séparerons l'étendue des latitudes en 3 classes

Faire une LDA pour les données environnementales d'acariens (deux premières variables) basé sur 4 groupes latitudinaux créés à partir des données mite.xy

>LDA.mite<-lda(mite.env[,1:2],mite.xy.group[,4])

>mite.class <- predict(LDA.mite)$class

>mite.post <- predict(LDA.mite)$posterior

>mite.table <- table(mite.xy.group[,4], mite.class)

>diag(prop.table(mite.table, 1))

#read in new sites

>classify.me<-read.csv("classifyme.csv", header = T)

#predict grouping of new data

>predict.group<-predict(LDA, newdata=classify.me)

#give classification for each new site

>group.new<-predict.group$class

Maintenant que nous avons la relation entre les variables environnementales et les groupements, nous pouvons tenter de classer de nouveaux sites dans les groupes en fonction de cette relation

Nous utiliserons les résultats de la LDA pour prédire le groupement des sites du tableau ClassifyMe

>mite.xy$site<-seq(1:70)

(max(mite.xy[,2])-min(mite.xy[,2]))/4

>mite.xy.group<-ddply(.data=mite.xy, .variables=.(x, y, site), .fun= summarise, group = if(y <= 2.5) 1 else if (y <= 4.9) 2 else if (y <= 7.3) 3 else 4)

>mite.xy.group<-mite.xy.group[with(mite.xy.group, order(site)), ]

  • L'analyse discriminante linéaire permet de déterminer la quantité de variation expliquée par un groupement a priori (ex. par latitude ou par écotype)

  • Peut être utilisé pour faire des prédictions sur la classification d'une nouvelle observation dans un des groupes

  • On pourrait tenter de classer un poisson dans une population de lacs ou d'océan en fonction de variables collectées sur ces deux types de populations

#classification des objets par la LDA

>spe.class <- predict(LDA)$class

#probabilités a posteriori que les objets appartiennent aux groupes

>spe.post <- predict(LDA)$posterior

data(mite.xy)

#run LDA

>LDA<-lda(env,spa.group[,4])

#ajout des numéros de sites

>numbers<-(1:30)

>numbers<-numbers[!numbers%in%8]

>spa$site<-numbers

#faire les groupes en fonction de la latitude y<82=group1, 82<y<156=group2, y>156=group3

>spa.group<-ddply(.data=spa, .variables=.(x, y, site), .fun= summarise, group = if(y <= 82) 1 else if (y <= 156) 2 else 3)

#ordonner par site

>spa.group<-spa.group[with(spa.group, order(site)), ]

L'ordre du tableau est conservé. Ainsi, le site 31 est dans le groupe 1, le site 32 est éaglement dans le groupe 1 et ainsi de suite...

La syntaxe de la LDA est:

(variable à utiliser pour expliquer le groupement, variables de groupement)

LDA output