library(readxl)
library(ggplot2)
library(Rmisc)
library(scales)
library(nlme)
library(lme4)
library(vegan)
library(multcomp)
library(car)
library(ggrepel)
library(devtools)
library(pairwiseAdonis)
library(grid)
library(reshape2)
library(DHARMa)
library(FD)
library(glmmTMB)
library(devtools)
install_github("zdealveindy/weimea")
library(weimea)
library(tidyr)



###Interpretation plant communities by dissimilarity

##abundance data of species (only species in 5% of the plots)
varespec5 <- read_excel("data_paper_AndreuCera.xlsx", sheet="vegetation_data_5%")
rownames(varespec5)<-varespec5$ID

##Soil feautres
#we added values of PC1, PC2, PC3 on environmental_data file
soil <- read_excel("data_paper_AndreuCera.xlsx", sheet="soil_data")
rownames(soil)<-soil$ID

pca<-rda(soil[,4:11], scale=TRUE)
summary(pca)
sites<-scores(pca, display = "sites",scaling=1, choices = 1:3)
species<-scores(pca, display = "species",scaling=1, choices = 1:3)
arrows<-data.frame(species,rownames(species))
colnames(arrows)<-c("PC1","PC2","PC3","soil")
grad<-data.frame(sites, soil)
rownames(grad)<-grad$ID

soil_melt<-melt(grad, id=c("ID","Locality","Grazing"))
soil_means <- summarySE(soil_melt, measurevar="value", groupvars=c("variable","Locality","Grazing"))


write.csv2(soil_means,"grad.csv")

##PERMANOVA
varechem <- read_excel("data_paper_AndreuCera.xlsx", sheet="environmental_data") ##can be grad (data frame)
varespec5 <- read_excel("data_paper_AndreuCera.xlsx", sheet="vegetation_data_5%")
varespec5<-data.frame(varespec5)
abs<-sqrt(varespec5[,5:34])

perm <- how(nperm = 1000)
setBlocks(perm) <- with(varespec5, locality)
adonis2(abs~PC1+grazing, data=varechem, perm=perm, method="bray") # locality as strata

##Correlaction environmental variables with plant composition
vare.mds <- metaMDS(varespec5[,5:34], trace = FALSE)
ef<- envfit(vare.mds, varechem[,5:12], permu = 999)
ef

#Figure S2
points<-data.frame(varespec5$grazing,vare.mds$points)
rownames(points)<-varespec5$ID
colnames(points)<-c("pert","NMDS1","NMDS2")
tgc1 <- summarySE(points, measurevar=c("NMDS1"), groupvars=c("pert"))
tgc2 <- summarySE(points, measurevar=c("NMDS2"), groupvars=c("pert"))
graphic<-data.frame(tgc1$pert,tgc1$NMDS1,tgc1$se,tgc2$NMDS2,tgc2$se)
colnames(graphic)<-c("pert","NMDS1","se1","NMDS2","se2")
vectors<-ef$vectors
var<-colnames(varechem[,5:12])
arrows<-data.frame(var,vectors$arrows)

graphic$pert<-c("High","Medium","Low")
m<-ggplot(graphic, aes(x=NMDS1,y=NMDS2)) #plot centroids
m+geom_segment(data=arrows,aes(x=0,xend=NMDS1/2.3,y=0,yend=NMDS2/2.3),arrow = arrow(length = unit(0.5, "cm")),colour="grey60")+
  annotate("text", x = ((arrows$NMDS1/2.2)), y = ((arrows$NMDS2/2.2)),label = arrows$var)+
  geom_point(aes(colour=pert), size=3)+
  geom_text_repel(aes(label=pert), size = 4,nudge_x=-0.01, nudge_y=0.02) +
  labs(title =element_blank(), x = "NMDS1", y = "NMDS2")+
  geom_errorbarh(aes(xmin=NMDS1-se1, xmax=NMDS1+se1), height=.01, colour="grey")+
  geom_errorbar(aes(ymin=NMDS2-se2, ymax=NMDS2+se2), width=.01, colour="grey")+
  scale_shape_manual(labels = c("High", "Medium", "Low"), values = c(0,15,1)) +
  scale_colour_grey() +
  theme( 
    legend.position="none",
    legend.title = element_blank(), 
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.background = element_rect(fill = "white", colour = "grey50"),
  )



