# Noms: Amanda Winegardner, Xavier Giroux-Bougard and Berenger Bourgeois  
# Date: Janvier 2015
# Description: Multivariate analysis
# Dataset: Fichiers "DoubsSpe.csv" et "DoubsEnv.csv"
# Notes: Matériel du script R fournis à partir de
#        Borcard D., Gillet F. et Legendre P., 2011. Numerical Ecology with R. Springer.
#************************************************************************************#


# 0. Chargez les librairies et les données ####

rm(list=ls())
install.packages("vegan")
install.packages("gclus")
install.packages("ape")
library(vegan)
library(gclus)
library(ape)
source(file.choose()) # fichier coldiss.R

# Matrice d'abondances d'espèces (abondances de poissons): "DoubsSpe.csv"
spe <- read.csv(file.choose(), row.names=1)
spe <- spe[-8,] # Retirer le site 8 car il ne contient aucun poisson.

# Données environnementales: "DoubsEnv.csv"
env <- read.csv(file.choose(), row.names=1)
env <- env[-8,]

#************************************************************************************#

# 1.Exploration des données ####

## 1.1 Données d'abondances d'espèces ####
names(spe)
dim(spe)
str(spe)
head(spe)
summary(spe) 

### Distribution des espèces dans les sites
(ab <- table(unlist(spe)))
barplot(ab, las=1, xlab="Abundance class", ylab="Frequency", col=grey(5:0/5))

### Nombre d'absences
sum(spe==0)

### Proportion de zéros dans les données de la communauté de poissons
sum(spe==0)/(nrow(spe)*ncol(spe))

### Co-occurrence des espèces
spe.pres <- colSums(spe>0) # Somme des sites où chaque espèce est présente
hist(spe.pres, main="Co-occurence des espèces", las=1, xlab="Fréquence", ylab="Nombre d'espèces", breaks=seq(0,30, by=5), col="grey")

### Richesse en espèces
site.pres <- rowSums(spe>0) # compute the number of species at ecah site
hist(site.pres, main="Richesse en espèces", las=1, xlab="Fréquence", ylab="Nombre d'espèces", breaks=seq(0,30, by=5), col="grey")

### Calcul d'indices de diversité
?diversity
(H <- diversity(spe, index="shannon")) # Indice de Shannon 
(N <- diversity(spe, index="simpson")) # Indice de Simpson 

## 1.2 Données environnementales ####
names(env)
dim(env)
str(env)
head(env)
summary(env)
pairs(env, main="Corrélations entre les variables environnementales" )

# Standardisation (centrage-réduction) des 11 variables environnementales 
env.z <- decostand(env, method="standardize")
apply(env.z, 2, mean) # Les données sont maintenant centrées (moyennes~0)
apply(env.z, 2, sd)   # et réduites (écart-types=1)

# Note: les variables environnementales mesurées en diffférentes unités doivent être
#  		standardisées avant le calcul de mesures de distances et les ordinations.

#***********************************************************************************#

# 2.Mesures d'association ####

## 2.1. Mesures de distances ####

### Données quantitatives des espèces
?vegdist # cette fonction calcule des indices de dissimilarité pour les données de communauté
(spe.db<-vegdist(spe, method="bray")) # distance de Bray-Curtis 
(spe.dj<-vegdist(spe, method="jac")) # distance de Jaccard
(spe.dg<-vegdist(spe, method="gower")) # distance de Gower

### Pour données de présence-absences
(spe.db.pa<-vegdist(spe, method="bray", binary=TRUE)) # distance de Bray-Curtis 
(spe.dj.pa<-vegdist(spe, method="jac", binary=TRUE)) # distance de Jaccard
(spe.dg.pa<-vegdist(spe, method="gower", binary=TRUE)) # distance de Gower

### Représentations graphiques des matrices d'association: carte des points chauds
windows() # utilisez x11() ou quartz() pour mac
coldiss(spe.db, byrank=FALSE, diag=TRUE) # Carte des points chauds - Bray-Curtis dissimilarity
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  

### Données environnementales quantitatives
#### Associations entre variables environnementales (Q mode)
?dist # cette fonction génère aussi des matrices de distance 
env.de <- dist(env.z, method = "euclidean") # matrice de distances euclidiennes des données environnementales centrées-réduites
windows() 
coldiss(env.de, diag=TRUE) 

