rm(list=ls(all=TRUE))
dir()

## Author: Rony Peterson Santos Almeida (PPGZOOL-UFPA-Brazil)
## Email: rony__peterson@hotmail.com
## Date: December 15, 2020

### Packages:
library(fitdistrplus)
library(gamlss)
library(vegan)
library(ggplot2)
library(ggrepel)
library(gridExtra)
library(grid)

### Colors and graphic parameters
cor_ESECAFLOR <- c("#0072B2", "#E69F00")
cor_pont <- "red"
n_alpha_jit <- 0.5
n_point_jit <- 2.8
col_jit <- 1

### Functions:
#Calculate mean and deviation:
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,conf.interval=.95, .drop=TRUE) {
  library(plyr)
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)}
  datac <- ddply(data, groupvars, .drop=.drop,.fun = function(xx, col) {c(N    = length2(xx[[col]], na.rm=na.rm), mean = mean   (xx[[col]], na.rm=na.rm),sd = sd(xx[[col]], na.rm=na.rm))},measurevar)
  datac <- rename(datac, c("mean" = measurevar))
  datac$se <- datac$sd / sqrt(datac$N)
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  return(datac)
}

#fitidist() for each species
func_fitdist <- function(x){
  fitL <- fitdist(x, "lnorm")
  fitG <- fitdist(x, "gamma")
  fitW <- fitdist(x, "weibull")
  objeto <- list(fitL, fitG, fitW)
  legend_fit <- c("Lognorm", "Gamma", "Weibull")
  aic_family <- gofstat(objeto, fitnames=legend_fit)
  aic_results <- data.frame(t(aic_family$aic), model=ssp_name)
  return(aic_results)
}


#### 1 - Detection distance----
#### 1.1 - Transect----
tab1 <- read.csv("comunity_transect.csv", h=T, sep = ";", stringsAsFactors = T)
str(tab1)
names(tab1)

#Selection of the distribution family
variable <- tab1$Distance
fitL <- fitdist(variable, "lnorm"); fitG <- fitdist(variable, "gamma"); fitW <- fitdist(variable, "weibull")

plot(fitL)
plot(fitG)
plot(fitW)

objeto <- list(fitL,fitG, fitW)
legend_fit <- c("Lognorm", "Gamma", "Weibull")

cdfcomp(objeto, legendtext=legend_fit)
denscomp(objeto, legendtext=legend_fit)
qqcomp(objeto, legendtext=legend_fit)
ppcomp(objeto, legendtext=legend_fit)

gofstat(objeto, fitnames=legend_fit)

(aic_family <- gofstat(objeto, fitnames=legend_fit))
tab_aic_family <- data.frame(t(aic_family$aic), model="Especie")

#Model:
m1 <- gamlss(Distance ~ Plot, data=tab1, family = WEI())
plot(m1)
var <- t(data.frame(summary(m1)[2,]))

tab_modelos <- data.frame(var, DF="1, 18", var="dist_median", level="transect", family = "Weibull")


### Figure:
tab_simp1 <- summarySE(tab1, measurevar="Distance", groupvars= "Plot")

f1 <- ggplot(data=tab1, aes(x=Plot , y=Distance, fill=Plot))+
  theme_classic()+
  geom_boxplot(width=0.75, size = 0.8, alpha=0.4, position=position_dodge(0.5)) +
  geom_jitter(position=position_jitter(0.3), alpha=n_alpha_jit, size=2, color=cor_pont)+
  scale_y_continuous(limits=c(0, 150), breaks = c(0,50,100,150))+
  scale_fill_manual(values=cor_ESECAFLOR)+
  labs(x = "Plot", y = "Distance median (Transect; cm)", size = 3)+
  geom_text(aes(x=1.5, y=150, label ="NS"), size=4)+
  theme(legend.position = "none",
        axis.line.x.bottom = element_line(color = "black", size = .8),
        axis.text.x.bottom = element_text(color = "black", size = 15),
        axis.line.y.left = element_line(color = "black", size = .8),
        axis.text.y.left = element_text(color = "black"),
        axis.title.y = element_text(color = "black", size = 15))