###Community indices
##How to calculate CWM_fixed of leaf traits
varespec <- read_excel("data_paper_AndreuCera.xlsx", sheet="vegetation_data")
traits<-read_excel("data_paper_AndreuCera.xlsx", sheet="leaf_traits")
traits$species<-as.factor(traits$species)
spp<-levels(traits$species)
veg_subset<-data.frame(subset(varespec, select=c("locality","grazing", "plot","ID")),subset(varespec,select=spp)) #only subset species with leaf traits

traits_melt<-melt(traits, measure.vars=c("C","N","S"))
traits_means <- summarySE(na.omit(traits_melt), measurevar="value", groupvars=c("variable","species"))

rownames(veg_subset)<-veg_subset$ID


nut<-c("C","N","S", "plot")
CWM=data.frame(matrix(ncol=4,nrow=315))
colnames(CWM)<-nut


for (i in 1:length(nut)){ 
  traits4 <- traits_means[which(traits_means$variable==nut[i]), ]
  B<-data.frame(traits4[,4])
  A<-veg_subset[,5:18]
  colnames(A)<-traits4$species
  rownames(A)<-veg_subset$ID
  rownames(B)<-traits4$species
  ex1 <- dbFD(B, A,w.abun = TRUE,calc.FRic=TRUE)
  tabla<-data.frame(ex1$CWM)
  ta2<-data.frame(tabla,rownames(tabla))
  colnames(ta2)<-c("aff","id")
  CWM[,i]<-ta2$aff
  CWM[,4]<-ta2$id
  rownames(CWM)<-ta2[,2]
}
  

CWM_fixed_leaf<-data.frame(CWM %>% separate(plot,c("locality","grazing","plot"), sep="_"))
    

##How to calculate CWM_specific of leaf traits
traits<-read_excel("data_paper_AndreuCera.xlsx", sheet="leaf_traits")
traits$species<-as.factor(traits$species)
spp<-levels(traits$species)
veg_subset<-data.frame(subset(varespec, select=c("locality","grazing", "plot","ID")),subset(varespec,select=spp)) #only subset species with leaf traits
traits_melt<-melt(traits, measure.vars=c("C","N","S"))
traits_means <- summarySE(na.omit(traits_melt), measurevar="value", groupvars=c("variable","species","grazing"))

#High
veg_subset_G<-veg_subset[which(veg_subset$grazing=="High"),] 
traits_means_G<-traits_means[which(traits_means$grazing=="High"),]
rownames(veg_subset_G)<-veg_subset_G$ID
nut<-c("C","N","S","plot")
CWM=data.frame(matrix(ncol=4,nrow=105))
colnames(CWM)<-nut


for (i in 1:length(nut)){ 
  traits4 <- traits_means_G[which(traits_means_G$variable==nut[i]), ]
  B<-data.frame(traits4[,5])
  A<-veg_subset_G[,5:18]
  colnames(A)<-traits4$species
  rownames(A)<-veg_subset_G$ID
  rownames(B)<-traits4$species
  ex1 <- dbFD(B, A,w.abun = TRUE,calc.FRic=TRUE)
  tabla<-data.frame(ex1$CWM)
  ta2<-data.frame(tabla,rownames(tabla))
  colnames(ta2)<-c("aff","id")
  CWM[,i]<-ta2$aff
  CWM[,4]<-ta2$id
  rownames(CWM)<-ta2[,2]
}

CWM_High_leaf<-data.frame(CWM %>% separate(plot,c("locality","grazing","plot"), sep="_"))

#Medium
veg_subset_G<-veg_subset[which(veg_subset$grazing=="Medium"),] 
traits_means_G<-traits_means[which(traits_means$grazing=="Medium"),]
rownames(veg_subset_G)<-veg_subset_G$ID
nut<-c("C","N","S", "plot")
CWM=data.frame(matrix(ncol=4,nrow=105))
colnames(CWM)<-nut