### Corrélations entre les variables environnementales (R mode)
(env.pearson<-cor(env)) # corrélation linéaire du r de Pearson 
round(env.pearson, 2)
(env.ken<-cor(env, method="kendall")) # corrélation de rang (tau de Kendall)
round(env.ken, 2)

### Types mixtes de variables environnementales
#### Association entre variables environnementales de types variés (mode Q)
##### Création d'un tableau de données fictif
var.g1 <- rnorm(30, 0, 1)
var.g2 <- runif(30, 0, 5)
var.g3 <- gl(3, 10)
var.g4 <- gl(2, 5, 30)
(dat2 <- data.frame(var.g1, var.g2, var.g3, var.g4))
str(dat2)
summary(dat2)
##### Générer l'indice de dissimilarité (distances) de Gower
?daisy
(dat2.dg <- daisy(dat2, metric="gower"))
coldiss(dat2.dg)

# Défi 1 Avancé ####

# Calculer à la main les distances de Bray Curtis et de Gower pour les abondances d'espèces CHA, TRU et VAI pour les sites 1, 2 et 3
# (réponse ci-dessous)
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# Distance de Bray-Curtis : d[jk] = (sum abs(x[ij]-x[ik]))/(sum (x[ij]+x[ik]))
(spe.challenge<-spe[1:3,1:3])
# Calcul du total d'abondance d'espèces par site
(Abund.s1<-sum(spe.challenge[1,]))
(Abund.s2<-sum(spe.challenge[2,]))
(Abund.s3<-sum(spe.challenge[3,]))
# Calculer les différences d'abondances de chaque espèce pour chaque paire de sites
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])) }
# Calcul de la distance de Bray-Curtis 
(db.s1s2<-Spec.s1s2/(Abund.s1+Abund.s2))
(db.s1s3<-Spec.s1s3/(Abund.s1+Abund.s3)) 
(db.s2s3<-Spec.s2s3/(Abund.s2+Abund.s3)) 
# Comparez vos résultats
(spe.db.challenge<-vegdist(spe.challenge, method="bray"))

# Distance de Gower :    d[jk] = (1/M) sum(abs(x[ij]-x[ik])/(max(x[i])-min(x[i])))
# Calculer le nombre de colonnes dans la base de données
M<-ncol(spe.challenge)
# Calculer les différences d'abondances de chaque espèce pour chaque paire de sites
Spe1.s1s2<-abs(spe.challenge[1,1]-spe.challenge[2,1])
Spe2.s1s2<-abs(spe.challenge[1,2]-spe.challenge[2,2])
Spe3.s1s2<-abs(spe.challenge[1,3]-spe.challenge[2,3])
Spe1.s1s3<-abs(spe.challenge[1,1]-spe.challenge[3,1])
Spe2.s1s3<-abs(spe.challenge[1,2]-spe.challenge[3,2])
Spe3.s1s3<-abs(spe.challenge[1,3]-spe.challenge[3,3])
Spe1.s2s3<-abs(spe.challenge[2,1]-spe.challenge[3,1])
Spe2.s2s3<-abs(spe.challenge[2,2]-spe.challenge[3,2])
Spe3.s2s3<-abs(spe.challenge[2,3]-spe.challenge[3,3])
# Calculer l'étendue d'abondance de chaque espèce entre les sites
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  
(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)))
# Comparez vos résultats
(spe.db.challenge<-vegdist(spe.challenge, method="gower"))


## 2.2. Transformation des données de composition des communautés ####

?decostand # cette fonction fournit des options de standardisation et de transformation pour ce type de données
# Transformer les données d'abondance en présences-absences (1-0)
(spe.pa<-decostand(spe, method="pa"))

# Transformation de Hellinger 
(spe.hel<-decostand(spe, method="hellinger"))

# Transformation du chi-carré
(spe.chi<-decostand(spe, method="chi.square"))


# Défi 2 Avancé ####

# Calculez à la main les distances de Hellinger et du chi-carré sur les données 
# (answer below)
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
#
# Transformation de Hellinger
# Premièrement, calculer l'abondance totale 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 site
(scale.spe<-spe/site.totals)
# Calculer la racine carrée des abondances d'espèces réduites
(sqrt.scale.spe<-sqrt(scale.spe))
# Comparer les résultats
sqrt.scale.spe
spe.hel
sqrt.scale.spe-spe.hel # ou: sqrt.scale.spe/spe.hel

