equipe_5_oral_3.Rmd

---
title: "Les vaches,<br>les amphibiens<br>et vous"
author: "Équipe 5 - Florent Déry, Morgane Le Goff, Myriam Trottier-Paquet"
date: "`r Sys.Date()`"
output:
  ioslides_presentation:
    widescreen: true
    smaller: true
editor_options: 
  chunk_output_type: console
---

<style>
pre {
    line-height: 1.2em;
    font-size: 10px;
}
</style>
## Contexte, objectif et question de recherche (étape 1)

Les pâturages occupent actuellement 26% des surfaces terrestres libres de glaces et cette proportion est appelée à augmenter. Il est essentiel de comprendre l'impact du bétail sur la faune et la flore sauvages afin d’atténuer la perte de biodiversité mondiale due à l’augmentation de la superficie de terres converties en pâturages pour le bétail.

*Objectif :*
La présente étude cherche à quantifier les impacts du pâturage du bétail sur l'abondance et la diversité des amphibiens dans les zones agricoles du centre-sud de la Floride.

*Questions :*
1. Est-ce que la communauté d’amphibiens diffère entre les étangs dont le bétail est exclu de ceux non-exclu?
2. Est-ce que l’intensité de l’utilisation par le bétail affecte la communauté d’amphibiens en prenant en compte les variables environnementales et la communauté animale?
3. Quelle sont les variables les plus influentes sur la communauté d’amphibiens et combien de variation ces variables expliquent-elles?

*Objectif :* 
Vérifier l'impact du pâturage par le bétail (vaches) sur les communautés d'amphibiens d'étangs du Centre-Sud de la Floride, É.-U. 

### Hypothèses et prédictions (étape 2) 
La présence et l'intensité d'utilisation du bétail modifient les conditions physico-chimiques des étangs, et ainsi l'abondance des invertébrés, poissons et écrevisses, ce qui détermine la composition des communautés amphibiennes. 
Prédictions: 

* Les conditions physico-chimiques sont différentes entre les étangs cloturés vs non-cloturés. 
    + Les vaches seront corrélées positivement aux nitrites
    + Les vaches seront corrélées positivement à la *temperature*?
    + Les vaches seront corrélées positivement à la conductivité
    + Les vaches seront corrélées négativement au pH
    + Les vaches seront corrélées négativement à l'O² 
* Les espèces généralistes seront présentes dans tous les sites car elles sont tolérantes à une grande étendue de conditions.  
    + À l'inverse, les espèces fragiles ou spécialistes occuperont des sites avec des conditions optimales
  
## Variables du système (étape 3)
Le jeu de données contient cinq *composantes*:

* les variables reliées au *pâturage*
* les caractéristiques physico-chimiques
* les facteurs biotiques 
* les caractéristiques topographiques
* la matrice d'abondance des amphibiens

Jettons un coup d'oeil à la structure (`glimpse()`) puis au `summary()` du jeu de données:

## Variables du système (étape 3)
```{r load_n_glimpse_dat, echo=F, warnings=F, message=F}
require(tidyverse)
dat=read.delim("data/jeu_ecologique.txt") %>% 
  mutate(etang=as.factor(as.character(etang)),
         cloture=factor(cloture)) %>% 
  rename(vache_std=vache) %>% 
  arrange(aire) 
glimpse(dat)
```
## Variables du système (étape 3)
```{r summary_dat_xs, echo=F, warnings=F, message=F}
summary(dat[,1:13])
```
## Variables du système (étape 3)
```{r summary_dat_ys, echo=F, warnings=F, message=F}
summary(dat[,14:22])
```
## Rôle, type et distribution des variables (étape 4)
La variable dépendante (Y) est représentée par une matrice d'abondances d'amphibiens. Les données d'abondance ont été standardisées par unité d'effort et par abondance totale pour chaque espèce. 

* Les valeurs d'abondance relative sont entre 0 et 0.86. 
* Les zéros (absences) composent 70.5% de la matrice. 

```{r format_Ys, echo=F, warnings=F, message=F}
amphi<-dat %>% 
  select(etang,A_Gr:N_Vi) %>% 
  column_to_rownames("etang")
```

```{r step_4_quick_look_to_Ys, echo=F, warnings=F, message=F, eval=F}
range(amphi)# Étendue
sum(amphi==0) ; sum(amphi==0)/(nrow(amphi)*ncol(amphi))# Nombre de zéros, puis proportion de zéros dans la communauté
```

## Étape 4 - distribution des Ys
```{r distribution_ys,  echo=F, warnings=F, message=F, fig.align="center", fig.dim=c(10,5) }
amphi %>% 
  pivot_longer(1:9,names_to = "espece", values_to = "prop") %>% 
  mutate(lng_nm=case_when(espece=="A_Gr"~ "Acris gryllus",
                          espece=="G_Ca"~ "Gastrophryne carolinensis",
                          espece=="S_La"~ "Siren lacertina",
                          espece=="H_Ci"~ "Hyla cinerea",
                          espece=="H_Fe"~ "Hyla Femoralis",
                          espece=="H_Sq"~ "Hyla squirella",
                          espece=="L_Gr"~ "Lithobates grylio",
                          espece=="L_Sp"~ "Lithobates sphenocephalus",
                          espece=="N_Vi"~ "Notophthalmus viridescens")) %>% 
ggplot(aes(prop))+
  geom_histogram()+
  geom_density()+
  facet_wrap(~lng_nm, scales="free_y")+
  theme_bw()
```

## Étape 4 - Variables indépendantes X

* La variable `cloture` est une variable catégorique à deux facteurs (Yes, No). 
* La variable `vache_std` (Intensité d'utilisation de l'étang par le bétail) est une variable numérique continue dont la distribution [min ; max = 0 ; 0.6] est biaisée à gauche.
* 5 étangs sont clôturés ; 35 étangs ne l'étaient pas. Parmi ces derniers, 3 étangs ne semblent pas avoir été utilisés par les vaches (`vache_std`=0)

```{r distribution_xs,  echo=F, warnings=F, message=F, fig.align="center", fig.dim=c(10,4) }
p1= ggplot(dat, aes(cloture, fill=cloture))+
  geom_bar()+
  xlab("Étang clôturé ou non?")+
  ylab("Nombre d'étangs")+
  scale_fill_manual(values=c("steelblue4", "darkgoldenrod2"))+
  theme_bw()

p2=ggplot(dat, aes(vache_std))+  
geom_histogram(aes(fill=cloture))+
  geom_density()+
  xlab("Nombre de fèces standardisé")+
  ylab("Nombre d'étangs")+
  scale_fill_manual(values=c("steelblue4", "darkgoldenrod2"))+
  theme_bw()+
  theme(legend.position = c(0.8, 0.8))

require(patchwork)

p1 + p2 + plot_layout(ncol=2, widths=c(0.33, .67))
```

