#### Packages --------
library(vegan)
library(FD)
library(car)
library(emmeans)
library(multcompView)
library(DHARMa)
library(MuMIn)
library(multcomp)
library(glmmTMB)
library(multifunc)


### Preparation of the table with botanical observations ----
quadrat <- read.table("/Bota.csv", header=T, as.is=T, stringsAsFactors=F, sep=";", dec=",")
rownames(quadrat) <- quadrat[,1]
quadrat <- quadrat[,-1]
quadrat[,3:ncol(quadrat)] <- as.data.frame(lapply(quadrat[,3:ncol(quadrat)], as.numeric))


### Preparation of the trait matrix ----
data <- read.csv("/Traits.csv", sep=";", header=T, as.is=T, stringsAsFactors=F)
rownames(data) <- data[,1]
data <- data[,-1]
colnames(data)

trait <- cbind(as.data.frame(as.factor(data$UV_pattern) ), 
               as.data.frame(as.factor(data$couleur)),  # flower colour 
               as.data.frame(data$debut_flo),           # calendar week of flower onset
               as.data.frame(log(data$duree_flo)),      # duration of flowering, in week numbers
               as.data.frame(log(data$hauteur_repro)),  # reproductve height in m
               as.data.frame(as.factor(data$distri_feuilles)),  # leaf distribution
               as.data.frame(as.numeric(data$NEF)),     # presence of extrafloral nectar
               as.data.frame(as.numeric(data$quantite_nectar)),   # index of nectar amount
               as.data.frame(log(as.numeric(data$taille_flo))),  # flower or infloresence size in mm
               as.data.frame(log(as.numeric(data$prof_nectar)+1)),   # nectar depth in mm
               as.data.frame(as.factor(data$pilosite_feuille)),    # leaf pilosity
               as.data.frame(as.numeric(data$SLA)),  # specific leaf area
               as.data.frame(as.numeric(data$LDMC))   # leaf dry matter content
              )

#### Calculation of functional diversity ----
trait <- trait[colnames(quadrat)[-c(1:2)],]  
fdiv <- dbFD(trait, quadrat[,3:ncol(quadrat)], corr=c("cailliez"), CWM.type="all")
fdiv <- as.data.frame(fdiv)



#### Null models to stardardize functional diversity with regards to species richness -----
# To generate null FDis values for each community by shuffling the names of species on the trait data matrix, we can first write this function to shuffle names and then calculate a null FDis value using the dbFD() function.
dbfd.shuff <- function(x) {
  rownames(x) <- sample(rownames(x), length(rownames(x)), replace=F  )
  x <- x[rownames(trait),]
  dbFD(x, quadrat[,3:ncol(quadrat)], corr=c("cailliez"))$FDis
}
# We now replicate this function several times and combine it with the observed output
nrep <- 1000
obs.null.output <- cbind(dbFD(trait, quadrat[,3:ncol(quadrat)], corr=c("cailliez"))$FDis,  replicate(nrep, dbfd.shuff(trait)))
ses.all <- (fdiv$FDis - apply(obs.null.output, MARGIN=1, mean)) / apply(obs.null.output, MARGIN=1, sd)
fdiv 
<- cbind(fdiv,ses.all)


#### Linear models ----

# Preparation of the data table
pucerons <- read.table("/Aphids.csv", header=T, as.is=T, stringsAsFactors=F, sep=";")
pucerons <- merge(pucerons, fdiv, by.x=2, by.y=0, all.x=T)
pucerons$Distance <- as.factor(pucerons$Distance)

# Selection of the crop and of the date
sub <- subset(pucerons, Crop=="wheat") # selection of the wheat crop
sub <- subset(sub, Date=="T1") # selection of the first counting date

# Effect of plant mixtures on aphid number
mod <- glmmTMB( Aphid_number ~ Mixture*Distance + (1|Block/Number), data=sub, na.action="na.fail", family=nbinom1) 
# Effect of plant mixtures on predator-prey ratios
mod <- glmmTMB( Pred_prey_ratio ~ Mixture * Distance + (1|Block/Number), data=sub, na.action="na.fail", family=nbinom1()) 
# Effect of plant mixtures on parasitism rates
mod <- glmmTMB( cbind(Mummies,Aphid_number) ~ Mixture * Distance + (1|Block/Number), data=sub, na.action="na.fail", family=binomial) 

# Diagnostic and summary of the models
plot(simulateResiduals(mod))
summary(mod)
Anova(mod, type="II")

# Post-hoc tests
moy <- lsmeans(mod, ~ Mixture) ; cld(moy, Letters = c(LETTERS))
moy <- lsmeans(mod, ~ Distance) ; cld(moy, Letters = c(LETTERS))
moy <- lsmeans(mod, ~ Mixture*Distance) ; cld(moy, Letters = c(LETTERS))


# Effect of plant composition and diversity within flower strips on predator-prey ratios
mod <- glmmTMB( Pred_prey_ratio ~ scale(Legumes)*Distance + scale(Nectar)*Distance + scale(ses.all)*Distance+ scale(nbsp)*Distance + scale(ses.all)*scale(Nectar) + (1|Block/Number), data=sub, na.action="na.fail", family=nbinom1)
# Effect of plant composition and diversity within flower strips on parasitism rates
mod <- glmmTMB( cbind(Mummies,Aphid_number) ~ scale(Legumes)*Distance + scale(Nectar)*Distance + scale(ses.all)*Distance+ scale(nbsp)*Distance + scale(ses.all)*scale(Nectar) + (1|Block/Number), data=sub, na.action="na.fail", family=binomial)
# Effect of plant composition and diversity of flower strips on the growth rates of aphid populations
mod <- glmmTMB( Growth_rate ~ scale(Legumes)*Distance + scale(Nectar)*Distance + scale(ses.all)*Distance + scale(nbsp)*Distance + scale(ses.all)*scale(Nectar) + (1|Block/Number), data=sub, na.action="na.fail", family=gaussian) 

# Diagnostic and summary of the models
plot(simulateResiduals(mod))
summary(mod)
Anova(mod, type="II")

# Multimodel inference and model averaging
select <- dredge(mod, beta="none", evaluate=T, rank="AIC") ; select
avg <- model.avg(select, delta<4) ; summary(avg) ; sw(avg)
r.squaredGLMM(mod)