# Transformation de chi-carré
# Premièrement, calculer l'abondance totale des espèces par site
(site.totals<-apply(spe, 1, sum))
# Calculer la racine carrée des totaux d'abondances d'espèces
(sqrt.spe.totals<-sqrt(apply(spe, 2, sum)))
# Réduire les abondances d'espèces en les divisant par les totaux par sites et par espèces
scale.spe2<-spe
for (i in 1:nrow(spe)) {
  for (j in 1:ncol(spe)) {
    (scale.spe2[i,j]<-scale.spe2[i,j]/(site.totals[i]*sqrt.spe.totals[j]))   }}
#Ajuster les abondances en les multipliant par la racine carrée du total de la matrice des espèces                                         
(adjust.scale.spe2=scale.spe2*sqrt(sum(rowSums(spe))))
# Comparer les résultats
adjust.scale.spe2
spe.chi
adjust.scale.spe2-spe.chi # ou: adjust.scale.spe2/spe.chi

## 2.3 Groupement ####

### Calculer la distance de Hellinger
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
head(spe.hel)# données d’abondances transformées Hellingerhead(spe.dhel)# matrice de distances de Hellinger entre les sites

### Comparaison du groupement à liens simples et à liens complets

# Faire le groupement à liens simples
spe.dhel.single<-hclust(spe.dhel, method=”single”)
plot(spe.dhel.single)

# Faire le groupement à liens complet
spe.dhel.complete<-hclust(spe.dhel, method=”complete”)
plot(spe.dhel.complete)

### Méthode de Ward

# Faire le groupement de Ward 
spe.dhel.ward<-hclust(spe.dhel, method=”ward.D2”)
plot(spe.dhel.ward)

# Refaire le dendrogramme en utilisant la racine carrée des distances
spe.dhel.ward$height<-sqrt(spe.dhel.ward$height)
plot(spe.dhel.ward)
plot(spe.dhel.ward, hang=-1) # hang=-1 permet d’afficher les objets sur la même ligne


#***********************************************************************************#

# 3.Ordination sans contrainte ####

## 3.1. Analyse en composantes principales (ACP) ####

### ACP sur les abondances d'espèces transformées avec Hellinger
?rda
#### Lancer l'ACP
spe.h.pca <- rda(spe.hel)
# Extraire les résultats
summary(spe.h.pca)

#### Extraire des parties spécifiques des résultats
summary(spe.h.pca) # résultats globaux
summary(spe.h.pca, display=NULL) # seulement les valeurs propres et leur contribution à la variance expliquée
eigen(cov(spe.hel)) # autre moyen de calculer les valeurs propres
(spe.scores <- scores(spe.h.pca, display="species", choices=c(1,2))) # scores des espèces sur les deux premiers axes de l'ACP
(site.scores <- scores(spe.h.pca, display="sites", choices=c(1,2))) # scores des sites sur les deux premiers axes de l'ACP
#Note: si vous ne spécifiez pas le nombre de composantes principales à extraire (ex. choices=c(1,2) ou choices=c(1:2)), alors les scores de toutes les composantes principales seront extraits.

# Identifier les axes significatifs avec le critère de Kaiser-Guttman
ev <- spe.h.pca$CA$eig
ev[ev>mean(ev)]
n <- length(ev)
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")

### ACP sur les variables environnementales
#### Lancer l'ACP
env.pca <- rda(env.z) # ou rda(env, scale=TRUE)

#### Extraire les résultats
summary(env.pca)
summary(env.pca, scaling=2)

# Identifier les axes significatifs avec le critère de Kaiser-Guttman 
ev <- env.pca$CA$eig
ev[ev>mean(ev)]
n <- length(ev)
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")


### Créer un biplot
#### Biplot de l'ACP sur les abondances d'espèces transformées (cadrage 1)
windows()
plot(spe.h.pca)
windows()
?biplot 
biplot(spe.h.pca)
windows()
plot(spe.h.pca, scaling=1, type="none", # cadrage 1 = biplot de distance : 
     # la distance entre les objets = approximation de leur distance euclidienne
     # les angles entre les descripteurs ne reflètent PAS leur corrélation
     xlab=c("PC1 (%)", round((spe.h.pca$CA$eig[1]/sum(spe.h.pca$CA$eig))*100,2)),
     ylab=c("PC2 (%)", round((spe.h.pca$CA$eig[2]/sum(spe.h.pca$CA$eig))*100,2)))