## Étape 4 - Variables | Caractéristiques physico-chimiques

* `conductivite`: numérique continue, distribution normale
* `nitrates` : numérique continue, distribution déviée à gauche
* `temperature` (Celcius): numérique continue, distribution déviée à gauche
    + Quelques valeurs extrèmes **-> on y reviendra!**

```{r distribution_covariables_chim_a,  echo=F, warnings=F, message=F, fig.align="center", fig.dim=c(8,4) }
par(mfrow=c(1,3))
hist(dat$conductivite,
     main="Conductivité",
     xlab="Conductivité",
     col="cornflowerblue")
hist(dat$nitrates,
     main="Nitrates",
     xlab="Nitrates",
     col="cornflowerblue")
hist(dat$temperature,
     main="Température",
     xlab="Température en degré C",
     col="cornflowerblue")
range(dat$pH)
```
## Étape 4 - Variables | Caractéristiques physico-chimiques

* Oxygène: numérique continue, distribution normale
    + Il y a des valeurs négatives!!!  **-> on y reviendra!**
* pH: numérique continue, distribution déviée à gauche
    + Des étangs sont plus acides que du jus de citron, et presqu'aussi basiques que de l'eau de Javel! **-> on y reviendra!**
```{r distribution_covariables_chim_b,  echo=F, warnings=F, message=F, fig.align="center", fig.dim=c(10,4.3) }
par(mfrow=c(1,2))
hist(dat$oxygene,
     main="Oxygène",
     xlab="Oxygène dissous",
     col="cornflowerblue")
hist(dat$pH,
     main="pH",
     xlab="pH",
     col="cornflowerblue")
```

## Étape 4 - Covariables | Caractéristiques topographiques

* Profondeur (cm): numérique continue, distribution normale
* Aire (m²): numérique continue, distribution *relativement* normale

```{r distribution_covariables_topo,  echo=F, warnings=F, message=F, fig.align="center", fig.dim=c(10,4.3) }
par(mfrow=c(1,2))
hist(dat$profondeur,
     main="Profondeur",
     xlab="Profondeur en cm",
     col="cornflowerblue")
hist(dat$aire,
     main="Aire",
     xlab="Aire en m^2",
     col="cornflowerblue")
```

## Étape 4 - Variables | Facteurs biotiques

* Invertébrés: numérique continue, distribution déviée à gauche
* Écrevisses: numérique continue, distribution déviée à gauche avec beaucoup de zéros **-> on y reviendra!**
* Poissons: numérique continue, distribution avec beaucoup de zéros **-> on y reviendra!**

```{r distribution_covariables_biotiques,  echo=F, warnings=F, message=F, fig.align="center", fig.dim=c(8,4) }
par(mfrow=c(1,3))
hist(dat$invertebres,
     main="Invertébrés",
     xlab="Abondance d'invertébrés prédateurs",
     ylab="Nombre d'étangs",
     col="cornflowerblue")
hist(dat$ecrevisses,
     main="Écrevisses",
     xlab="Abondance d'écrevisses",
     ylab="Nombre d'étangs",
     col="cornflowerblue")
hist(dat$poissons,
     main="Poissons",
     xlab="Abondance de poissons",
     ylab="Nombre d'étangs",
     col="cornflowerblue")

```

## Interactions potentielles à tester (étape 5)

Aucune. 

## Exploration des données (étape 6) - Données aberrantes en X (6.1) 

* Les valeurs de vache_std sont entre 0 et 0.575
* Il ne semble pas y avoir de données aberrantes, mais 3 étangs non-clôturés ne semblent pas avoir été visités par les vaches. 

```{r abberent_Xs_vache, echo=F, warnings=F, message=F, fig.align="center", fig.dim=c(8,4) }
range(dat$vache_std) 
p3=ggplot(dat)+
  geom_histogram(aes(vache_std))+
  xlab("")+
  ylab("Fréquence")+
    theme_bw()
p4=ggplot(dat)+
  geom_point(aes(y=rownames(dat), x=vache_std))+
  xlab("")+
  ylab("Index (No de ligne)")+
    theme_bw()
p5=ggplot(dat)+
  geom_boxplot(aes(x=vache_std))+
  xlab("")+
  ylab("")+
  theme_bw()+
  theme(axis.line.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())
require(cowplot)
p6=plot_grid(p3,p5, nrow=2, rel_heights = c(3.5, 1))
plot_grid(p6,p4, ncol=2)+
  labs(caption="Nombre de fèces standardisé")
```

## 6.1 - Données aberrantes en X 

* O²: Il y a des valeurs négatives!!!  
    + pour la suite de nos analyses, nous avons décidé de retirer les sites ayant des valeurs d'oxygène inférieur à 0. Ces valeurs sont 
- pH: Des étangs sont plus acides que du jus de citron, et presqu'aussi basiques que de l'eau de Javel! 

```{r abberent_Xs_physico, echo=F, warnings=F, message=F, fig.align="center", fig.dim=c(8,4) }
Chimie<-dat[ , c(5:9)]
l1 <- function() {
  par(
    mar=c(1,1,1,1),
    mfrow=c(1,5)
  )
for(i in 1:ncol(Chimie)) {
  hist(Chimie[,i], main=names(Chimie)[i], breaks=30)
}
}

l2 <- function() {
  par(
    mar=c(1,1,1,1),
    mfrow=c(1,5)
  )
  for(i in 1:ncol(Chimie)) {
    boxplot(Chimie[,i],horizontal = T, main=NULL)
  }
}
library(gridGraphics); library(cowplot)
plot_grid(ggdraw(l1) , ggdraw(l2), ncol=1, rel_heights = c(.8, .2))

```
## 6.1 - Données aberrantes en X 

* O²
    + Pour la suite de nos analyses, nous avons décidé de retirer les sites (n=2) ayant des valeurs d'oxygène inférieures à 0. Ces valeurs sont vraissemblablement des erreurs de mesures, on ne peut pas assumer qu'elles étaient à 0, ou à une autre valeur. 