#### 1.2 - Species----
tab2 <- read.csv("species_plot_CLAM.csv", h=T, sep = ";", stringsAsFactors = T)
str(tab2)
names(tab2)

tab2.1 <- tab2[1:12,] #Espécies com mais de 4 registros por Plot

hist(tab2.1$Distance)
shapiro.test(tab2.1$Distance)
t.test(Distance ~ Plot, data = tab2.1, paired = TRUE)

mod <- t.test(Distance ~ Plot, data = tab2.1, paired = TRUE)

values <- data.frame(NA, NA, mod$statistic, mod$p.value, DF= "6, 11", var="dist_median", level="species", family = "Gausian")
colnames(values) <- colnames(tab_modelos)
tab_modelos <- rbind(tab_modelos, values)


## Figure
con_cor <- tab2.1[which(tab2.1$Plot=="Control"),]
exp_cor <- tab2.1[which(tab2.1$Plot=="Experiment"),]
data2 <- data.frame(rbind(con_cor, exp_cor))
data2$x_cor <- c(rep(1, dim(con_cor)[1]), rep(2, dim(exp_cor)[1]))

data3 <- data.frame(con_cor,exp_cor)
data3$x_s <- rep(1, dim(con_cor)[1])
data3$x_v <- rep(2, dim(exp_cor)[1])

ssp <- ggplot(data2, aes(x=Plot, y=Distance, fill=Plot)) +
  geom_boxplot(width=0.5, size = 0.8, alpha=0.4, position=position_dodge(0.25))+
  geom_segment(data = data3, aes(x = x_s, xend = x_v, y = Distance, yend = Distance.1), colour = "gray80", size = 0.6)+
  geom_text_repel(data = exp_cor, aes(y = Distance, x = 2, label = Specie), nudge_x = 0.6, direction = "y", hjust = 1, segment.size = 0.2, fontface="italic", segment.color = "gray50", colour="gray50", size = 3) +
  geom_point(aes(x = x_cor, y = Distance), colour = "black", size = 2)+
  geom_text(aes(x=1.5, y=90, label ="NS"), size=4)+
  labs(y = "Distance median (Species; cm)")+
  scale_fill_manual(values=cor_ESECAFLOR)+
  scale_y_continuous(limits=c(18, 90), breaks = c(18,42,66,90))+
  theme_classic()+
  theme(legend.position = "none",
        axis.line.x.bottom = element_line(color = "black", size = .8),
        axis.text.x.bottom = element_text(color = "black", size = 15),
        axis.line.y.left = element_line(color = "black", size = .8),
        axis.text.y.left = element_text(color = "black"),
        axis.title.y = element_text(color = "black", size = 15))



#### 1.3 - Habitat association----
#### _1.3.1 - CLAM----
names(tab1) 
com <- tab1[,7:ncol(tab1)] #Community

result_CLAM <- clamtest(com, tab1$Plot, alpha = 0.05, specialization = 0.5)
summary(result_CLAM)


## Figure
dados <- transform(result_CLAM,
                 Total_Controle=1+(get("Total_Control")),
                 Total_Experimental=1+(get("Total_Experiment")))
minval<-summary(result_CLAM)$minv
lineYvsgen<-minval[[1]]+1
lineXvsgen<-minval[[2]]+1

Ymin<-minval[[1]][1,2]+1; Xmin<-minval[[2]][1,1]+1

esp_s <- rbind(dados[which(dados$Classes == "Specialist_Control"), ],
               dados[which(dados$Classes == "Specialist_Experiment"), ],
               dados[which(dados$Classes == "Generalist"), ])

#Colors
cor_CLAM <- c(1,cor_ESECAFLOR,"gray60")

