---
title: "R Script for 'The contrasted impacts of grasshoppers on soil microbial activities in function of ecosystem productivity and herbivore diet' "
date: "July 2022"
author: "Sébastien Ibanez"
output:
  html_document:
    number_sections: yes
    toc: true  
  pdf_document:
    number_sections: yes
    toc: true  
  
---

```{r, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(ade4)
library(lmerTest); lmer.df <- "satterthwaite"
library(vegan)
```

# Data preparation

```{r}
data <- read.csv("data.csv")
colnames(data)[9] <- "EEP"
data$EEP <- data$EEP + 1e-6 # detection threshold
```

Preliminary analysis showed that several variables require log-transformation. 

```{r}
data$log_PNM <- log(data$nmp)
data$log_EEC <- log(data$C_ez)
data$log_EEN <- log(data$N_ez)
data$log_EEP <- log(data$EEP)
data$log_ratio <- log(data$ratio)

data$log_som <- log(data$som+1)
data$log_tdn <- log(data$tdn+1)
data$log_pho <- log(data$pho+1)
data$log_soil_CN <- log(data$CN+1)
```


The following sections correspond to the results section of the article, except the section "Leaf quality and consumption" which is not detailed here since it concerns preliminary results.

# Relations between vegetation, soil microbial and abiotic characteristics

```{r}
var_bio <- c("log_PNM","mbn","log_EEC","log_EEN","log_EEP","log_ratio")
pca_bio <- dudi.pca(df = data[, var_bio], scannf = F)
var_abio <- c("hum","log_som","pH","log_tdn","log_pho","log_soil_CN")
pca_abio <- dudi.pca(df = data[, var_abio], scannf = F)
var_veget <- c("biom","plantCN","N_B","shan","forbs")
pca_veget <- dudi.pca(df = data[, var_veget], scannf = F)
```

## Coinertia between microbial and abiotic variables

```{r}
coin_bio_abio <- coinertia(dudiX=pca_bio, dudiY=pca_abio, scannf = F)
rvtest_bio_abio <- RV.rtest(pca_bio$tab, pca_abio$tab, 99999)
rvtest_bio_abio
```

## Coinertia between abiotic and vegetation variables

```{r}
coin_abio_veget <- coinertia(dudiX=pca_abio, dudiY=pca_veget, scannf = F)
rvtest_abio_veget <- RV.rtest(pca_abio$tab, pca_veget$tab, 99999)
rvtest_abio_veget
```

## Coinertia between microbial and vegetation variables

```{r}
coin_bio_veget <- coinertia(dudiX=pca_bio, dudiY=pca_veget, scannf = F)
rvtest_bio_veget <- RV.rtest(pca_bio$tab, pca_veget$tab, 99999)
rvtest_bio_veget
```

## RDA

```{r}
varpart(Y=data[, var_bio], ~ site, ~ com, data=data)
varpart(Y=data[, var_abio], ~ site, ~ com, data=data)
varpart(Y=data[, var_veget], ~ site, ~ com, data=data)
```


# Effect of herbivory on microbial characteristics

## Calculation of standardized mean difference

```{r}
smd_data <- data[, c("meso","site","herb","com", var_bio)]

for(i in 1:nrow(data)){
  inter1 <- data[i,"site"]
  inter2 <- subset(data, site==inter1 & herb=="Control")
  smd_data[i,var_bio] <- (data[i,var_bio]-inter2[var_bio])/apply(data[data$site==inter2$site, var_bio],2,sd)
}

smd_data <- subset(smd_data, herb!="Control")

colnames(smd_data)[5:10] <- paste("smd_",colnames(smd_data)[5:10], sep="")
smd_data <- merge(smd_data, data[,c("meso","site","herb","com", var_abio, var_veget)])
rownames(smd_data) <- rownames(data[data$herb!="Control" , ])
```

## Coinertia

First, since the standardized responses are calculated for cages containing herbivores, control cages need to be excluded from the pcas.

```{r}
dudi_abio_herb <- dudi.pca(data[data$herb!="Control" , var_abio], scannf = F)
dudi_veget_herb <- dudi.pca(data[data$herb!="Control" , var_veget], scannf = F)
```

Second, the pca of the smd are calculated

```{r}
dudi_bio_smd <- dudi.pca(smd_data[ , 5:10], scannf = F)
```

Finally, the coinertias are performed.

### Coinertia response to herbivory vs abiotic soil characteristics

```{r}
coin_smdbio_abio <- coinertia(dudiX=dudi_abio_herb, dudiY=dudi_bio_smd, scannf = F)
rvtest_smdbio_abio <- RV.rtest(dudi_abio_herb$tab, dudi_bio_smd$tab, 99999)
rvtest_smdbio_abio
```

### Coinertia response to herbivory vs vegetation characteristics

```{r}
coin_smdbio_veget <- coinertia(dudiX=dudi_veget_herb, dudiY=dudi_bio_smd, scannf = F)
rvtest_smdbio_veget <- RV.rtest(dudi_veget_herb$tab, dudi_bio_smd$tab, 99999)
rvtest_smdbio_veget
```

## Univariate regressions

### Predictors : soil abiotic characteristics

```{r}
anova(lmer(smd_log_PNM~log_som*herb+(1|site), data=smd_data))
anova(lmer(smd_log_EEN~log_som*herb+(1|site), data=smd_data))
anova(lmer(smd_log_EEC~log_som*herb+(1|site), data=smd_data))
anova(lmer(smd_log_EEP~log_som*herb+(1|site), data=smd_data))
anova(lmer(smd_log_EEP~pH*herb+(1|site), data=smd_data))
anova(lmer(smd_mbn~log_som*herb+(1|site), data=smd_data))
anova(lmer(smd_mbn~log_tdn*herb+(1|site), data=smd_data))
```

### Predictors : vegetation characteristics

```{r}
anova(lmer(smd_log_EEC~N_B*herb+(1|site), data=smd_data))
anova(lmer(smd_log_EEN~N_B*herb+(1|site), data=smd_data))
anova(lmer(smd_log_PNM~N_B*herb+(1|site), data=smd_data))
anova(lmer(smd_log_EEP~N_B*herb+(1|site), data=smd_data))

anova(lmer(smd_log_EEC~biom*herb+(1|site), data=smd_data))
anova(lmer(smd_log_EEN~biom*herb+(1|site), data=smd_data))
anova(lmer(smd_log_PNM~biom*herb+(1|site), data=smd_data))
anova(lmer(smd_log_EEP~biom*herb+(1|site), data=smd_data))
```