* pH
    + Dans la section "Matériel supplémentaire", nous avons vérifié les valeurs de pH supérieures à 9 et inférieures à 4. 
      - Rien ne semblait anormal au niveau des autres paramètres mesurés à ces étangs, sauf peut-être la présence de poissons mêmne lorsque le pH était de 1.88. Ceci nous pousse à croire que la variable pH n'est probablement pas fiable. 
      - *Pour la suite, nous allons ignorer les variables de pH pour les analyses*

## 6.1 - Données aberrantes en X 

Pas de données particulièrement aberrantes en ce qui concerne la topographie des étangs

```{r aberrant_Xs_topo, echo=F, warnings=F, message=F, fig.align="center", fig.dim=c(8,4) }
Top<-dat[ , c(4, 10)]

l1=function()  {
  par(
  mfrow=c(1,2),
    mar=c(1,1,1,1)
  )
for(i in 1:ncol(Top)) {
  hist(Top[,i], main=names(Top)[i], breaks=30)
}
}
l2=function()  {
  par(
  mfrow=c(1,2),
    mar=c(1,1,1,1)
  )
for(i in 1:ncol(Top)) {
  boxplot(Top[,i], main=NULL, horizontal=T)
}
}

plot_grid(ggdraw(l1) , ggdraw(l2), ncol=1, rel_heights = c(.8, .2))
```

## 6.1 - Données aberrantes en X 

* En ce qui concerne les poissons et les écrevisses, il semble y avoir beaucoup de zéros. Il faudra peut-être transformer les variables poissons et écrevisses. comme on a tenté de le faire, sans grand succès, à la section "Matériel supplémentaire". 
* Les écrevisses étant des prédateurs des grenouilles (tétards), on pourrait peut-être travailler en absence-présence d'écrevisses dans nos analyses, si cette variable est problématique. 
* On a vérifié l'influence de travailler avec des variables binaires pour les comptes de poissons et d'écrevisses à la section "Matériel supplémentaire". 

*Pour la suite, nous travaillerons avec les données originales*, qu'on remet graphiquement ici:
```{r aberrant_Xs_biotiques, echo=F, warnings=F, message=F, fig.align="center", fig.dim=c(8,3) }
living<-dat[ , c("invertebres", "poissons", "ecrevisses")]

l1=function () {
par(
  mfrow=c(1,3),
    mar=c(1,1,1,1)
  )
for(i in 1:ncol(living)) {
  hist(living[,i], main=names(living)[i], breaks=30)
} 
}

l2=function () {
par(
  mfrow=c(1,3),
    mar=c(1,1,1,1)
  )
for(i in 1:ncol(living)) {
  boxplot(living[,i], main=NULL, horizontal=T)
}
}

plot_grid(ggdraw(l1) , ggdraw(l2), ncol=1, rel_heights = c(.8, .2))
```

## 6.1 - Données aberrantes en Y

* On a déjà jetté un coup d'oeil à la distribution des abondances relatives d'amphibien à l'étape 4. Le voici en rappel, mais seulement avec des boxplots. 
* Pas vraiment de données aberrantes (beaucoup de 0, mais pas aberrant en soi), si ce n'est que de quelques espèces rares. 

```{r aberrant_Ys, echo=F, warnings=F, message=F, fig.align="center", fig.dim=c(10,5) }
amphi<-dat %>% 
  select(etang,A_Gr:N_Vi) %>% 
  column_to_rownames("etang") 
amphi %>% 
  pivot_longer(1:9,names_to = "espece", values_to = "prop") %>% 
  mutate(lng_nm=case_when(espece=="A_Gr"~ "Acris gryllus",
                          espece=="G_Ca"~ "Gastrophryne carolinensis",
                          espece=="S_La"~ "Siren lacertina",
                          espece=="H_Ci"~ "Hyla cinerea",
                          espece=="H_Fe"~ "Hyla Femoralis",
                          espece=="H_Sq"~ "Hyla squirella",
                          espece=="L_Gr"~ "Lithobates grylio",
                          espece=="L_Sp"~ "Lithobates sphenocephalus",
                          espece=="N_Vi"~ "Notophthalmus viridescens")) %>% 
ggplot()+
  geom_boxplot(aes(prop), varwidth = F, width=.1)+
  facet_wrap(~lng_nm, scales="free_x")+
  theme_bw()+
  theme(axis.text.y = element_blank())

```

## 6.1 - Données aberrantes en Y

* On peut toutefois inspecter le compte d'espèces à chaque site (convertir les abondances relatives en 0-1 pour chaque site). On peut aussi jetter un coup d'oeil à la richesse spécifique, avec le fait d'être cloturé ou pas (pourquoi pas?!). 

* On peut voir quelles espèces sont rares dans notre jeu de données à gauche, et que les seuls sites mono-spcifiques sont aussi sans clôture, à droite. 

```{r amphi_barplot, echo=F, warnings=F, message=F, fig.align="center", fig.dim=c(8,3)}
amphi<-dat %>% 
  select(etang,A_Gr:N_Vi) %>% 
  column_to_rownames("etang") %>% 
  vegan::decostand(amphi, method = "pa") %>% 
  mutate(cloture=dat$cloture)

tmp=data.frame(table(rowSums(amphi[amphi$cloture=="No", 1:9]))) %>% 
  bind_rows(data.frame(table(rowSums(amphi[amphi$cloture=="Yes", 1:9])))) %>% 
  mutate(cloture=c(rep("No", 6), rep("Yes", 4)))
colnames(tmp) =  c("especes_presentes", "nombre_de_sites", "cloture")

plot_grid(
  ggplot()+
    geom_col(aes(y=sort(colSums(amphi[,1:9])),
                 x=factor(names(sort(colSums(amphi[,1:9]))), 
                          levels=names(sort(colSums(amphi[,1:9])))))
             )+ 
    xlab("Espèce")+
    ylab("Nombre de présences\n(tous les sites confondus)")+
    theme_classic(),

  ggplot(tmp)+
    geom_col(aes(especes_presentes,nombre_de_sites, fill=cloture))+
    xlab("Richesse spécifique à chaque site")+
  scale_fill_manual(values=c("steelblue4", "darkgoldenrod2"))+
    scale_y_continuous("Nombre de sites", 
                       breaks=seq(0, 12, 2), 
                       limits=c(0,12.1))+
    theme_classic()+
    theme(legend.position = c(0.8, .8))
)

```

## 6.2 - Homogénéité de la variance et 6.3 - Normalité des Ys

* Dans notre cas, on ne s'attend pas à ce que la variance des Ys soit homogène (on ne fait pas des modèles linéaires où les résidus devraient être homogènes, etc.); qui plus est, ce n'est pas une supposition à respecter dans le cadre d'une analyse multivariée...