graphClam <-
  ggplot(dados, aes(x=Total_Experimental, y=Total_Controle, colour=Classes))+
  geom_curve(aes(x=Xmin, y=1, xend=1, yend=Ymin), size=1.5, colour=cor_CLAM[4], alpha=0.2, curvature = .55)+
  geom_line(data=lineYvsgen, aes(x=x, y=y), linetype=2, size=1.5, colour=cor_CLAM[2], alpha=0.5)+
  geom_line(data=lineXvsgen, aes(x=x, y=y), linetype=2, size=1.5, colour=cor_CLAM[3], alpha=0.5)+
  geom_text_repel(data=esp_s, aes(label=gsub("_"," ",esp_s[,1])), show.legend=F, fontface="italic", colour="gray30", direction="y", hjust=1.3, segment.size = 0.3, size=4)+
  geom_point(aes(colour=Classes, shape=Classes), size=4)+
  scale_shape_manual(values=c(16,18,18,17), breaks=c("Generalist", "Specialist_Control", "Specialist_Experiment", "Too_rare"), labels=c("No preference", "Control", "Experiment", "Too rare"))+
  scale_colour_manual(values = cor_CLAM, breaks=c("Generalist", "Specialist_Control", "Specialist_Experiment", "Too_rare"), labels=c("No preference", "Control", "Experiment", "Too rare"))+
  scale_x_continuous(trans = 'log10', breaks = c(1,4,8,16)) +
  scale_y_continuous(trans = 'log10', breaks = c(1,5,10,20))+
  labs(x="Experiment (abundance + 1)", y="Control (abundance + 1)")+
  theme_classic()+
  theme(legend.position=c(0.1, 0.1),
        legend.title=element_blank(),
        legend.background = element_rect(color = "black", size = .5, linetype = "solid"),
        legend.key.size = unit(0, 'lines'))



#### _1.3.2 - Model Species CLAM----
names(tab2)

tab2.2 <- tab2[13:25,c(2,4,5)] #Espécies classificadas pelo CLAM
as.factor(esp_s$Species[order(esp_s$Species)]) == as.factor(tab2.2$Specie[order(tab2.2$Specie)])

variable <- tab2.2$Distance

fitL <- fitdist(variable, "lnorm")
fitG <- fitdist(variable, "gamma")
fitW <- fitdist(variable, "weibull")

objeto <- list(fitL, fitG, fitW)
legend_fit <- c("Lognorm", "Gamma", "Weibull")


(aic_family <- gofstat(objeto, fitnames=legend_fit))
aic_family2 <- data.frame(t(aic_family$aic), model="ssp_Classes")
tab_aic_family <- rbind(tab_aic_family, aic_family2)

m4 <- gamlss(Distance ~ Group , data = tab2.2, family = LOGNO())
var <- t(data.frame(summary(m4)[4,]))

values <- data.frame(var, DF="4,13", var="dist_median", level="species_CLAM", family = "Lognorm")

tab_modelos <- rbind(tab_modelos, values)

## Pairwise Test:
tab2.2 <- tab2.2[order(tab2.2$Group),]
summary(glm(Distance ~ Group, data = tab2.2[1:8,])) #Control x Experiment
summary(glm(Distance ~ Group, data = tab2.2[c(1:3, 9:13),])) #Control x Generalist
summary(glm(Distance ~ Group, data = tab2.2[4:13,])) #Experiment X Generalist


## Figure:
laters <- aggregate(tab2.2$Distance, list(tab2.2$Group), max)
laters$x2 <- laters$x+15
laters$let <- c("a","a","b")
colnames(laters)[1] <- "Group"