for (i in 1:length(nut)){ 
  traits4 <- traits_means_G[which(traits_means_G$variable==nut[i]), ]
  B<-data.frame(traits4[,5])
  A<-veg_subset_G[,5:18]
  colnames(A)<-traits4$species
  rownames(A)<-veg_subset_G$ID
  rownames(B)<-traits4$species
  ex1 <- dbFD(B, A,w.abun = TRUE,calc.FRic=TRUE)
  tabla<-data.frame(ex1$CWM)
  ta2<-data.frame(tabla,rownames(tabla))
  colnames(ta2)<-c("aff","id")
  CWM[,i]<-ta2$aff
  CWM[,4]<-ta2$id
  rownames(CWM)<-ta2[,2]
}

CWM_Medium_leaf<-data.frame(CWM %>% separate(plot,c("locality","grazing","plot"), sep="_"))

#Low
veg_subset_G<-veg_subset[which(veg_subset$grazing=="Low"),] 
traits_means_G<-traits_means[which(traits_means$grazing=="Low"),]
rownames(veg_subset_G)<-veg_subset_G$ID
nut<-c("C","N","S", "plot")
CWM=data.frame(matrix(ncol=4,nrow=105))
colnames(CWM)<-nut


for (i in 1:length(nut)){ 
  traits4 <- traits_means_G[which(traits_means_G$variable==nut[i]), ]
  B<-data.frame(traits4[,5])
  A<-veg_subset_G[,5:18]
  colnames(A)<-traits4$species
  rownames(A)<-veg_subset_G$ID
  rownames(B)<-traits4$species
  ex1 <- dbFD(B, A,w.abun = TRUE,calc.FRic=TRUE)
  tabla<-data.frame(ex1$CWM)
  ta2<-data.frame(tabla,rownames(tabla))
  colnames(ta2)<-c("aff","id")
  CWM[,i]<-ta2$aff
  CWM[,4]<-ta2$id
  rownames(CWM)<-ta2[,2]
}

CWM_Low_leaf<-data.frame(CWM %>% separate(plot,c("locality","grazing","plot"), sep="_"))

# Bind High, Medium and low together in the same data.frame

CWM_Specific_leaf<-rbind(CWM_High_leaf,CWM_Medium_leaf,CWM_Low_leaf)


##CWM_ITV of leaf traits
#important to reorder columns to avoid mistakes
attach(CWM_Specific_leaf)
CWM_Specific_leaf<-CWM_Specific_leaf[order(locality, grazing, plot),]
detach(CWM_Specific_leaf)
attach(CWM_fixed_leaf)
CWM_fixed_leaf<-CWM_fixed_leaf[order(locality, grazing, plot),]
detach(CWM_fixed_leaf)
CWM_ITV_leaf<-data.frame(CWM_Specific_leaf$C-CWM_fixed_leaf$C,CWM_Specific_leaf$N-CWM_fixed_leaf$N,CWM_Specific_leaf$S-CWM_fixed_leaf$S, CWM_Specific_leaf$locality, CWM_Specific_leaf$grazing,CWM_Specific_leaf$plot)
colnames(CWM_ITV_leaf)<-colnames(CWM_Specific_leaf)

##RAO
#How to calculate Rao_fixed of leaf traits
varespec <- read_excel("data_paper_AndreuCera.xlsx", sheet="vegetation_data")
traits<-read_excel("data_paper_AndreuCera.xlsx", sheet="leaf_traits")
traits$species<-as.factor(traits$species)
spp<-levels(traits$species)
veg_subset<-data.frame(subset(varespec, select=c("locality","grazing", "plot","ID")),subset(varespec,select=spp)) #only subset species with leaf traits
traits_melt<-melt(traits, measure.vars=c("C","N","S"))
traits_means <- summarySE(na.omit(traits_melt), measurevar="value", groupvars=c("variable","species"))
rownames(veg_subset)<-veg_subset$ID
nut<-c("C","N","S")
Rao=data.frame(matrix(ncol=3,nrow=315))
colnames(Rao)<-nut