points(scores(spe.h.pca, display="sites", choices=c(1,2), scaling=1),
       pch=21, col="black", bg="steelblue", cex=1.2)
text(scores(spe.h.pca, display="species", choices=c(1), scaling=1),
     scores(spe.h.pca, display="species", choices=c(2), scaling=1),
     labels=rownames(scores(spe.h.pca, display="species", scaling=1)),
     col="red", cex=0.8)       

#### Biplot de l'ACP sur les variables environnementales (cadrage 2)
windows()
plot(env.pca)
windows()
plot(env.pca, scaling=2, type="none", # cadrage 2 = biplot de corrélations : 
     # la distance entre les objets N'EST PAS une approximation de leur distance euclidienne
     # mais les angles entre les vecteurs des descripteurs reflètent leur corrélation
     xlab=c("PC1 (%)", round((env.pca$CA$eig[1]/sum(env.pca$CA$eig))*100,2)),
     ylab=c("PC2 (%)", round((env.pca$CA$eig[2]/sum(env.pca$CA$eig))*100,2)),
     xlim=c(-1,1), ylim=c(-1,1))
points(scores(env.pca, display="sites", choices=c(1,2), scaling=2), 
       pch=21, col="black", bg="darkgreen", cex=1.2) 
text(scores(env.pca, display="species", choices=c(1), scaling=2),
     scores(env.pca, display="species", choices=c(2), scaling=2),
     labels=rownames(scores(env.pca, display="species", scaling=2)),
     col="red", cex=0.8)

# Défi 3 ####
# Lancer une ACP sur les données d'abondances d'acariens
# Quels sont les axes significatifs ?
# Quels groupes pouvez-vous identifier ?
# Quelles espèces sont sont liées à chaque groupe de site
data(mite)
mite.spe <- mite # les données proviennent de la librairie vegan
# (réponse ci-dessous)
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
#
# 
mite.spe.hel<-decostand(mite.spe, "hel")
mite.spe.h.pca<-rda(mite.spe.hel)        
ev<-mite.spe.h.pca$CA$eig
ev[ev>mean(ev)]
n=length(ev)
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")       
summary(mite.spe.h.pca, display=NULL)
windows()
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)  

## 3.2. Analyse en coordonnées principales (PCoA)
### Lancer une PCoA sur les abondances d'espèces transformées
?cmdscale
cmdscale(dist(spe.hel), k=(nrow(spe)-1), eig=TRUE)
?pcoa
spe.h.pcoa<-pcoa(dist(spe.hel))
### Extraire les résultats
spe.h.pcoa
### Créer le biplot
windows()
biplot.pcoa(spe.h.pcoa, spe.hel)

### Lancer une PCoA sur la matrice de distances Bray-Curtis
spe.db
spe.bray.pcoa<-pcoa(spe.db, correction="cailliez")
spe.bray.pcoa
windows()
biplot.pcoa(spe.bray.pcoa)
# La mesure de distance choisie influence grandement les résultats de la PCoA

# Défi 4 ####
# Lancer une PCoA sur les données d'acariens transformées Hellinger
# Quels sont les axes significatifs ?
# Quels groupes pouvez-vous identifier ?
# Quelles espèces sont sont liées à chaque groupe de site
# Comparez les résultats avec ceux de l'ACP
# (réponse ci-dessous)
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
#
# 
mite.spe.h.pcoa<-pcoa(dist(mite.spe.hel))
mite.spe.h.pcoa
windows()
biplot.pcoa(mite.spe.h.pcoa, mite.spe.hel)

#**********************************************************************************#
# 4. Analyses canoniques
## 4.1. Analyse canonique de redondances (ACR)

### Préparer les données
env <- subset(env, select = -das) # enlever la variable "distance from the source"

### Lancer l'ACR avec toutes les variables explicatives sur les abondances d'espèces
?rda
spe.rda <- rda(spe.hel~., data=env)

### Extraire les résultats
summary(spe.rda, display=NULL)