p_Classes <-
  ggplot(data=tab2.2, aes(x=Group, y=Distance, fill=Group))+
  geom_boxplot(width=0.5, size = 0.8, alpha=0.4, position=position_dodge(0.25))+
  geom_text_repel(aes(label = Specie), nudge_x = 0.55, direction = "y", hjust = 1, segment.size = 0.2, fontface="italic", segment.color = "gray50", colour="gray50", size = 3) +
  geom_point(colour = "black", size = 2)+
  geom_text(aes(x=2, y=220, label ="**"), size=7)+
  geom_text(data=laters, aes(x=Group, y=x2, label=let), size=6, colour="gray50")+
  scale_x_discrete(breaks=c("Control_Specialist", "Experiment_Specialist",  "Generalist"), labels=c("Control\nSpecialist", "Experimen\nSpecialist", "Generalist"), limits=c("Control_Specialist", "Generalist", "Experiment_Specialist"))+
  scale_y_continuous(limits=c(0, 220), breaks = c(0,110,220))+
  scale_fill_manual(values=c(cor_ESECAFLOR[1], cor_ESECAFLOR[2], 1))+
  labs(x = "Groups", y = "Distance median (Species; cm)", size = 3)+
  theme_classic()+
  theme(legend.position = "none",
        axis.line.x.bottom = element_line(color = "black", size = .8),
        axis.text.x.bottom = element_text(color = "black", size = 15),
        axis.line.y.left = element_line(color = "black", size = .8),
        axis.text.y.left = element_text(color = "black"),
        axis.title.y = element_text(color = "black", size = 15))



#### 1.4 - Specie----
tab3 <- read.csv("species_colony.csv", h=T, sep = ";", stringsAsFactors = T)
str(tab3)
levels(tab3$Specie)


## Each species separately:
ssp_name <- c("Cr.lev")
BD <- tab3[which(tab3$Specie==ssp_name),]
variable <- BD$Distance

(aic_family2 <- func_fitdist(variable))
tab_aic_family <- rbind(tab_aic_family, aic_family2) #"Lognorm"

ms1 <- gamlss(BD$Distance ~ BD$Plot, family = LOGNO())
var <- t(data.frame(summary(ms1)[2,]))
values <- data.frame(var, DF="2,19", var="dist_colony", level=ssp_name, family = "Lognorm")
tab_modelos <- rbind(tab_modelos, values)

ssp1 <- ggplot(BD, aes(x=Plot, y=Distance, fill=Plot)) +
  geom_boxplot(width=0.5, size = 0.8, alpha=0.4) +
  geom_jitter(position=position_jitter(0.2), alpha=n_alpha_jit, size=2, color=col_jit)+
  geom_text(aes(x=1.5, y=400, label ="NS"), size=4)+
  theme_classic()+
  labs(x = "Plot", y = "Distance (cm)", title = "Cr.lev")+
  scale_fill_manual(values=cor_ESECAFLOR)+
  theme(legend.position = "none",
        plot.title = element_text(face="italic"),
        axis.line.x.bottom = element_line(color = "black", size = .8),
        axis.text.x.bottom = element_text(color = "black", size = 15),
        axis.line.y.left = element_line(color = "black", size = .8),
        axis.text.y.left = element_text(color = "black"),
        axis.title.y = element_text(color = "black", size = 15))



ssp_name <- c("Pa.bug")
BD <- tab3[which(tab3$Specie==ssp_name),]
variable <- BD$Distance

(aic_family2 <- func_fitdist(variable))
tab_aic_family <- rbind(tab_aic_family, aic_family2) #"Lognorm"

ms1 <- gamlss(BD$Distance ~ BD$Plot, family = LOGNO())
var <- t(data.frame(summary(ms1)[2,]))
values <- data.frame(var, DF="2,12", var="dist_colony", level=ssp_name, family = "Lognorm")
tab_modelos <- rbind(tab_modelos, values)

ssp2 <- ggplot(BD, aes(x=Plot, y=Distance, fill=Plot)) +
  geom_boxplot(width=0.5, size = 0.8, alpha=0.4) +
  geom_jitter(position=position_jitter(0.2), alpha=n_alpha_jit, size=2, color=col_jit)+
  geom_text(aes(x=1.5, y=100, label ="NS"), size=4)+
  theme_classic()+
  labs(x = "Sites", y = "Distance (cm)", title = "Pa.bug")+
  scale_fill_manual(values=cor_ESECAFLOR)+
  theme(legend.position = "none",
        plot.title = element_text(face="italic"),
        axis.line.x.bottom = element_line(color = "black", size = .8),
        axis.text.x.bottom = element_text(color = "black", size = 15),
        axis.line.y.left = element_line(color = "black", size = .8),
        axis.text.y.left = element_text(color = "black"),
        axis.title.y = element_text(color = "black", size = 15))