for (i in 1:length(nut)){ 
  traits4 <- traits_means[which(traits_means$variable==nut[i]), ]
  samp<-data.frame(veg_subset[,5:18])
  colnames(samp)<-traits4$species
  rownames(samp)<-veg_subset$ID
  dis<-vegdist(traits4$value, method="gower")
  df.dist=as.matrix(dis, labels=TRUE)
  colnames(df.dist) <- rownames(df.dist) <- traits4$species
  a<-melodic(samp,df.dist, type="both")  # melodic functions was extracted from De Bello et al. 2016
  b<-data.frame(a$abundance)
  Rao[,i]<-b$rao
  Rao[,4]<-veg_subset$ID
  
}

Rao_fixed_leaf<-data.frame(Rao %>% separate(V4,c("locality","grazing","plot"), sep="_"))



# How to calculate Rao_Specific of leaf traits
traits<-read_excel("data_paper_AndreuCera.xlsx", sheet="leaf_traits")
traits$species<-as.factor(traits$species)
spp<-levels(traits$species)
veg_subset<-data.frame(subset(varespec, select=c("locality","grazing", "plot","ID")),subset(varespec,select=spp)) #only subset species with leaf traits
traits_melt<-melt(traits, measure.vars=c("C","N","S"))
traits_means <- summarySE(na.omit(traits_melt), measurevar="value", groupvars=c("variable","species","grazing"))

#High
veg_subset_G<-veg_subset[which(veg_subset$grazing=="High"),] 
traits_means_G<-traits_means[which(traits_means$grazing=="High"),]
rownames(veg_subset_G)<-veg_subset_G$ID
nut<-c("C","N","S")
Rao=data.frame(matrix(ncol=3,nrow=105))
colnames(Rao)<-nut


for (i in 1:length(nut)){ 
  traits4 <- traits_means_G[which(traits_means_G$variable==nut[i]), ]
  samp<-data.frame(veg_subset_G[,5:18])
  colnames(samp)<-traits4$species
  rownames(samp)<-veg_subset_G$ID
  dis<-vegdist(traits4$value, method="gower")
  df.dist=as.matrix(dis, labels=TRUE)
  colnames(df.dist) <- rownames(df.dist) <- traits4$species
  a<-melodic(samp,df.dist, type="both")  # melodic functions was extracted from De Bello et al. 2016
  b<-data.frame(a$abundance)
  Rao[,i]<-b$rao
  Rao[,4]<-veg_subset_G$ID
  
}

Rao_High_leaf<-data.frame(Rao %>% separate(V4,c("locality","grazing","plot"), sep="_"))

#Medium
veg_subset_G<-veg_subset[which(veg_subset$grazing=="Medium"),] 
traits_means_G<-traits_means[which(traits_means$grazing=="Medium"),]
rownames(veg_subset_G)<-veg_subset_G$ID
nut<-c("C","N","S")
Rao=data.frame(matrix(ncol=3,nrow=105))
colnames(Rao)<-nut


for (i in 1:length(nut)){ 
  traits4 <- traits_means_G[which(traits_means_G$variable==nut[i]), ]
  samp<-data.frame(veg_subset_G[,5:18])
  colnames(samp)<-traits4$species
  rownames(samp)<-veg_subset_G$ID
  dis<-vegdist(traits4$value, method="gower")
  df.dist=as.matrix(dis, labels=TRUE)
  colnames(df.dist) <- rownames(df.dist) <- traits4$species
  a<-melodic(samp,df.dist, type="both")  # melodic functions was extracted from De Bello et al. 2016
  b<-data.frame(a$abundance)
  Rao[,i]<-b$rao
  Rao[,4]<-veg_subset_G$ID
  
}

Rao_Medium_leaf<-data.frame(Rao %>% separate(V4,c("locality","grazing","plot"), sep="_"))

#Low
veg_subset_G<-veg_subset[which(veg_subset$grazing=="Low"),] 
traits_means_G<-traits_means[which(traits_means$grazing=="Low"),]
rownames(veg_subset_G)<-veg_subset_G$ID
nut<-c("C","N","S")
Rao=data.frame(matrix(ncol=3,nrow=105))
colnames(Rao)<-nut