* Dans notre cas, on ne s'attend pas à ce que la variance des Ys soit normale non plus. 

## 6.4 - Problèmes de 0 dans les Y

Comme on a vu plus tôt à l'étape 4 les zéros (absences) composent 70.5% de la matrice. Il se pourrait d'ailleurs que la faible présence de certaines espèces (H_Fe, N_Vi et A_GR) pose des problèmes pour l'analyse visée. 

## 6.5 - Collinéarité des X

```{r, echo=F, message=F, warnings=F, fig.width=10, fig.height=6}
dat %>% 
  select(vache_std:poissons) %>%
  relocate(aire, .after = profondeur) %>% 
  psych::pairs.panels()

```

## 6.5 - Collinéarité des X

En ce qui concerne les variables numériques, on a quelques corrélations, mais rien de révoltant, surtout si on choisi une analyse qui n'est invalidée par ces collinéarités. 

Pêut-être un peu plus *d'intérêt* pour les corrélations avec `vache_std`: 

* Corrélé positivement aux `nitrates` et à la `conductivite`, 
* Corrélé négativement au `pH`

## 6.5 - Collinéarité des X
Variable cloture vs. variables *physico-chimiques*: rien d'alarmant, considérant qu'on a seulement 5 étangs clôturés. Semble y avoir une tendance entre le fait d'être clôturé et les nitrates.  

```{r, echo=F, message=F, warnings=F, fig.width=7, fig.height=4.3, fig.align="center"}
dat %>% 
  select(cloture, conductivite:pH) %>% 
  pivot_longer(conductivite:pH, names_to = "parametre", values_to = "mesure") %>% 
  ggplot(aes(y= mesure, x=cloture )) +
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(height=0, width=.1)+
  facet_wrap(~parametre, scales = "free_y")+
  theme_bw()
```

## 6.5 - Collinéarité des X
Variable cloture vs. variables *topographiques* : rien d'alarmant, considérant qu'on a seulement 5 étangs clôturés.

```{r, echo=F, message=F, warnings=F, fig.width=7, fig.height=4.3, fig.align="center"}
dat %>% 
  select(cloture, profondeur, aire) %>% 
  pivot_longer(profondeur:aire, names_to = "caracteristique", values_to = "mesure") %>% 
  ggplot(aes(y= mesure, x=cloture )) +
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(height=0, width=.1)+
  facet_wrap(~caracteristique, scales = "free_y")+
  theme_bw()

```

## 6.5 - Collinéarité des X
Variable clôture vs. variables *biotiques*

```{r, echo=F, message=F, warnings=F, fig.width=7, fig.height=4.3, fig.align="center"}
dat %>% 
  select(cloture, invertebres:poissons) %>% 
  pivot_longer(invertebres:poissons, names_to = "taxon", values_to = "compte") %>% 
  ggplot(aes(y= compte, x=cloture )) +
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(height=0, width=.1)+
  facet_wrap(~taxon, scales = "free_y")+
  theme_bw()

```

## 6.6 - Liens X et Y
Abondances relatives en fonction des paramètres physico-chimiques

  > Ligne + interval = smooth ('loess'; 'y ~ x').
    
```{r echo=F, fig.align="center", fig.height=4.3, fig.width=7, message=FALSE, warnings=F}
set.seed(1)
pal=paletteer::paletteer_d("colorBlindness::paletteMartin") %>%  sample(9)

dat%>% 
  pivot_longer(A_Gr:N_Vi, names_to = "espece", values_to = "prop") %>% 
  pivot_longer(conductivite:pH, names_to = "parametre", values_to = "mesure") %>% 
  ggplot(aes(x= mesure, y= prop))+
  geom_point(aes(color=espece), size=2) +
    geom_smooth()+
   scale_color_manual(values=pal)+
 facet_wrap(~parametre, scales= "free_x")+
  ylab("Proportion par site")+
  xlab("Mesure du paramètre")+
  theme_bw()

 # %>% 
 #  mutate(prof_cat=factor(
 #    case_when(
 #      profondeur < 181 ~ "very_shallow", 
 #      profondeur %in% c(181:242)~ "shallow",
 #      profondeur %in% c(243:330)~ "medium",
 #      profondeur %in% c(331:460)~ "deep",
 #      profondeur > 460 ~ "very_deep"    
 #    ), levels= c("very_shallow", 
 #                 "shallow",
 #                 "medium",
 #                 "deep",
 #                 "very_deep")))
```

## 6.6 - Liens X et Y
Abondances relatives en fonction des caractéristiques topographiques 

> Ligne + interval = smooth ('loess'; 'y ~ x').

```{r, echo=F, message=F, warnings=F, fig.width=7, fig.height=4.3, fig.align="center"}
dat%>% 
  pivot_longer(A_Gr:N_Vi, names_to = "espece", values_to = "prop") %>% 
  select(prop, profondeur, aire, espece) %>% 
  pivot_longer(profondeur:aire, names_to = "caracteristique", values_to = "mesure") %>% 
  ggplot(aes(x= mesure, y= prop))+
  geom_point(aes(color=espece), size=2) +
    geom_smooth()+
   scale_color_manual(values=pal)+
 facet_wrap(~caracteristique, scales= "free_x")+
  ylab("Proportion par site")+
  xlab("Mesure du paramètre")+
  theme_bw()

```

## 6.6 - Liens X et Y
Abondances relatives en fonction des facteurs biotiques

> Ligne + interval = smooth ('loess'; 'y ~ x').

```{r, echo=F, message=F, warnings=F, fig.width=10, fig.height=4.3, fig.align="center"}
dat%>% 
  pivot_longer(A_Gr:N_Vi, names_to = "espece", values_to = "prop") %>% 
  pivot_longer(invertebres:poissons, names_to = "facteur_biotique", values_to = "compte") %>% 
  ggplot(aes(x= compte, y= prop))+
  geom_point(aes(color=espece), size=2) +
    geom_smooth()+
   scale_color_manual(values=pal)+
 facet_wrap(~facteur_biotique, scales= "free_x")+
  ylab("Proportion par site")+
  xlab("Mesure du paramètre")+
  theme_bw()

```

## 6.6 - Liens X et Y
Abondances relatives en fonction de l'intensité des vaches et de la présence de cloture.

> Ligne + interval = smooth ('loess'; 'y ~ x').