ssp_name <- c("Ph.cur")
BD <- tab3[which(tab3$Specie==ssp_name),]
variable <- BD$Distance

(aic_family2 <- func_fitdist(variable))
tab_aic_family <- rbind(tab_aic_family, aic_family2) #"Weibull"

ms1 <- gamlss(BD$Distance ~ BD$Plot, family = WEI())
var <- t(data.frame(summary(ms1)[2,]))
values <- data.frame(var, DF="2,18", var="dist_colony", level=ssp_name, family = "Weibull")
tab_modelos <- rbind(tab_modelos, values)

ssp3 <- ggplot(BD, aes(x=Plot, y=Distance, fill=Plot)) +
  geom_boxplot(width=0.5, size = 0.8, alpha=0.4) +
  geom_jitter(position=position_jitter(0.2), alpha=n_alpha_jit, size=2, color=col_jit)+
  geom_text(aes(x=1.5, y=160, label ="NS"), size=4)+
  theme_classic()+
  labs(x = "Sites", y = "Distance (cm)", title = "Ph.cur")+
  scale_fill_manual(values=cor_ESECAFLOR)+
  theme(legend.position = "none",
        plot.title = element_text(face="italic"),
        axis.line.x.bottom = element_line(color = "black", size = .8),
        axis.text.x.bottom = element_text(color = "black", size = 15),
        axis.line.y.left = element_line(color = "black", size = .8),
        axis.text.y.left = element_text(color = "black"),
        axis.title.y = element_text(color = "black", size = 15))



ssp_name <- c("Ph.fra")
BD <- tab3[which(tab3$Specie==ssp_name),]
variable <- BD$Distance

(aic_family2 <- func_fitdist(variable))
tab_aic_family <- rbind(tab_aic_family, aic_family2) #"Lognorm"

ms1 <- gamlss(BD$Distance ~ BD$Plot, family = LOGNO())
var <- t(data.frame(summary(ms1)[2,]))
values <- data.frame(var, DF="2,28", var="dist_colony", level=ssp_name, family = "Lognorm")
tab_modelos <- rbind(tab_modelos, values)

ssp4 <- ggplot(BD, aes(x=Plot, y=Distance, fill=Plot)) +
  geom_boxplot(width=0.5, size = 0.8, alpha=0.4) +
  geom_jitter(position=position_jitter(0.2), alpha=n_alpha_jit, size=2, color=col_jit)+
  geom_text(aes(x=1.5, y=100, label ="NS"), size=4)+
  theme_classic()+
  labs(x = "Sites", y = "Distance (cm)", title = "Ph.fra")+
  scale_fill_manual(values=cor_ESECAFLOR)+
  theme(legend.position = "none",
        plot.title = element_text(face="italic"),
        axis.line.x.bottom = element_line(color = "black", size = .8),
        axis.text.x.bottom = element_text(color = "black", size = 15),
        axis.line.y.left = element_line(color = "black", size = .8),
        axis.text.y.left = element_text(color = "black"),
        axis.title.y = element_text(color = "black", size = 15))



ssp_name <- c("Ph.sch")
BD <- tab3[which(tab3$Specie==ssp_name),]
variable <- BD$Distance

(aic_family2 <- func_fitdist(variable))
tab_aic_family <- rbind(tab_aic_family, aic_family2) #"Weibull"

ms1 <- gamlss(BD$Distance ~ BD$Plot, family = WEI())
var <- t(data.frame(summary(ms1)[2,]))
values <- data.frame(var, DF="2,7", var="dist_colony", level=ssp_name, family = "Weibull")
tab_modelos <- rbind(tab_modelos, values)