for (i in 1:length(nut)){ 
  traits4 <- traits_means_G[which(traits_means_G$variable==nut[i]), ]
  samp<-data.frame(veg_subset_G[,5:18])
  colnames(samp)<-traits4$species
  rownames(samp)<-veg_subset_G$ID
  dis<-vegdist(traits4$value, method="gower")
  df.dist=as.matrix(dis, labels=TRUE)
  colnames(df.dist) <- rownames(df.dist) <- traits4$species
  a<-melodic(samp,df.dist, type="both")  # melodic functions was extracted from De Bello et al. 2016
  b<-data.frame(a$abundance)
  Rao[,i]<-b$rao
  Rao[,4]<-veg_subset_G$ID
  
}

Rao_Low_leaf<-data.frame(Rao %>% separate(V4,c("locality","grazing","plot"), sep="_"))

# Bind High, Medium and low together in the same data frame

Rao_Specific_leaf<-rbind(Rao_High_leaf,Rao_Medium_leaf,Rao_Low_leaf)


##Rao_ITV of leaf traits
attach(Rao_Specific_leaf)
Rao_Specific_leaf<-Rao_Specific_leaf[order(locality, grazing, plot),]
detach(Rao_Specific_leaf)
attach(Rao_fixed_leaf)
Rao_fixed_leaf<-Rao_fixed_leaf[order(locality, grazing, plot),]
detach(Rao_fixed_leaf)
Rao_ITV_leaf<-data.frame(Rao_Specific_leaf$C-Rao_fixed_leaf$C,Rao_Specific_leaf$N-Rao_fixed_leaf$N,Rao_Specific_leaf$S-Rao_fixed_leaf$S, Rao_Specific_leaf$locality, Rao_Specific_leaf$grazing,Rao_Specific_leaf$plot)
colnames(Rao_ITV_leaf)<-colnames(Rao_Specific_leaf)


###How to calculate CWM Gypsophilic value
varespec <- read_excel("data_paper_AndreuCera.xlsx", sheet="vegetation_data")
GV<-read_excel("data_paper_AndreuCera.xlsx", sheet="gypsophilic_value")

#CWM
B<-data.frame(GV[,2])
A<-varespec[,5:56]
colnames(A)<-GV$spp
rownames(A)<-varespec$ID
rownames(B)<-GV$spp
ex1 <- dbFD(B, A,w.abun = TRUE,calc.FRic=TRUE)
CWM_GV<-data.frame(ex1$CWM, varespec[,1:4])
attach(CWM_GV)
CWM_GV<-CWM_GV[order(locality, grazing, plot),]
detach(CWM_GV)

###How to calculate Simpson index
samp<-data.frame(varespec[,5:56])
colnames(samp)<-GV$spp
rownames(samp)<-varespec$ID
dis<-vegdist(GV[,2], method="gower")
df.dist=as.matrix(dis, labels=TRUE)
colnames(df.dist) <- rownames(df.dist) <- GV$spp
a<-melodic(samp,df.dist, type="both")
b<-data.frame(a$abundance)

Simpson<-data.frame(b$simpson,varespec[,1:4])
attach(Simpson)
Simpson<-Simpson[order(locality, grazing, plot),]
detach(Simpson)

###Barplots vegetation properties, CWM and FD specific traits
#Vegetation properties Figure 1
varechem <- read_excel("data_paper_AndreuCera.xlsx", sheet="environmental_data")
attach(varechem)
varechem<-varechem[order(locality, grazing, plot),]
detach(varechem)
prop<-data.frame(CWM_GV,Simpson$b.simpson,varechem$cover, varechem$canopy_height)
colnames(prop)<-c("GV", "locality","grazing","plot", "ID","Simpson","Vegetation_cover","Canopy_height")

data1<-melt(prop, measure.vars=c("GV","Simpson","Vegetation_cover","Canopy_height"))
tgc1 <- summarySE(na.omit(data1), measurevar="value", groupvars=c("variable","grazing"))

traits<-c(
  "GV"="'Gypsophilic value'",
  "Simpson"="'Simpson'",
  "Vegetation_cover"="'Vegetation cover (%)'",
  "Canopy_height"="'Canopy height (cm)'")