### Sélectionner les variables explicatives significatives par sélection progressive
?ordiR2step
ordiR2step(rda(spe.hel~1, data=env), scope= formula(spe.rda), direction= "forward", R2scope=TRUE, pstep=1000)
env.signif <- subset(env, select = c("alt", "oxy", "dbo"))

### Refaire l'ACR avec les variables significatives sur les abondances d'espèces
spe.rda.signif <- rda(spe.hel~., data=env.signif)

### Extraire les résultats
summary(spe.rda.signif, display=NULL)

### Calculer le R^2 ajusté de l'ACR
(R2adj <- RsquareAdj(spe.rda.signif)$adj.r.squared)

### Tester la signification du modèle et des axes
?anova.cca
anova(spe.rda.signif, step=1000)
anova(spe.rda.signif, step=1000, by="axis")

### Représenter les résultats graphiquement
#### Graphiques (cadrage 1 et cadrage 2) facile à faire
windows()
plot(spe.rda.signif, scaling=1, main="Triplot RDA (cadrage 1)")
windows()
plot(spe.rda.signif, scaling=2, main="Triplot RDA (cadrage 2)")

#### Graphique avancé - cadrage 1
windows()
plot(spe.rda.signif, scaling=1, main="Triplot RDA - cadrage 1", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
points(scores(spe.rda.signif, display="sites", choices=c(1,2), scaling=1),
       pch=21, col="black", bg="steelblue", cex=1.2)
text(scores(spe.rda.signif, display="species", choices=c(1), scaling=1),
     scores(spe.rda.signif, display="species", choices=c(2), scaling=1),
     labels=rownames(scores(spe.rda.signif, display="species", scaling=1)),
     col="grey", cex=0.8)  
arrows(0,0,
       scores(spe.rda.signif, display="bp", choices=c(1), scaling=1),
       scores(spe.rda.signif, display="bp", choices=c(2), scaling=1),
       col="red")
text(scores(spe.rda.signif, display="bp", choices=c(1), scaling=1)+0.05,
     scores(spe.rda.signif, display="bp", choices=c(2), scaling=1)+0.05,
     labels=rownames(scores(spe.rda.signif, display="bp", choices=c(2), scaling=1)),
     col="red", cex=1) 

#### Graphique avancé - cadrage 2
windows()
plot(spe.rda.signif, scaling=2, main="Triplot RDA - cadrage 2", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
points(scores(spe.rda.signif, display="sites", choices=c(1,2), scaling=2),
       pch=21, col="black", bg="steelblue", cex=1.2)
text(scores(spe.rda.signif, display="species", choices=c(1), scaling=2),
     scores(spe.rda.signif, display="species", choices=c(2), scaling=2),
     labels=rownames(scores(spe.rda.signif, display="species", scaling=2)),
     col="grey", cex=0.8)  
arrows(0,0,
       scores(spe.rda.signif, display="bp", choices=c(1), scaling=2),
       scores(spe.rda.signif, display="bp", choices=c(2), scaling=2),
       col="red")
text(scores(spe.rda.signif, display="bp", choices=c(1), scaling=2)+0.05,
     scores(spe.rda.signif, display="bp", choices=c(2), scaling=2)+0.05,
     labels=rownames(scores(spe.rda.signif, display="bp", choices=c(2), scaling=2)),
     col="red", cex=1)  

# Défi 5 ####
# Lancer une ACR des données environnementales des acariens sur leurs abondances
data(mite.env) # données de la librairie vegan
# Quelles sont les variables explicatives significatives ?  
# Quelle proportion de variance est expliquée par les variables significatives?
# Quels sont les axes significatifs de l'ACR ?
# Quels groupes de sites pouvez-vous identifier ? 
# Quelles espèces sont liées à chaque groupes de sites ?
# (réponses ci-dessous)
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
#
# 

mite.spe.rda <- rda(mite.spe.hel~., data=mite.env)
ordiR2step(rda(mite.spe.hel~1, data=mite.env), 
           scope= formula(mite.spe.rda), direction= "forward", R2scope=TRUE, pstep=1000)
mite.env.signif <- subset(mite.env, 
                          select = c("WatrCont", "Shrub", "Substrate", "Topo", "SubsDens"))
mite.spe.rda.signif <- rda(mite.spe~., data=mite.env.signif) 
summary(mite.spe.rda.signif, display=NULL)
(R2adj <- RsquareAdj(mite.spe.rda.signif)$adj.r.squared)
anova(mite.spe.rda.signif, step=1000)
anova(mite.spe.rda.signif, step=1000, by="axis")          
windows()
plot(mite.spe.rda.signif, scaling=1, main="Triplot RDA - scaling 1", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
points(scores(mite.spe.rda.signif, display="sites", choices=c(1,2), scaling=1),
       pch=21, col="black", bg="steelblue", cex=1.2)
text(scores(mite.spe.rda.signif, display="species", choices=c(1), scaling=1),
     scores(mite.spe.rda.signif, display="species", choices=c(2), scaling=1),
     labels=rownames(scores(mite.spe.rda.signif, display="species", scaling=1)),
     col="grey", cex=0.8)  
arrows(0,0,
       scores(mite.spe.rda.signif, display="bp", choices=c(1), scaling=1),
       scores(mite.spe.rda.signif, display="bp", choices=c(2), scaling=1),
       col="red")
text(scores(mite.spe.rda.signif, display="bp", choices=c(1), scaling=1)+0.05,
     scores(mite.spe.rda.signif, display="bp", choices=c(2), scaling=1)+0.05,
     labels=rownames(scores(mite.spe.rda.signif, display="bp", choices=c(2), scaling=1)),
     col="red", cex=1) 

## 4.2  ACR partielle : effets de la chimie de l'eau sur les abondances d'espèces en contrôlant l'effet des variables liées à la topographie
### Diviser la base de données sur l'environnement en deux tableaux
envtopo <- env[, c(1:3)] # Topographie : matrice explicative 1
names(envtopo)
envchem <- env[, c(4:10)] # Chimie de l'eau : matrice explicative 2
names(envchem)
### Lancer l'ACR partielle
spechem.physio <- rda(spe.hel, envchem, envtopo)
# or : 
spechem.physio2 <- rda(spe.hel ~ pH + dur + pho + nit + amm + oxy + dbo +
                         Condition(alt + pen + deb), data=env)
### Extraire les résultats
summary(spechem.physio, display=NULL)
### Calculer le R^2 ajusté de l'ACR partielle
(R2adj <- RsquareAdj(spechem.physio)$adj.r.squared)
### Test de la signification de l'ACR partielle et des axes
anova(spechem.physio, step=1000)
anova(spechem.physio2, step=1000, by="axis")
### Créer les triplots
#### Cadrage 1
windows(title="Partial RDA scaling 1")
plot(spechem.physio, scaling=1, main="Triplot ACR partielle - cadrage 1", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
points(scores(spechem.physio, display="sites", choices=c(1,2), scaling=1),
       pch=21, col="black", bg="steelblue", cex=1.2)
text(scores(spechem.physio, display="species", choices=c(1), scaling=1),
     scores(spechem.physio, display="species", choices=c(2), scaling=1),
     labels=rownames(scores(spechem.physio, display="species", scaling=1)),
     col="grey", cex=0.8)  
arrows(0,0,
       scores(spechem.physio, display="bp", choices=c(1), scaling=1),
       scores(spechem.physio, display="bp", choices=c(2), scaling=1),
       col="red")
text(scores(spechem.physio, display="bp", choices=c(1), scaling=1)+0.05,
     scores(spechem.physio, display="bp", choices=c(2), scaling=1)+0.05,
     labels=rownames(scores(spechem.physio, display="bp", choices=c(2), scaling=1)),
     col="red", cex=1) 
#### Cadrage 2
windows(title="Partial RDA scaling 2")
plot(spechem.physio, scaling=2, main="Triplot ACR partielle - cadrage 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)  

# Défi 6 ####
# Lancer l'ACR partielle des variables environnementales des acariens sur leur abondance
# en controlant pour l'effet des variables du substrat
# Est-ce que le modèle est sifnigicatif ?
# Quels sont les axes significatifs ?
# Interpréter le triplot.
# (réponse ci-dessous)
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
#
# 
mite.spe.subs<-rda(mite.spe.hel ~ Shrub + Topo
                  + Condition(SubsDens + WatrCont + Substrate), data=mite.env)
summary(mite.spe.subs, display=NULL)
(R2adj <- RsquareAdj(mite.spe.subs)$adj.r.squared)
anova(mite.spe.subs, step=1000)
anova(mite.spe.subs, step=1000, by="axis")
#### Triplot cadrage 1
windows(title="ACR partielle - cadrage 1")
plot(mite.spe.subs, scaling=1, main="ACR partielle - cadrage 1", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
points(scores(mite.spe.subs, display="sites", choices=c(1,2), scaling=1),
       pch=21, col="black", bg="steelblue", cex=1.2)
text(scores(mite.spe.subs, display="species", choices=c(1), scaling=1),
     scores(mite.spe.subs, display="species", choices=c(2), scaling=1),
     labels=rownames(scores(mite.spe.subs, display="species", scaling=1)),
     col="grey", cex=0.8)  
arrows(0,0,
       scores(mite.spe.subs, display="bp", choices=c(1), scaling=1),
       scores(mite.spe.subs, display="bp", choices=c(2), scaling=1),
       col="red")
text(scores(mite.spe.subs, display="bp", choices=c(1), scaling=1)+0.05,
     scores(mite.spe.subs, display="bp", choices=c(2), scaling=1)+0.05,
     labels=rownames(scores(mite.spe.subs, display="bp", choices=c(2), scaling=1)),
     col="red", cex=1) 
#### Triplot cadrage 2
windows(title="ACR partielle - cadrage 2")
plot(mite.spe.subs, scaling=2, main="ACR partielle - cadrage 2 2", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
points(scores(mite.spe.subs, display="sites", choices=c(1,2), scaling=2),
       pch=21, col="black", bg="steelblue", cex=1.2)
text(scores(mite.spe.subs, display="species", choices=c(1), scaling=2),
     scores(mite.spe.subs, display="species", choices=c(2), scaling=2),
     labels=rownames(scores(mite.spe.subs, display="species", scaling=2)),
     col="grey", cex=0.8)  
arrows(0,0,
       scores(mite.spe.subs, display="bp", choices=c(1), scaling=2),
       scores(mite.spe.subs, display="bp", choices=c(2), scaling=2),
       col="red")
text(scores(mite.spe.subs, display="bp", choices=c(1), scaling=2)+0.05,
     scores(mite.spe.subs, display="bp", choices=c(2), scaling=2)+0.05,
     labels=rownames(scores(mite.spe.subs, display="bp", choices=c(2), scaling=2)),
     col="red", cex=1)  

# 4.3. Partitionnement de la variation avec l'ACR partielle
?varpart
vegandocs("partitioning.pdf")
## Partitionnement de la variation avec toutes les variables explicatives
spe.part.all <- varpart(spe.hel, envchem, envtopo)
spe.part.all
windows(title="Partitionnement de la variation - toutes les variables")
plot(spe.part.all, digits=2)  
## Partitionnement de la variation avec les variables explicatives significatives
### Sélection progressive des variables dans chaque sous-ensemble de variables
spe.chem <- rda(spe.hel~., data=envchem)
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 )])
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)]
### Partitionnement de la variation
(spe.part <- varpart(spe.hel, envchem.pars, envtopo.pars))
windows(title="Partitionnement de la variation - sous-ensembles parcimonieux")
plot(spe.part, digits=2)
# Test de toutes les fractions pouvant être testées
anova(rda(spe.hel, envchem.pars), step=1000) # test des fractions [a+b]
anova(rda(spe.hel, envtopo.pars), step=1000) # test des fractions [b+c]
env.pars <- cbind(envchem.pars, envtopo.pars) 
anova(rda(spe.hel, env.pars), step=1000) # test des fractions [a+b+c]
anova(rda(spe.hel, envchem.pars, envtopo.pars), step=1000) # test des fractions [a]
anova(rda(spe.hel, envtopo.pars, envchem.pars), step=1000) # test des fractions [c]

# Défi 7 ####
# Faire le partitionnement de la variation pour les abondances d'acariens
# avec une matrice incluant les variables significatives du substrat (SubsDens, WaterCont and Substrate)
# puis avec une matrice incluant les autres variables significatives (Shrud and Topo)
# Quelle proportion de variance est expliquée par chaque matrice?
# Quelles sont les fractions significatives ?
# (réponse ci-dessous)
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
#
# 
str(mite.env)
(mite.subs<-mite.env[,c(1,2,3)])
(mite.other<-mite.env[,c(4,5)])
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.mite.other <- rda(mite.spe.hel~., data=mite.other)
R2a.all.chem <- 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)])
### Partionnement de la variation
(mite.spe.part <- varpart(mite.spe.hel, ~WatrCont+Substrate, ~Shrub+Topo,
                          data=mite.env))