ssp5 <- ggplot(BD, aes(x=Plot, y=Distance, fill=Plot)) +
  geom_boxplot(width=0.5, size = 0.8, alpha=0.4) +
  geom_jitter(position=position_jitter(0.2), alpha=n_alpha_jit, size=2, color=col_jit)+
  geom_text(aes(x=1.5, y=80, label ="NS"), size=4)+
  theme_classic()+
  labs(x = "Sites", y = "Distance (cm)", title = "Ph.sch")+
  scale_fill_manual(values=cor_ESECAFLOR)+
  theme(legend.position = "none",
        plot.title = element_text(face="italic"),
        axis.line.x.bottom = element_line(color = "black", size = .8),
        axis.text.x.bottom = element_text(color = "black", size = 15),
        axis.line.y.left = element_line(color = "black", size = .8),
        axis.text.y.left = element_text(color = "black"),
        axis.title.y = element_text(color = "black", size = 15))



ssp_name <- c("Ph.sco")
BD <- tab3[which(tab3$Specie==ssp_name),]
variable <- BD$Distance

(aic_family2 <- func_fitdist(variable))
tab_aic_family <- rbind(tab_aic_family, aic_family2) #"Weibull"

ms1 <- gamlss(BD$Distance ~ BD$Plot, family = WEI())
var <- t(data.frame(summary(ms1)[2,]))
values <- data.frame(var, DF="2,7", var="dist_colony", level=ssp_name, family = "Weibull")
tab_modelos <- rbind(tab_modelos, values)

ssp6 <- ggplot(BD, aes(x=Plot, y=Distance, fill=Plot)) +
  geom_boxplot(width=0.5, size = 0.8, alpha=0.4) +
  geom_jitter(position=position_jitter(0.2), alpha=n_alpha_jit, size=2, color=col_jit)+
  geom_text(aes(x=1.5, y=41, label ="NS"), size=4)+
  theme_classic()+
  labs(x = "Sites", y = "Distance (cm)", title = "Ph.sco")+
  scale_fill_manual(values=cor_ESECAFLOR)+
  theme(legend.position = "none",
        plot.title = element_text(face="italic"),
        axis.line.x.bottom = element_line(color = "black", size = .8),
        axis.text.x.bottom = element_text(color = "black", size = 15),
        axis.line.y.left = element_line(color = "black", size = .8),
        axis.text.y.left = element_text(color = "black"),
        axis.title.y = element_text(color = "black", size = 15))



#### 2 - Colony number, richness ans composition----
#### 2.1 - Colony number----
names(tab1)
variable <- tab1$Colony
fitP <- fitdist(variable, "pois")
fitN <- fitdist(variable, "nbinom")

objeto <- list(fitP, fitN)
legend_fit <- c("Poisson", "Neg.Binom.")

(aic_family <- gofstat(objeto, fitnames=legend_fit)) #Neg.Binom.
aic_family2 <- data.frame(t(aic_family$aic), NA, model="Number_colonys")
aic_family2 <- rbind(colnames(aic_family2),aic_family2)
colnames(aic_family2) <- colnames(tab_aic_family)
tab_aic_family <- rbind(tab_aic_family, aic_family2)

m3 <- gamlss(Colony ~ Plot, data=tab1, family = NBI())
var <- t(data.frame(summary(m3)[2,]))

values <- data.frame(var, DF="2,17", var="number_colonys", level="transect", family = "Neg.Binom.")

tab_modelos <- rbind(tab_modelos, values)

## Figure:
tab_simp3 <- summarySE(tab1, measurevar="Colony", groupvars=c("Plot"))

f3 <- ggplot(data=tab_simp3, aes(x=Plot, y=Colony, fill=Plot))+
  theme_classic()+
  geom_bar(stat="identity", color="gray20", size = 0.8, alpha=0.4)+
  geom_jitter(data=tab1, aes(x=Plot, y=Colony), height = 0.0001, width = 0.35, alpha=n_alpha_jit, size=n_point_jit, color=cor_pont)+
  geom_errorbar(aes(ymin=Colony-ci, ymax=Colony+ci), width=.2, position=position_dodge(.9), size=.9)+
  scale_fill_manual(values=cor_ESECAFLOR)+
  labs(x = "Plot", y = "Colony number (mean±IC)", size = 3)+
  geom_text(aes(x=1.5, y=30, label ="***"), size=7)+
  theme(legend.position = "none",
        axis.line.x.bottom = element_line(color = "black", size = .8),
        axis.text.x.bottom = element_text(color = "black", size = 15),
        axis.line.y.left = element_line(color = "black", size = .8),
        axis.text.y.left = element_text(color = "black"),
        axis.title.y = element_text(color = "black", size = 15))