s<-ggplot(data=tgc1, aes(x=grazing, y=value, fill=grazing))
s+geom_bar(stat="identity", colour="black", position=position_dodge()) +
  scale_x_discrete(limits=c("Low","Medium","High"))+
  labs(y=element_blank(), x=element_blank())+
  scale_fill_grey()+
  scale_y_continuous(expand = expand_scale(mult = c(0, .18))) +
  geom_errorbar(aes(ymin=value-se, ymax=value+se), width=.5,position=position_dodge(.9))+
  theme(
    legend.position="none",
    strip.text.x = element_text(size=15, face="bold"),
    strip.background = element_rect(colour="black", fill="white"),
    axis.text=element_text(size=15),
    legend.title = element_blank(), 
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.background = element_rect(fill = "white", colour = "grey50"))+
  facet_wrap(~variable, scales = "free_y", nrow=1, labeller=as_labeller(traits,label_parsed))


#CWM Figure 2a
data_CWM<-data.frame(CWM_Specific_leaf)
colnames(data_CWM)<-c("C","N","S","locality","grazing","plot")
data_CWM$C<-data_CWM$C*10 #to convert from % to mg/g
data_CWM$N<-data_CWM$N*10
data_CWM$S<-data_CWM$S*10


data1<-melt(data_CWM, measure.vars=c("C","N","S"))
tgc1 <- summarySE(na.omit(data1), measurevar="value", groupvars=c("variable","grazing"))

traits<-c(
  "C"="'Leaf C (mg/g)'",
  "N"="'Leaf N (mg/g)'",
  "S"="'Leaf S (mg/g)'")

s<-ggplot(data=tgc1, aes(x=grazing, y=value, fill=grazing))
s+geom_bar(stat="identity", colour="black", position=position_dodge()) +
  scale_x_discrete(limits=c("W","M","H"),labels=c("Low","Medium","High"))+
  labs(y=element_blank(), x=element_blank())+
  scale_fill_grey()+
  scale_y_continuous(expand = expand_scale(mult = c(0, .18))) +
  geom_errorbar(aes(ymin=value-se, ymax=value+se), width=.5,position=position_dodge(.9))+
  theme(
    legend.position="none",
    strip.text.x = element_text(size=15, face="bold"),
    strip.background = element_rect(colour="black", fill="white"),
    axis.text=element_text(size=15),
    legend.title = element_blank(), 
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.background = element_rect(fill = "white", colour = "grey50"))+
  facet_wrap(~variable, scales = "free_y", nrow=1, labeller=as_labeller(traits,label_parsed))

#FD Figure 2b
data_Rao<-data.frame(Rao_Specific_leaf)
colnames(data_Rao)<-c("C","N","S","locality","grazing","plot")

data1<-melt(data_Rao,measure.vars=c("C","N","S"))
tgc1 <- summarySE(na.omit(data1), measurevar="value", groupvars=c("variable","grazing"))

traits<-c(
  "C"="'Leaf C'",
  "N"="'Leaf N'",
  "S"= "'Leaf S'")

s<-ggplot(data=tgc1, aes(x=grazing, y=value, fill=grazing))
s+geom_bar(stat="identity", colour="black", position=position_dodge()) +
  scale_x_discrete(limits=c("W","M","H"),labels=c("Low","Medium","High"))+
  labs(y=element_blank(), x=element_blank())+
  scale_fill_grey()+
  scale_y_continuous(expand = expand_scale(mult = c(0, .18))) +
  geom_errorbar(aes(ymin=value-se, ymax=value+se), width=.5,position=position_dodge(.9))+
  theme(
    legend.position="none",
    strip.text.x = element_text(size=15, face="bold"),
    strip.background = element_rect(colour="black", fill="white"),
    axis.text=element_text(size=15),
    legend.title = element_blank(), 
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.background = element_rect(fill = "white", colour = "grey50"))+
  facet_wrap(~variable, scales = "free_y", nrow=1, labeller=as_labeller(traits,label_parsed))