```{r, echo=F, message=F, warnings=F, fig.width=10, fig.height=4.3, fig.align="center"}
p1=dat%>% 
  pivot_longer(A_Gr:N_Vi, names_to = "espece", values_to = "prop") %>% 
  ggplot(aes(x= cloture, y= prop))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(aes(color=espece),height=0, width=.1)+
   scale_color_manual(values=pal)+
  ylab("Proportion par site")+
  xlab("Cloturé?")+
  theme_bw()+
  theme(legend.position = "none")

p2= dat%>% 
  pivot_longer(A_Gr:N_Vi, names_to = "espece", values_to = "prop") %>% 
  ggplot(aes(x= vache_std, y= prop))+
  geom_point(aes(color=espece), size=2) +
    geom_smooth()+
   scale_color_manual(values=pal)+
  ylab("Proportion par site")+
  xlab("Intensité d'utilisation par les vaches")+
  theme_bw()

require(patchwork)
p1+p2+plot_layout(ncol=2, guides = "collect")
```

## 6.7 - Indépendance des Y

Richesse spécifique à chaque site:
* Les données semblent indépendantes. 

Plus ou moins pour vérifier l'indépendance des sites, mles étangs avec plus d'espèces semblent être plus creux et plus gros. 

```{r, echo=F, message=F, warnings=F, fig.width=10, fig.height=4.3, fig.align="center"}
# Richesse spécifique vs sites
amphi<-dat %>% 
  select(etang,A_Gr:N_Vi) %>% 
  column_to_rownames("etang") %>% 
  vegan::decostand(amphi, method = "pa") 

sit.pres <- apply(amphi > 0, 1, sum)

sit.pres= data.frame(sit.pres) %>% 
  rownames_to_column(var = "etang") %>% 
  inner_join(., dat[, c("etang", "profondeur","aire")]) %>% 
  rename(richesse= sit.pres)

require(ggplot2)

plot_grid(
ggplot(sit.pres)+
  #geom_point( aes(x=etang, y=richesse))+
  geom_step(aes(x=etang, y=richesse, group=1))+
  xlab("") +
  scale_y_continuous("Richesse spécifique",
                     breaks= 1:6) +
    ggtitle("Ordre dans les données")+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90)),

sit.pres %>% 
  mutate(etang=factor(etang, levels=etang[order(.$profondeur)], ordered = T)) %>% 
ggplot()+
  #geom_point( aes(x=etang, y=richesse))+
  geom_step(aes(x=factor(etang), y=richesse, group=1))+
  xlab("Site") +
  scale_y_continuous("",
                     breaks= 1:6) +
  ggtitle("Peu profond ---> très creux")+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90)),

sit.pres %>% 
  mutate(etang=factor(etang, levels=etang[order(.$aire)], ordered = T)) %>% 
ggplot()+
  #geom_point( aes(x=etang, y=richesse))+
  geom_step(aes(x=factor(etang), y=richesse, group=1))+
  xlab("") +
  scale_y_continuous("",
                     breaks= 1:6) +
  ggtitle("Moins gros ---> plus gros")+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90)),
ncol=3)
```

## 6.8 - Collinéarité des Y

```{r, echo=F, message=F, warnings=F, fig.width=9, fig.height=5.5, fig.align="center"}
dat%>% 
  select(A_Gr:N_Vi) %>% 
  psych::pairs.panels()
```

## Faits saillants de l'exploration (étape 6):
* Les données d'abondance ont été standardisées par unité d'effort et par abondance totale pour chaque espèce (somme pour chaque espèce est égale à 1). 
    + Près de 71% des données de la matrice d'abondance sont égales à 0.
* Les variables physico-chimiques devront être standardisées
    + Certaines variables physico-chimiques semblent avoir des mesures erronnées.
    + L'oxygène et le pH, sont plus problématique. On va ignorer les valeurs de pH, et retirer les sites d'O² négatifs. 
    + On a vérifié les valeurs extrèmes de temperature (section matériel supplémentaire), rien à signaler. La transformation log ne colle pas particulièrement les valeurs extrèmes au reste des données.
* La transformation log des poissons ou des écrevisses n'atténue pas la sureprésentation des zéros vus dans la distribution des données brutes.  
* La topographie des étangs pourrait influencer d'autres facteurs physico-chimiques, et donc l'aire et la profondeur devraient être considérés comme des covariables. 
    + Les étangs plus creux et plus grands semblent avoir une richesse spécifique supérieure. 

## Étape 7. Choix de l'analyse

**Pourquoi une analyse multivariée?** 
* Analyse de données complexes en terme de dimensionalité: Plusieurs variables avec plusieurs dimensions de manière simultanée - On veut analyser la réponse de la communauté d’amphibiens (abondance de plusieurs espèces) à différents facteurs: conditions physico-chimiques, abondance d’autres espèces animales et topographie.

**Quelle analyse multivariée?**
* Nous avons décidé d'utiliser une analyse de redondance (RDA) partielle pour analyser comment la présence de clôture, l'utlisation des étangs par les vaches et les variables physico-chimiques en plus des facteurs biotiques font varier la communauté d'amphibiens, en prennant en compte l'effet de l'aire et de la profondeur des étangs comme co-variables.

**Pourquoi?**
* Approprié pour mettre en lien des variables liées aux espèces (matrice Y) et des variables environnementales (matrice X) comme dans notre cas. 
* Test non paramétrique donc pas besoin de supposer une distribution normale multivariée: la majorité de nos variables dans la RDA partielle ne sont pas distribués de manière normale.
* RDA partielle parce qu’on a des covariables qui ne nous intéressent pas mais qui vont influencer les relations entre les autres variables.
* Combine une régression linéaire multiple à une ordination classique (PCA).

**Pourquoi pas MANOVA? Pourquoi pas PCA/analyse non-contrainte?**

## Étape 7. Choix de transformation de la matrice d'abondance: 
Communauté d’amphibiens :

* Problème de doubles zéros
* Données standardisées par abondance totale d’espèce
Transformation nécessaire: mesure de **dissimilarité asymétrique**. Nous avons  2 options:

1. Utiliser les données d’abondance standardisées et les transformer
2. Utiliser les données de présence-absence (binaire) et les transformer

Nous choisissons la 1re option.

* On conserver davantage d’information sur les co-fluctuations des abondances et ça respecte l’idée des investigateurs de donner davantage de poids aux espèces rares dans l’analyse de la communauté amphibienne de ces étangs. 
* Nous utiliserons la transformation de Hellinger, recommandée par de nombreux auteurs pour l’ordination des données d’abondance d’espèces (Legendre et Legendre, 2012) et qui obtient le meilleur R² pour représenter la distance géographique réelle (Legendre et Gallagher, 2001).