#### 2.2 - Richness----
variable <- tab1$Richness
fitP <- fitdist(variable, "pois")
fitN <- fitdist(variable, "nbinom")

objeto <- list(fitP, fitN)
legend_fit <- c("Poisson", "Neg.Binom.")

(aic_family <- gofstat(objeto, fitnames=legend_fit)) #Poisson
aic_family2 <- data.frame(t(aic_family$aic), NA, model="Richnes")
colnames(aic_family2) <- colnames(tab_aic_family)
tab_aic_family <- rbind(tab_aic_family, aic_family2)

m4 <- gamlss(Richness ~ Plot, data=tab1, family = PO())
var <- t(data.frame(summary(m4)[2,]))

values <- data.frame(var, DF="2,17", var="richnes", level="transect", family = "Poisson")

tab_modelos <- rbind(tab_modelos, values)
tab_modelos <- data.frame(round(tab_modelos[,1:4],4), tab_modelos[,5:8])


## Figure:
tab_simp4 <- summarySE(tab1, measurevar="Richness", groupvars=c("Plot"))

f4 <- ggplot(data=tab_simp4, aes(x=Plot, y=Richness, fill=Plot))+
  theme_classic()+
  geom_bar(stat="identity", color="gray20", size = 0.8, alpha=0.4)+
  geom_jitter(data=tab1, aes(x=Plot, y=Richness), height = 0.0001, width = 0.35, alpha=n_alpha_jit, size=n_point_jit, color=cor_pont)+
  geom_errorbar(aes(ymin=Richness-ci, ymax=Richness+ci), width=.2, position=position_dodge(.9), size=.9)+
  scale_fill_manual(values=cor_ESECAFLOR)+
  labs(x = "Sites", y = "Ants richnes (mean±IC)", size = 3)+
  geom_text(aes(x=1.5, y=20, label ="***"), size=7)+
  theme(legend.position = "none",
        axis.line.x.bottom = element_line(color = "black", size = .8),
        axis.text.x.bottom = element_text(color = "black", size = 15),
        axis.line.y.left = element_line(color = "black", size = .8),
        axis.text.y.left = element_text(color = "black"),
        axis.title.y = element_text(color = "black", size = 15))



#### 2.3 - Composition----
#### _2.3.1 - PERMANOVA----
result <- adonis(com ~ tab1$Plot ,permutations=9999, method = "bray")
result


#### _2.3.2 - PCoA----
ssp_dist <- vegdist(com, method="bray")

pcoa <- cmdscale(ssp_dist, k=2, eig=T)
autovalor <- (pcoa$eig) # retira o autovalor da análise
explicacao <- round(((pcoa$eig/sum(pcoa$eig))*100),4) # % de explicação de cada eixo.

#Species
especie <- wascores(pcoa$points[,1:2], com)
colnames(especie) <- c("PCC1","PCC2")
especie2 <- data.frame(especie, names_ssp=row.names(especie))

#Points
points1 <- data.frame(pcoa$points[,1:2])#Autovetor
points <- data.frame(points1, tab1$Plot)
colnames(points) <- c("PC1","PC2", "Plot")

#Hull
data_con <- points[1:10,]
hull_con <- data_con[chull(data_con$PC1, data_con$PC2),]
data_exp <- points[11:20,]
hull_exp <- data_exp[chull(data_exp$PC1, data_exp$PC2),]
hull <- rbind(hull_con, hull_exp)
hull <- data.frame(hull, Plot=c(rep("Control", nrow(hull_con)), rep("Experiment", nrow(hull_exp))))