windows(title="Partitionnement de la variation - sous-ensembles parcimonieux")
plot(mite.spe.part, digits=2)
# Test de toutes les fractions pouvant être testées
anova(rda(mite.spe.hel~ WatrCont+Substrate, data=mite.env), step=1000) #  test des fractions [a+b]
anova(rda(mite.spe.hel~Shrub+Topo, data=mite.env), step=1000) # test des fractions [b+c]
(env.pars <- cbind(mite.env[,c(2,3,4,5)])) 
anova(rda(mite.spe.hel~ WatrCont+Substrate+Shrub+Topo, data=env.pars), step=1000) # test des fractions [a+b+c]
anova(rda(mite.spe.hel~WatrCont+Substrate + Condition(Shrub+Topo), data=env.pars), step=1000) # test des fractions [a]
anova(rda(mite.spe.hel~Shrub+Topo+ Condition(WatrCont+Substrate ), data=env.pars), step=1000) # test des fractions [c]



## 4.4. Courbes de réponses (PRC)
?prc
# van den Brink, P.J. & ter Braak, C.J.F. (1999). 
# Courbes de réponses principales: Analyse de la réponse temporelle de communautés biologiques au stress
# Environmental Toxicology and Chemistry, 18, 138-148
# Expérience sur Chlorpyrifos et design expérimental
data(pyrifos)
week <- gl(11, 12, labels=c(-4, -1, 0.1, 1, 2, 4, 8, 12, 15, 19, 24))
dose <- factor(rep(c(0.1, 0, 0, 0.9, 0, 44, 6, 0.1, 44, 0.9, 0, 6), 11))
# PRC
mod <- prc(pyrifos, dose, week)
mod            # ACR
summary(mod)   # PRC
logabu <- colSums(pyrifos)
windows()
plot(mod, select = logabu > 100)
# Les permutations devraient s'effectuer sur une semaine 
# et nous ne sommes intéressés que par le permier axe
anova(mod, strata = week, first=TRUE, perm.max = 100)