## Étape 8.0 Exécuter l'analyse 
A- On commence par transformer les matrices d'abondance. 

```{r hellinger_tranformation, message=FALSE}
library("vegan")
#library(adespatial)
# Faire une matrice de distance Hellinger

# Alternate, two-step computation in vegan: 
amphi.hel= dat%>% 
    subset(oxygene>0) %>% 
  # on enlève les sites avec l'oxygène à 0 ^^
  select(A_Gr:N_Vi) %>% 
  decostand("hellinger")
rownames(amphi.hel) <- subset(dat, oxygene>0)$etang
# On change le nom des lignes pour se retrouver plus tard ^^

head(amphi.hel) %>% round(digits = 2)
```

## Étape 8.0 Exécuter l'analyse 
 
B- On standardize les covariables.

```{r standardize_covariables,message=FALSE}
# Co-variables
topo=dat %>% 
    subset(oxygene>0) %>% 
  select(aire, profondeur) %>% 
  decostand("standardize")

rownames(topo) <- subset(dat, oxygene>0)$etang
# On change le nom des lignes pour se retrouver plus tard ^^
head(topo) %>% round(digits = 2)
```
## Étape 8.0 Exécuter l'analyse 

C- Les variables.
```{r standardize_variables, message=FALSE}
env=dat %>% 
  select( vache_std, conductivite:temperature, invertebres:poissons) %>%
  select(-pH) %>% # on enlève ceci. pas assez fiable. 
  subset(oxygene>0) %>% # erreurs de mesure. pas assez fiable.
  decostand("standardize") %>% 
  mutate(cloture=dat[dat$oxygene>0,]$cloture) %>% data.frame()
# on rattache la variable cloture, qui n'est pas standardisée (catégorique)^^
rownames(env) <- subset(dat, oxygene>0)$etang
# On change le nom des lignes pour se retrouver plus tard ^^
glimpse(env)
```

# RDA partielle | modèle global 

## Étape 8.1 Exécuter l'analyse 
On lance la RDA partielle: 

```{r RDA_launched, message=FALSE}
require(vegan)
(prda=rda(amphi.hel~ . + Condition(topo$aire +topo$profondeur), env) )
```

## Étape 8.1 Exécuter l'analyse 

Donc, *avant* de valider si l'analyse a bien fonctionné, on peu voir dans le output de la RDA que près de 10% de la variance dans l'abondance est expliquée par les covariables, et environ 37% par les variables environnementales. <br> Par contre, ces pourcentages de variance expliqué sont biaisés. Effectivement,lorsqu'on utilise un R² ajusté, la variation expliquée par les variables environnementales tombent à :

```{r RDA_code, message=FALSE}
 (R2adj <- round(RsquareAdj(prda)$adj.r.squared, dig=3))
```
Plus modeste! Aussi, ça concorde un peu plus avec la comparaison des premières valeurs propres de la RDA et de la PCA sur la variation résiduelle, qui sont pas mal similaires. 

## Étape 8.2 Valider l'analyse 

Le test de permutations du modèle global nous indique que le modèle a quand même une certaine valeur, il permet d'expliquer une proportion de la variance dans les communautés d'abondances :
```{r RDA_check_model}
#source("panelutils.R")
anova(prda, permutations = how(nperm = 999))
```

## Étape 8.2 Valider l'analyse 

Le test de permutations sur les axes nous indique que seulement le premier axe de la RDA est *significatif*

```{r RDA_check_axis}
anova(prda, by="axis", permutations = how(nperm = 999))
```

## Étape 8.2 Valider l'analyse 
Le tests de permutations sur les variables explicatives montre que les variables `cloture`, `nitrates`, `vache_std` et `poissons` semblent mieux expliquer la variance dans la matrice d'abondances des amphibiens. C'est aussi envisageable de réduire le nombre de variables dans la RDA dans de prochaines étapes. 

```{r RDA_checks_terms}
anova(prda, by="term", permutations = how(nperm = 999))
```

## Étape 8.2 Valider l'analyse 

Vérifier s'il y a de la collinéarité entre les variables explicatives:
```{r RDA_checks_vif}
vif.cca(prda)
```
Tout est bien. <br>
Jettons un coup d'oeil au triplots maintenant

## Étape 8.3 - Interpréter les résultats : triplots 

```{r RDA_triplots_scale_1, echo=F, message=F, warnings=F, results="hide", fig.dim=c(10, 5), fig.align="center"}

source("triplot.rda.R")

par(mfrow=c(1,2))

triplot.rda(prda, scaling=1, 
            move.origin=c(0, 0), 
            mar.percent=0,
            silent=T, site.sc = "lc", optimum = T)

triplot.rda(prda, scaling=2, 
            move.origin=c(0, 0), 
            mar.percent=0.05,
            silent=T,site.sc = "lc", optimum = T)
```

## 8.3 - Interpréter les résultats : triplots 
* Différences dans la communauté d’amphibiens entre les étangs clôturés et non clôturés
    + Les espèces d’amphibiens rares dans notre système sont favorisées par la présence de clôture
* L’intensité d’utilisation des étangs par les vaches n’influence pas la communauté d’amphibiens

Ces triplots sont très chargés, et beaucoup de variables sont près de l'origine.<br>On peut surement réduire les variables dans la RDA!


## - à compléter ! Modèle parcimonieux 
On prend le backward ou le forward ? On choisi d'inclure les variables issues du backward selection ? <br>
On lance la RDA partielle pour le modèle parcimonieux. 

```{r RDA.parci_launched, message=FALSE, echo=T, eval=F}
#require(vegan)

(rda.parci=rda(amphi.hel~ cloture + poissons + Condition(topo$aire +topo$profondeur), env) )
(R2adj.parci <- round(RsquareAdj(rda.parci)$adj.r.squared, dig=3))

set.seed(1)
anova(rda.parci, permutations = how(nperm = 999))
anova(rda.parci, by="axis", permutations = how(nperm = 999))
anova(rda.parci, by="term", permutations = how(nperm = 999))
vif.cca(rda.parci)

#source("triplot.rda.R")

par(mfrow=c(1,3), 
    mar=c(4,2,1,1))
triplot.rda(rda.parci, scaling=1, 
            move.origin=c(0, 0), 
            mar.percent=0.05,
            silent=T)
triplot.rda(rda.parci, scaling=2, 
            move.origin=c(0, 0), 
            mar.percent=0.05,
            silent=T)
```