p_pcoa1 <- ggplot(data=points, aes(x = PC1, y = PC2))+
  geom_hline(yintercept = 0, colour="gray", linetype="dashed")+geom_vline(xintercept = 0, colour="gray", linetype="dashed")+
  geom_polygon(data = hull, aes(x=PC1, y=PC2, fill=Plot), colour="gray20", size =.5, alpha=.2)+
  geom_point(aes(shape=Plot, colour=Plot), size = 4)+
  geom_text_repel(data = especie2, aes(x = PCC1, y = PCC2, label = names_ssp), size=2.4, show.legend=FALSE, fontface="italic", colour="gray50")+
  geom_text(aes(x = .05, y = .35, label = "***"), size = 7)+
  scale_color_manual(values = cor_ESECAFLOR)+
  scale_fill_manual(values = cor_ESECAFLOR)+
  scale_shape_manual(values = c(15, 16))+
  labs(x = "PCoA 1 (24.034%)", y = "PCoA 2 (12.202%)", size = 2)+
  theme_classic()+
  theme(legend.position = c(0.9, 0.9),
        legend.background = element_rect(color = "black", size = .5, linetype = "solid"),
        panel.border = element_rect(colour = "black", fill=NA, size=1.2),
        axis.text.x.bottom = element_text(color = "black"),
        axis.text.y.left = element_text(color = "black"))



##3 - Save results and Figures----

write.csv(tab_modelos,"RES_Models.csv")
write.csv(tab_aic_family,"RES_AIC_family.csv")

lay1 <- rbind(c(0,0,0,1,1,1,1),
              c(2,2,2,2,2,2,2))

png(file = "Figure2.png", width = 3000, height = 4000, units = "px", res = 400)
grid.arrange(arrangeGrob(f1 +ggtitle("(a) - Transect")+labs(y = NULL)+ theme(plot.margin = unit(c(.1, .5, 0.8, 0), "cm")), #topo, direita, abaixo, esquerda
                         ssp +ggtitle("(b) - Species")+labs(y = NULL)+ theme(plot.margin = unit(c(.1, .5, 0.8, 0), "cm")),
                         p_Classes +ggtitle("(c) - Habitat association")+labs(y = NULL)+ theme(plot.margin = unit(c(.1, .5, 0.2, 0), "cm")),
                         layout_matrix = lay1, left = textGrob("Distance median (cm)", rot = 90, vjust = 0.3, gp=gpar(fontsize=15))))
dev.off()


png(file = "Figure3.png", width = 2800, height = 3600, units = "px", res = 400)
grid.arrange(arrangeGrob(ssp1+labs(x = NULL, y = NULL)+scale_x_discrete(labels=NULL),ssp2+labs(x = NULL, y = NULL)+scale_x_discrete(labels=NULL), ssp3+labs(x = NULL, y = NULL)+scale_x_discrete(labels=NULL), ssp4+labs(x = NULL, y = NULL)+scale_x_discrete(labels=NULL), ssp5+labs(x = "", y = NULL), ssp6+labs(x = "", y = NULL),
                         ncol=2, heights = c(4.5, 4.5, 5), widths = c(1.1,1),
                         bottom = textGrob("Plot", vjust = -.5 , gp=gpar(fontsize=15)),
                         left = textGrob("Distance (cm)", rot = 90, vjust = 1, gp=gpar(fontsize=15))))
dev.off()


png(file = "Figure4.png", width = 2200, height = 3500, units = "px", res = 400)
grid.arrange(arrangeGrob(f3+ggtitle("(a)")+labs(x = NULL)+scale_x_discrete(labels=NULL),
                         f4+ggtitle("(b)")+labs(x = NULL),
                         ncol=1, heights = c(5.5,6),
                         bottom = textGrob("Plot", vjust = -.5 , gp=gpar(fontsize=15))))
dev.off()


png(file = "Figure5.png", width = 2200, height = 1800, units = "px", res = 300)
p_pcoa1
dev.off()


png(file = "FigureS2.png", width = 2000, height = 1500, units = "px", res = 300)
graphClam
dev.off()


#### CITATION----
citation()
citation("fitdistrplus")
citation("gamlss")
citation("vegan")
citation("ggplot2")
citation("ggrepel")
citation("gridExtra")