# Les PRC peuvent être utiles en écologie des communautés
# pour comparer l'évolution temporelle d'une communauté dans des sites expérimentaux (ex. écosystèmes restaurés)
# à l'évolution temporelle d'une communauté dans des sites contrôles (ex. écosystèmes de référence)

# Défi 8 ####
# Faire une PRC avec les données de poissons de la rivière Doubs en enlevant le dernier site, considérant que:
# les 14 premiers sites correspondent à 4 écosystèmes de référence
# les 14 derniers sites correspondent à 4 écosystèmes restaurés
# les premières sept lignes correspondent à l'échantillonnage successif dans les sites de références 1 durant les années 1,2,3,4,5,6 et 7
# les lignes 8 à 14 correpondent à l'échantillonnage successif dans les sites de références 2 durant les années 1,2,3,4,5,6 et 7
# etc.
# Est-ce que la restauration permet le rétablissement des communautés de référence?  
# Quelles espèces sont liées aux écosystèmes de référence, et aux écosystèmes restaurés ?
# (réponse ci-dessous)
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
#
# 
dim(spe.hel)
spe.hel<-spe.hel[-29,]
reference<-rep(c(0), 14)
restored<-rep(c(1), 14)
type.site<-as.factor(c(reference, restored))
year<-as.factor(rep(c(1,2,3,4,5,6,7), 4))
spe.hel.prc<-prc(spe.hel, type.site, year)
spe.hel.prc
summary(spe.hel.prc)
plot(spe.hel.prc)

#***************************************************************************************#

# 5. Autres méthodes d'ordination utiles

?cca # Analyse canonique des correspondances (CCA)

?metaMDS # Cadrage multidimensionnel non-métrique (nMDS)

help(lda, package=MASS) # Analyse linéaire discriminante

?CCorA # Analyse canonique de corrélations

help(coinertia, package=ade4) # Analyse de co-inertie

help(mfa, package=ade4) # Analyse factorielle multiple 

https://r-forge.r-project.org/R/?group_id=195  # les librairies AEM et PCNM contiennes des fonctions pour effectuer des analyses spatiales