## - à compléter. triplot super cool qui donne pas la meme chose qu'avec triplot.rda().

```{r, echo=FALSE, eval=F}

# Triplot joli

perc <- round(100*(summary(rda.parci)$cont$importance[2, 1:2]), 2)

## extract scores - these are coordinates in the RDA space
sc_si <- scores(rda.parci, display=c("lc"), choices=c(1,2))
sc_sp <- scores(rda.parci, display="species", choices=c(1,2))
sc_bp <- scores(rda.parci, display="bp", choices=c(1, 2))


par(mfrow=c(1,3), 
    mar=c(4,2,1,1))
triplot.rda(rda.parci, scaling=1, 
            move.origin=c(0, 0), 
            mar.percent=0.05,
            silent=T)
triplot.rda(rda.parci, scaling=2, 
            move.origin=c(0, 0), 
            mar.percent=0.05,
            silent=T)
# Set up a blank plot with scaling, axes, and labels
plot(rda.parci,
     scaling = 1, # set scaling type 
     type = "none", # this excludes the plotting of any points from the results
     frame = T,
     # set axis limits
     ylim = c(-1,1),
     xlim = c(-.6,.6),
     # label the plot (title, and axes)
     main = "Triplot RDA - scaling 1",
     xlab = paste0("RDA1 (", perc[1], "%)"), 
     ylab = paste0("RDA2 (", perc[2], "%)") , 
     bty="o"
)

# add points for site scores
points(sc_si, 
       pch = 21, # set shape (here, circle with a fill colour)
       col = "black", # outline colour
       bg = "lightgreen", # fill colour
       cex = 1.2) # size

points(sc_si[c("48","90", "126", "150", "209"), ], 
       pch = 21, # set shape (here, circle with a fill colour)
       col = "black", # outline colour
       bg = "darkgreen", # fill colour
       cex = 1.2) # size


# add points for species scores
points(sc_sp, 
       pch = 24, # set shape (here, square with a fill colour)
       col = "orange",
       bg = "orange", 
       cex = 1.2)
# add text labels for species abbreviations
text(sc_sp + c(0.03, 0.09), # adjust text coordinates to avoid overlap with points 
     labels = rownames(sc_sp), 
     col = "orange", 
     font = 2, # bold
     cex = 1)
# add arrows for effects of the expanatory variables
arrows(0,0, # start them from (0,0)
       sc_bp[,1], sc_bp[,2], # end them at the score value
       col = "red", 
       lwd = 1)
# add text labels for arrows
text(x = sc_bp[,1] -0.1, # adjust text coordinate to avoid overlap with arrow tip
     y = sc_bp[,2] - 0.03, 
     labels = rownames(sc_bp), 
     col = "red", 
     cex = 0.9, 
     font = 1)

```

## 8.3 - à compléter ! Partitionner la variance : modèle global vs parcimonieux

```{r var-par-modele-global, echo=F, message=F, warnings=F, results="hide", fig.dim=c(10.5, 5.3), fig.align="left"}
# Explanation of fraction labels (two, three and four explanatory matrices) with optional colours)
#par(mfrow = c(1, 3), mar = c(1, 1, 1, 1))
# showvarparts(2, bg = c("red", "blue"))
# showvarparts(3, bg = c("red", "blue", "yellow")) 
# showvarparts(4, bg = c("red", "blue", "yellow", "green"))

## Partitionner la variance avec toutes les variables du modèle 
vache<-env %>% select(cloture, vache_std)
anim<-env %>% select(invertebres:poissons)
chem<-env %>% select(conductivite:temperature)
topo=topo

amphi.part.all <- varpart(amphi.hel, vache, chem, topo, anim)

par(mfrow= c(1,2),
        mar=c(1,0,2,.2))
plot(amphi.part.all, 
     digits = 3,
     bg = c("red",
            "blue", 
            "yellow", "orange"), 
     Xnames = c("Bétail", "Physico-chimie", "Topo", "Facteurs\nbiotiques"), 
     id.size=0.9, cex=.8)+
    title("Modèle global")



```

## À compléter : interprétation variance 


## Matériel supplémentaire | Effet de quelques variables moins discipliniés

On a voulu vérifier la validité de valeurs plus élevées ou plus basses pour le pH et la température.  <br>
- Premièrement, les valeurs extrèmes de température ont lieu surtout dans les petits étangs peu profonds, donc ce ne sont probablement pas des erreurs de mesures.
- Deuxièmement, les valeurs de pH très élevées et très basses sont louches un peu. P. ex. même avec un pH de 1.88, des poissons ont été observés dans l'étang 155.


```{r extreme_values, eval=FALSE, results="hide", message=F, warnings=F, echo=F}
# Valeurs extrèmes dans les paramètres ####

# temperature
dat %>% 
  subset(temperature>30) %>% 
  select(etang, vache_std, cloture,
         conductivite:temperature,
         aire, profondeur)

dat %>% 
  subset(temperature>30) %>% 
  select(etang, temperature, invertebres:poissons) 

#Donc les valeurs de temperatures élevées sont dans les étangs les plus petits et moins profonds. Rassurant, et ce ne sont probablement pas des valeurs erronées.  Ces deux étangs ont aussi des valeurs de pH dans  le 1er quantile et des concentrations d'oxygène dans le 3e quantile.

# pH

dat %>% 
  subset(pH<4|pH>9) %>% 
  select(etang, vache_std, cloture,
         conductivite:temperature,
         aire, profondeur)

dat %>% 
  subset(pH<4|pH>9) %>% 
  select(etang, pH, invertebres:poissons) 
# Effectivement, dur de dire ce qu'il s'est passé, mais même avec un pH de 1.38, des poissons ont quand même été comptés dans l'étang  # 155!

```

Ensuite, on épargne le code un peu, pour des raisons visuelles, mais on a testé différentes transformation (plusieurs saveurs de log, et aussi la racine carré) sur la température, les poissons et les ecrevisses.  
- Visuellement, ce n'était pas très concluant. On a gardé la valeur de température sur son échelle originale. 
- Par contre, on a modélisé les poissons et écrevisses comme des valeurs binaires, question de vérifier si les conclusions varient beaucoup entre deux modèles globaux.  