###GLMM for leaf traits and Gypsophilic values (CWM and FD), Simpson, Vegetation cover and canopy height
#CWM
data_CWM
#FD
data_Rao
#GV
CWM_GV
#Simpson
Simpson
#Vegetation cover and canopy height
varechem

glm<-lmer(N~grazing+(1|locality)+(1|locality:plot), data=data_Rao)
res<-residuals(glm)
a<-simulateResiduals(glm)
plot(a)


glm<-glmer(S~grazing+(1|locality)+(1|locality:plot), data=data_Rao, family=Gamma)
glm<-glmer.nb(b.simpson~grazing+(1|locality)+(1|locality:plot), data=Simpson)
glm<-glmmTMB(canopy_height~grazing+(1|locality)+(1|locality:plot),family=Gamma(link="identity"),data=varechem)
AIC(glm)

varechem$grazing<-as.factor(varechem$grazing)
Anova(glm,test="Chi")
c<-glht(glm,linfct = mcp(grazing = "Tukey"))
summary(c, test = adjusted("bonferroni"))
tuk.cld <- cld(c)
plot(tuk.cld)


###Correction type I error
varespec <- read_excel("data_paper_AndreuCera.xlsx", sheet="vegetation_data")
GV<-read_excel("data_paper_AndreuCera.xlsx", sheet="gypsophilic_value")
varespec$grazing<-as.factor(varespec$grazing)
CWM_traits <- cwm (com = varespec[,5:56], traits = GV[,2])
re_traits <- test_cwm (cwm = CWM_traits, env= varespec[,2], method = 'aov', adjustP = TRUE, test="all", perm=999)
re_traits
plot(re_traits)



###Disentangling effects of ITV and turnover species following Leps et al. 2011's method.
CWM_fixed_leaf
CWM_Specific_leaf
CWM_ITV_leaf

Rao_fixed_leaf
Rao_Specific_leaf
Rao_ITV_leaf


adonis(S~pert,data=CWM_fixed_leaf, method="euclidean")

#data saved on excel to calculate as Leps et al. 2011 proposed


##graph Figure 3
data <- read_excel("data_paper_AndreuCera.xlsx", sheet="data_Leps")
data<-data.frame(data)
s<-ggplot(data=data, aes(x=trait, y=value, fill=factor(variable)))
s+geom_bar(stat = "identity") + labs(y="% explained variability", x=element_blank(), size=15)+
  scale_x_discrete(limits=c("S","N","C"), labels=c("Leaf S","Leaf N", "Leaf C"))+
  geom_boxplot(data=data,aes(x=trait, y=spec))+
  theme(
    legend.position="none",
    strip.text.x = element_text(size=15),
    strip.background = element_rect(colour="black", fill="white"),
    axis.text=element_text(size=15, colour="black"),
    axis.title.x = element_text(size=15),
    legend.title = element_blank(), 
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.background = element_rect(fill = "white", colour = "grey50"))+
  scale_fill_manual(values=c('gray91','gray55'))+
  facet_wrap(~indices)+ 
  coord_flip()


###Figure S3 Differences gypsophilic value between localities
CWM_GV
tgc1 <- summarySE(na.omit(CWM_GV), measurevar="Gypsophilic_value", groupvars=c("locality","grazing"))

s<-ggplot(data=tgc1, aes(x=grazing, y=Gypsophilic_value, fill=grazing))
s+geom_bar(stat="identity", colour="black", position=position_dodge()) +
  scale_x_discrete(limits=c("Low","Medium","High"))+
  labs(y="Gypsophilic value", x=element_blank())+
  scale_fill_grey()+
  geom_errorbar(aes(ymin=Gypsophilic_value-se, ymax=Gypsophilic_value+se), width=.5,position=position_dodge(.9))+
  theme(
    legend.position="none",
    strip.text.x = element_text(size=15, face="bold"),
    strip.background = element_rect(colour="black", fill="white"),
    axis.text=element_text(size=15),
    axis.title.y = element_text(size = 20, angle = 90),
    legend.title = element_blank(), 
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.background = element_rect(fill = "white", colour = "grey50"))+
  facet_wrap(~locality, nrow=1)