```{r check_if_transformations_help, results="hide", message=F, warnings=F, echo=F, eval=F}
# Est-ce que les variables bénéficieraient d'une transformation ? ####

# la température, les poissons et les écrevisses ont beaucoup de 0. on pourrait vérifier l'effet de transformer la variable peut-être ?

par(mfrow=c(2,3))
dat$temperature %>% 
hist(main="temperature")
log(dat$temperature) %>% 
hist(main="ln(temperature)")

log(dat$temperature+0.001) %>% 
hist(main="ln(temperature+0.001)")
log(dat$temperature+0.01) %>% 
hist(main="ln(temperature+0.01)")
log(dat$temperature+0.1) %>% 
hist(main="ln(temperature+0.1)")
log(dat$temperature+1) %>% 
hist(main="ln(temperature+1)")

par(mfrow=c(1,1))

sqrt(dat$temperature) %>% 
hist(main="sqrt(temperature)")
 # temperature, pas évident, mais peu importe les transformations la distribution la plus belle est l'originale. 
# restons avec les valeurs originales.

# Poissons:
par(mfrow=c(2,3))

dat$poissons %>% 
hist(main="poissons")

log(dat$poissons+0.001) %>% 
hist(main="ln(poissons+0.001)")

log(dat$poissons+0.01) %>% 
hist(main="ln(poissons+0.01)")

log(dat$poissons+0.1) %>% 
hist(main="ln(poissons+0.1)")

log(dat$poissons+1) %>% 
hist(main="ln(poissons+1)")

par(mfrow=c(1,1))
sqrt(dat$poissons) %>% 
hist(main="sqrt(poissons)")

#Poissons, pas évident, il y a beaucoup de 0. Peu importe les transformations la distribution la plus belle est l'originale. 
# restons avec les valeurs originale, ou encore on peut créer une variable binaire. 

# ecrevisses:
par(mfrow=c(2,3))

dat$ecrevisses %>% 
hist(main="ecrevisses")

log(dat$ecrevisses+0.001) %>% 
hist(main="ln(ecrevisses+0.001)")

log(dat$ecrevisses+0.01) %>% 
hist(main="ln(ecrevisses+0.01)")

log(dat$ecrevisses+0.1) %>% 
hist(main="ln(ecrevisses+0.1)")

log(dat$ecrevisses+1) %>% 
hist(main="ln(ecrevisses+1)")

sqrt(dat$ecrevisses) %>% 
hist(main="sqrt(ecrevisses)")

 # ecrevisses, pas évident, il y a beaucoup de 0. Peu importe les transformations la distribution ne s'améliore pas beaucoup. restons avec les valeurs originale, ou encore on peut créer une variable binaire.
```

## Matériel supplémentaire | RDA avec variables binaires
```{r scnd_RDA_prepare_data, results="hide", message=F, warnings=F, echo=F}
topo=dat %>% 
    subset(oxygene>0) %>% 
  select(aire, profondeur) %>% 
  decostand("standardize")

rownames(topo) <- subset(dat, oxygene>0)$etang
# On change le nom des lignes pour se retrouver plus tard ^^
env2=dat %>% 
  select( vache_std, 
          conductivite:temperature,
          invertebres) %>%
  mutate(temperature2=as.numeric(ifelse(temperature >29.2, 30.2, temperature))) %>%
  select(-temperature)%>% 
  rename(temperature=temperature2) %>% 
  subset(oxygene>0) %>% 
  decostand("standardize") 

tmp=dat %>% 
  select(cloture, poissons, ecrevisses, oxygene) %>% 
    subset(oxygene>0) %>% 
  mutate(poiss_bin=case_when(poissons>0~ "Yes", 
                             T~"No"),
         ecre_bin=case_when(ecrevisses>0~"Yes", 
                         T~"No")) %>%  
  mutate_at(.vars=c("poiss_bin", "ecre_bin"), factor) %>% 
  select(-poissons, -ecrevisses, -oxygene) 
env2=cbind(tmp, 
          env2) 
rm(tmp)
  # on rattache les variables catégoriques ^^
rownames(env2) <- subset(dat, oxygene>0)$etang
# On change le nom des lignes pour se retrouver plus tard ^^
```
On lance la RDA partielle, avec des variables binaires poiur les poissons et les écrevisses. `env2` a été fait plus haut mais le code est caché. 

```{r scnd_RDA_launched, message=FALSE}
require(vegan)
(prda2=rda(amphi.hel ~ . + Condition(topo$aire +topo$profondeur), env2) )
```

```{r scnd_RDA_checks, eval=F, message=FALSE}
 (R2adj_2 <- round(RsquareAdj(prda2)$adj.r.squared, dig=3))

anova(prda2, permutations = how(nperm = 999))
anova(prda2, by="axis", permutations = how(nperm = 999))
anova(prda2, by="term", permutations = how(nperm = 999))
vif.cca(prda2)

source("triplot.rda.R")

par(mfrow=c(1,2))

triplot.rda(prda2, scaling=1, 
            move.origin=c(0, 0), 
            mar.percent=0.05,
            silent=T)
triplot.rda(prda2, scaling=2, 
            move.origin=c(0, 0), 
            mar.percent=0.05,
            silent=T)

```

## Matériel supplémentaire | Selection backward de variables
En faisant la selection **backward** de la RDA partielle (modèle global), ces variables sont celles qui ressortent: 
```{r, echo=F, eval=T, message=FALSE, results="hide", warning=FALSE}
set.seed(1)
step.back<- ordistep(prda, permutations = how(nperm=499))
```

```{r, echo=TRUE, eval=FALSE}
step.back<- ordistep(prda, permutations = how(nperm=499))
step.back$terminfo$terms %>% attr(., "variables")

# list(nitrates, temperature, poissons, cloture, topo$aire, topo$profondeur)
```
avec un R² ajusté de 
```{r, echo=F, eval=T, message=FALSE, warning=FALSE}
RsquareAdj(step.back)$adj.r.squared
```
et le test de permutations sur les variables enlevées à l'air de ceci:
```{r, echdo=F, eval=T, message=FALSE, warning=FALSE}
step.back$anova
```

et le partitionnement de la variance pour ce modèle : 
```{r}
cloture<-env[,1]
poiss=env$poissons
physico=env[,c("temperature", "nitrates")]

amphi.parsi <- varpart(amphi.hel, cloture, physico, poiss, topo)
par(mar=c(1,.2,2,0))

plot(amphi.parsi, 
     digits=3, bg = c("red", 
                      "orange", 
                      "yellow", "blue"),
     Xnames = c("Clôture","Physico-chimie", "Poissons", "Topo"), 
     id.size=0.9, cex=.8)+
  title("Modèle le plus parsimonieux")
```