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

##Choose color:
#This color will be used for the plots (Control and Experiment, respectively)
cor_ESECAFLOR <- c("#0072B2", "#E69F00") #Colors
#or
#cor_ESECAFLOR <- c("black", "gray60") #B&W

##Functions----
#Functions for correlation figures

#Continuous figure function:
my_fn <- function(data, mapping, ...){
  p <- ggplot(data = data, mapping = mapping) + 
    geom_point(aes(colour=Plot, shape=Plot), size=2) +
    geom_smooth(method=lm, formula = y ~ x, color = 1, se = F)+
    scale_color_manual(values = cor_ESECAFLOR)+
    scale_shape_manual(values = c(16, 17))
  p
}

#Diagonal figure function:
my_fn2 <- function(data, mapping, ...){
  p <- ggplot(data = data, mapping = mapping) + 
    geom_histogram(bins = 10, colour="gray20")
  p
}

#Correlation figure function:
cor_fun <- function(data, mapping, method="pearson", ndp=2, sz=5, stars=FALSE, ...){
  x <- eval_data_col(data, mapping$x)
  y <- eval_data_col(data, mapping$y)
  
  corr <- cor.test(x, y, method=method)
  est <- corr$estimate
  lb.size <- sz
  
  if(stars){
    stars <- c("***", "**", "*", "")[findInterval(corr$p.value, c(0, 0.001, 0.01, 0.05, 1))]
    lbl <- paste0(round(est, ndp), stars)
  }else{
    lbl <- round(est, ndp)
  }
  ggplot(data=data, mapping=mapping) + 
    annotate("text", x=mean(x, na.rm=TRUE), y=mean(y, na.rm=TRUE), label=lbl, size=lb.size,...)+
    theme(panel.grid = element_blank())
}


#####################
##1 - Analysis months----

tab_m <- read.table("ants_months_ESECAFLOR_2011-2012.txt",h=T,sep="\t")

##_Correlation of variables ----
var_pairs <- tab_m[,c(4,5,6,1)]

library(GGally)
library(ggplot2)

p1 <- ggpairs(var_pairs, columns = 1:3,
              lower = list(continuous = my_fn),
              diag = list(continuous = my_fn2), 
              upper = list(continuous = cor_fun))

cor_m <- p1+theme_minimal()+
  theme(legend.position = "none", 
        panel.grid.minor = element_blank(),
        axis.line.x.bottom = element_line(color = "black", size = .8),
        axis.text.x.bottom = element_text(color = "black"),
        axis.line.y.left = element_line(color = "black", size = .8),
        axis.text.y.left = element_text(color = "black"),
        panel.border = element_rect(linetype = "solid", colour = "gray50", fill = NA),
        strip.text.x = element_text(size=12),
        strip.text.y = element_text(size=12),
        strip.background = element_rect(colour="black", fill="white"))


##_1.1 - Richness ----
##_1.1.1 - Best model and selection AIC ----
library(nlme)
mod1 <- nlme::lme(Sobs ~ (Biomass + Litt.Div. + Moisture) * Plot, random=~1|Plot, method = "REML", data=tab_m)

options(na.action = "na.fail")
library(MuMIn)
selec <- dredge(mod1, rank="AICc", extra = "adjR^2")
(mod_acaic <-summary(model.avg(selec))) #Resultado do Akaik
importance(selec) #Substituto para a importância das variáveis no modelo quando há interação. Refazer usando isso em detrimento do hier.part


##_1.1.2 - Importance of variables (Plot)----
plot(importance(selec)); abline(h=0.7, col=2)

import_var_a <- data.frame(importance(selec))
import_var_a$nomes_var_a <- row.names(import_var_a)
import_var_a$import_a <- round(import_var_a$importance.selec., 3)
import_var_a2 <- import_var_a[order(import_var_a$import_a, decreasing = T),]

imp_m <-
ggplot(data=import_var_a2, aes(x=reorder(nomes_var_a, +import_a), y=import_a)) +
  geom_bar(position="dodge",stat="identity", colour="black") + 
  geom_hline(aes(yintercept=0.7), color="red", linetype="dashed", size=1)+
  scale_y_continuous(breaks = c(0.1,0.4,0.7,1))+
  labs(x = "Temporal variables", y = "", size = 20)+
  theme_classic()+
  coord_flip()+
  theme(axis.line = element_line (size=.4),
        axis.line.x.bottom = element_line(color = "black", size = .6),
        axis.text.x.bottom = element_text(color = "black", size = 12),
        axis.title.x = element_text(size=14),
        axis.line.y.left = element_line(color = "black", size = .6),
        axis.text.y.left = element_text(color = "black", size = 12),
        axis.title.y = element_text(size=14))


##_1.1.3 - Variables (Plot)----
library(ggrepel)
#Ant richness and diversity of items in the litter, with interaction:
lit_m <- 
ggplot(tab_m, aes(x = Litt.Div., y = Sobs, fill = Plot)) +
  stat_smooth(aes(color=Plot), formula = y ~ x, method = "lm", fullrange=T, alpha = 0.2, size = 1.5)+
  geom_point(aes(color=Plot, shape=Plot), size=3)+
  geom_text_repel(aes(label = Months), direction = "y", hjust = 1, segment.size = 0.2, segment.color = 1, colour=1, size = 4) +
  scale_x_continuous(breaks=c(1.5,2,2.5,3,3.5), limits = c(1.5, 3.5))+
  scale_fill_manual(values= cor_ESECAFLOR)+
  scale_color_manual(values= cor_ESECAFLOR)+
  labs(x = "", y = "")+
  theme_classic()+
  theme(legend.position = "none",
        axis.line.x.bottom = element_line(color = "black", size = .6),
        axis.text.x.bottom = element_text(color = "black", size = 12),
        axis.title.x = element_text(size=14),
        axis.line.y.left = element_line(color = "black", size = .6),
        axis.text.y.left = element_text(color = "black", size = 12),
        axis.title.y = element_text(size=14))


##_1.2 - Composition----
names(tab_m)
ssp <- tab_m[,7:ncol(tab_m)]
var <- tab_m[,c(5,6,4)]

##_1.2.2 - PERMANOVA ----
library(vegan)
result <- adonis(ssp ~ (tab_m$Moisture +
                          tab_m$Litt.Div. +
                          tab_m$Plot +
                          tab_m$Biomass) ,
                 permutations=9999, method = "bray")
result


##_1.2.3 - Redundancy Analysis (RDA)----
ssp2 <- decostand(ssp, method="hellinger")
var2 <- decostand(var, method = "standardize")

RDA1 <- rda(ssp2, var2)

variaveis <- t(intersetcor(RDA1)[,1:3])
compo1 <- as.data.frame(summary(RDA1)$cont)
eigenv1 <- t(as.vector(compo1[1,1:3]))
cont1 <- t(as.vector(compo1[2,1:3]*100))
acum1 <- t(as.vector(compo1[3,1:3]*100))
teste <- as.data.frame(anova(RDA1))
p_valor <- as.data.frame(c(teste[1,4],NA,NA))

(tabela <- t(cbind(variaveis, eigenv1, cont1, acum1, p_valor)))
#write.table(tabela,"RDA_formigas_tabela.txt")

plot(RDA1, choices = c(1, 2), display = c("sp", "lc", "cn"))

##Figure ggplot2
value1 <- summary(RDA1)

#Units:
rda_units <- data.frame(value1$sites)[,1:2]
rda_units$Months <-tab_m$Months #Insert months
rda_units$Plot <-tab_m$Plot #Insert months

#Species:
rda_ssp <- data.frame(value1$species)[,1:2]
rda_ssp$species_names <-gsub("_","",rownames(rda_ssp))
#Adjust for scale species
rda_ssp$RDA1_nor <- scales::rescale(rda_ssp$RDA1, to=c(-1,1))
rda_ssp$RDA2_nor <- scales::rescale(rda_ssp$RDA2, to=c(-1,1))
1/mean(abs(range(rda_ssp$RDA1))) #Adjust x2 in RDA1 species = 6.28648
1/mean(abs(range(rda_ssp$RDA2))) #Adjust y2 in RDA2 species = 7.843531

#Arrows:
rda_arrows <- data.frame(value1$biplot)[,1:2]
rda_arrows$vari <- rownames(rda_arrows)
#Adjust for arrows
mult <-attributes(plot(RDA1)$biplot)$arrow.mul

#Colors:
cor_ssp <- "gray50" #Color species
cor_name_seta <- "red" #Color names variables


library(ggrepel)
rda_m <-
ggplot(rda_units, aes(x = RDA1, y = RDA2))+
  geom_hline(yintercept = 0, colour="gray", linetype="dashed")+ geom_vline(xintercept = 0, colour="gray", linetype="dashed")+
  stat_ellipse(aes(x=RDA1, y=RDA2, fill=Plot), geom="polygon", level=0.95, alpha=0.23)+
  geom_point(aes(colour=Plot, shape=Plot), size = 3)+
  #geom_text(data = rda_ssp, aes(x = RDA1_nor, y = RDA2_nor, label = species_names), check_overlap = F, colour = cor_ssp, size = 2.6, fontface="italic", alpha=.6)+ #Species name
  geom_point(data = rda_ssp, aes(x = RDA1_nor, y = RDA2_nor-.05), size = 1, shape=3, colour = cor_ssp, alpha=.6)+ #Species point
  geom_text_repel(aes(label = Months), direction = "y", hjust = 1, segment.size = 0.2, segment.color = 1, colour=1, size = 3.8) +
  geom_segment(data = rda_arrows, aes(x = 0, xend = mult * RDA1, y = 0, yend = mult * RDA2), arrow = arrow(length = unit(0.25, "cm")), colour = 1, size = 1)+
  geom_text(data = rda_arrows, aes(x = (mult + mult/10)*RDA1, y = (mult + mult/10)*RDA2, label = vari), size = 4.5, hjust = 0.7 , colour = cor_name_seta)+
  coord_cartesian(x = c(-1.1, 1), y = c(-.9, 1.1))+
  labs(x = "RDA 1 (20.02%)", y = "RDA 2 (9.30%)")+
  scale_x_continuous(sec.axis = sec_axis(~./6.28648))+
  scale_y_continuous(sec.axis = sec_axis(~./7.843531))+
  scale_fill_manual(values = cor_ESECAFLOR)+
  scale_color_manual(values = cor_ESECAFLOR)+
  scale_shape_manual(values = c(16, 17))+
  theme_light()+
  theme(legend.position = c(0.12, 0.92),
        legend.background = element_rect(color = "black", size = .5, linetype = "solid"),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        axis.title = element_text(size = 12),
        axis.line.x.bottom = element_line(color = "black", size = .8),
        axis.text.x.bottom = element_text(color = "black"),
        axis.line.y.left = element_line(color = "black", size = .8),
        axis.text.y.left = element_text(color = "black"),
        axis.line.x.top = element_line(color = cor_ssp, size = .8),
        axis.text.x.top = element_text(color = cor_ssp),
        axis.line.y.right = element_line(color = cor_ssp, size = .8),
        axis.text.y.right = element_text(color = cor_ssp))


##_1.2.4 - Multinomial Species Classification Method (CLAM) ----
library(vegan)

result_CLAM <- clamtest(ssp, tab_m$Plot, alpha = 0.05)
#write.table(result_CLAM,"CLAM_formigas_tabela_06-11-2019.txt")

summary(result_CLAM)
plot(result_CLAM)

# Abundance + 1
dados <- transform(result_CLAM,
                 Total_Con=1+(get("Total_Control")),
                 Total_Exp=1+(get("Total_Experiment")))

#Lines in figure
minval <- summary(result_CLAM)$minv
lineYvsgen <- minval[[1]]+1 #Line delimit specialist species Y
lineXvsgen <- minval[[2]]+1 #Line delimit specialist species X
Ymin <- minval[[1]][1,2]+1 #Line delimit rare species in Y
Xmin <- minval[[2]][1,1]+1 #Line delimit rare species in X

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

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

graphClam <-
ggplot(dados, aes(x=Total_Exp, y=Total_Con, colour=Classes))+
  geom_curve(alpha=0.2, aes(x=Xmin, y=1, xend=1, yend=Ymin), size=1.3, colour=cor_CLAM[4])+
  geom_line(data=lineYvsgen, aes(x=x, y=y), linetype=2, size=1.3, colour=cor_CLAM[2], alpha=0.5)+
  geom_line(data=lineXvsgen, aes(x=x, y=y), linetype=2, size=1.3, colour=cor_CLAM[3], alpha=0.5)+
  geom_point(aes(colour=Classes, shape=Classes), size=3.2, alpha=.8)+
  geom_text_repel(data=esp_s, aes(label=gsub("_","",Species)), show.legend=FALSE, fontface="italic", colour="gray30", box.padding = unit(0.2, 'lines'), segment.color = '#cccccc', size=3)+
  scale_shape_manual(values=c(18,19,17,15), 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,10,40,100)) +
  scale_y_continuous(trans = 'log10', breaks = c(1,4,10,40,100))+
  labs(x="Experiment (abundance + 1)", y="Control (abundance + 1)")+
  theme_classic()+
  theme(legend.position=c(0.1, 0.9),
        legend.title=element_blank(),
        legend.background = element_rect(color = "black", size = .5, linetype = "solid"),
        legend.key.size = unit(0, 'lines'),
        axis.title = element_text(size = 12),
        axis.line = element_line(color = "black", size=.8),
        axis.text.x.bottom = element_text(color = "black", size = 10),
        axis.text.y.left = element_text(color = "black", size = 10))


#####################
##2 - Analyses subplots----
tab_s <- read.table("ants_subplots_ESECAFLOR_2011-2012.txt",h=T,sep="\t")

##_Evaluating spatial autocorrelation----
geo <- cbind(tab_s$x_cor[1:25], tab_s$y_cor[1:25])
samples.dist <- as.matrix(dist(geo))

#Reverse
samples.dist.inv <- 1/samples.dist
diag(samples.dist.inv) <- 0

library(ape)
Moran.I(tab_s$Sobs[1:25], samples.dist.inv ,alternative="greater") #Site Control
Moran.I(tab_s$Sobs[26:50], samples.dist.inv ,alternative="greater") #Site Experiment
#None of the sites has autocorrelation
#View plot:

library(ggplot2)
library(scales)

ggplot(tab_s, aes(x = x_cor, y = y_cor, fill = Sobs)) +
  geom_tile() +
  geom_text(aes(label = Sobs), size = 7, colour = "black") +
  facet_wrap(~ Plot)+
  labs(x="Meters", y="Meters", fill="Richness")+
  scale_fill_gradient2(low = "blue", high = "red", limits = c(min(tab_s$Sobs), max(tab_s$Sobs)), midpoint = median(tab_s$Sobs)) +
  theme(panel.background = element_rect(fill="transparent"),
        axis.ticks = element_blank(),
        panel.border = element_blank(),
        strip.text.x = element_text(size=16),
        strip.background = element_rect(colour="black", fill="white"))+
  coord_equal()


##_Correlation of variables ----
var_pairs <- tab_s[,c(6:11, 1)]

library(GGally)
p1 <- ggpairs(var_pairs, columns = 1:6,
              lower = list(continuous = my_fn),
              diag = list(continuous = my_fn2), 
              upper = list(continuous=cor_fun))

cor_s <- p1+theme_minimal()+
  theme(legend.position = "none", 
        panel.grid.minor = element_blank(),
        axis.line.x.bottom = element_line(color = "black", size = .8),
        axis.text.x.bottom = element_text(color = "black"),
        axis.line.y.left = element_line(color = "black", size = .8),
        axis.text.y.left = element_text(color = "black"),
        panel.border = element_rect(linetype = "solid", colour = "gray50", fill = NA),
        strip.text.x = element_text(size=12),
        strip.text.y = element_text(size=12),
        strip.background = element_rect(colour="black", fill="white"))


##_2.1 - Richness ----
##_2.1.1 - Best model and selection AIC ----
library(nlme)
mod2 <- nlme::lme(Sobs ~ (Litt.Div. + Biomass + Moisture + Trees + DBH) * Plot, random=~1|Plot, method = "REML", data=tab_s)

options(na.action = "na.fail")
library(MuMIn)
selec <- dredge(mod2, rank="AICc", extra = "adjR^2")
(mod_acaic <- summary(model.avg(selec))) #Resultado do Akaik
importance(selec) #Substituto para a importância das variáveis no modelo quando há interação. Refazer usando isso em detrimento do hier.part


##_2.1.2 - Importance of variables (Plot)----
plot(importance(selec)); abline(h=0.7, col=2)

import_var_a <- data.frame(importance(selec))
import_var_a$nomes_var_a <- row.names(import_var_a)
import_var_a$import_a <- round(import_var_a$importance.selec., 3)
import_var_a2 <- import_var_a[order(import_var_a$import_a, decreasing = T),]

imp_s <-
  ggplot(data=import_var_a2, aes(x=reorder(nomes_var_a, +import_a), y=import_a)) +
  geom_bar(position="dodge",stat="identity", colour="black") + 
  geom_hline(aes(yintercept=0.7), color="red", linetype="dashed", size=1)+
  scale_y_continuous(breaks = c(0.1,0.4,0.7,1))+
  labs(x = "Temporal variables", y = "Importance", size = 20)+
  theme_classic()+
  coord_flip()+
  theme(axis.line = element_line (size=.4),
        axis.line.x.bottom = element_line(color = "black", size = .6),
        axis.text.x.bottom = element_text(color = "black", size = 12),
        axis.title.x = element_text(size=14),
        axis.line.y.left = element_line(color = "black", size = .6),
        axis.text.y.left = element_text(color = "black", size = 12),
        axis.title.y = element_text(size=14))


##_2.1.3 - Variables (Plot)----
library(ggrepel)
#Ant richness and diversity of items in the litter, with interaction:
lit_s <- 
  ggplot(tab_s, aes(x = Litt.Div., y = Sobs, fill = Plot)) +
  stat_smooth(aes(color=Plot), formula = y ~ x, method = "lm", fullrange=T, alpha = 0.2, size = 1.5)+
  geom_point(aes(color=Plot, shape=Plot), size=3)+
  scale_x_continuous(breaks=c(1.5,2,2.5,3,3.5), limits = c(1.5, 3.5))+
  scale_fill_manual(values= cor_ESECAFLOR)+
  scale_color_manual(values= cor_ESECAFLOR)+
  labs(x = "Litter diversity", y = "")+
  theme_classic()+
  theme(legend.position = "none",
        axis.line.x.bottom = element_line(color = "black", size = .6),
        axis.text.x.bottom = element_text(color = "black", size = 12),
        axis.title.x = element_text(size=14),
        axis.line.y.left = element_line(color = "black", size = .6),
        axis.text.y.left = element_text(color = "black", size = 12),
        axis.title.y = element_text(size=14))


#Ant richness and tree richness
tre_s <-
  ggplot(tab_s, aes(x = Trees, y = Sobs)) +
  stat_smooth(method = "lm", formula = y ~ x, fullrange=T, alpha = 0.3, size = 1.5, color=1)+
  geom_point(aes(color=Plot, shape=Plot), size=3)+
  labs(x = "Tree richness", y = "Ants richness", size = 20)+
  scale_x_continuous(breaks=c(13,15,17,19,21), limits = c(12, 23))+
  scale_fill_manual(values= cor_ESECAFLOR)+
  scale_color_manual(values= cor_ESECAFLOR)+
  theme_classic()+
  theme(legend.position = "none",
        axis.line.x.bottom = element_line(color = "black", size = .6),
        axis.text.x.bottom = element_text(color = "black", size = 12),
        axis.title.x = element_text(size=14),
        axis.line.y.left = element_line(color = "black", size = .6),
        axis.text.y.left = element_text(color = "black", size = 12),
        axis.title.y = element_text(size=14))


##_2.2 - Composition----
names(tab_s)
ssp <- tab_s[,12:ncol(tab_s)]
var <- tab_s[,c(6:11)]


##_2.2.2 - PERMANOVA ----
library(vegan)
result <- adonis(ssp ~ (var$Moisture +
                          var$Litt.Div. +
                          tab_s$Plot +
                          var$Biomass +
                          var$Resp. +
                          var$Trees +
                          var$DBH),
                 permutations=9999, method = "bray") #Importância biológica
result


##_2.2.3 - Redundancy Analysis (RDA)----
ssp2 <- decostand(ssp, method="hellinger")
var2 <- decostand(var, method = "standardize")

RDA1 <- rda(ssp2, var2)

variaveis <- t(intersetcor(RDA1)[,1:3])
compo1 <- as.data.frame(summary(RDA1)$cont)
eigenv1 <- t(as.vector(compo1[1,1:3]))
cont1 <- t(as.vector(compo1[2,1:3]*100))
acum1 <- t(as.vector(compo1[3,1:3]*100))
teste <- as.data.frame(anova(RDA1))
p_valor <- as.data.frame(c(teste[1,4],NA,NA))

(tabela <- t(cbind(variaveis, eigenv1, cont1, acum1, p_valor)))
#write.table(tabela,"RDA_formigas_tabela.txt")

plot(RDA1, choices = c(1, 2), display = c("sp", "lc", "cn"))

##Figure ggplot2
value1 <- summary(RDA1)

#Units:
rda_units <- data.frame(value1$sites)[,1:2]
rda_units$Subplots <-tab_s$Subplots #Insert subplots
rda_units$Plot <-tab_s$Plot #Insert Sites

#Species:
rda_ssp <- data.frame(value1$species)[,1:2]
rda_ssp$species_names <- gsub("_","",rownames(rda_ssp))
#Adjust for scale species
rda_ssp$RDA1_nor <- scales::rescale(rda_ssp$RDA1, to=c(-1,1))
rda_ssp$RDA2_nor <- scales::rescale(rda_ssp$RDA2, to=c(-1,1))
1/mean(abs(range(rda_ssp$RDA1))) #Adjust x2 in RDA1 species = 5.867226
1/mean(abs(range(rda_ssp$RDA2))) #Adjust y2 in RDA2 species = 5.759028

#Arrows:
rda_arrows <- data.frame(value1$biplot)[,1:2]
rda_arrows$vari <- rownames(rda_arrows)
#Adjust for arrows
mult <- attributes(plot(RDA1)$biplot)$arrow.mul

#Colors:
cor_ssp <- "gray50" #Color species
cor_name_seta <- "red" #Color names variables


library(ggrepel)
rda_s <-
ggplot(rda_units, aes(x = RDA1, y = RDA2))+
  geom_hline(yintercept = 0, colour="gray", linetype="dashed")+ geom_vline(xintercept = 0, colour="gray", linetype="dashed")+
  stat_ellipse(aes(x=RDA1, y=RDA2, fill=Plot), geom="polygon", level=0.95, alpha=0.23)+
  geom_point(aes(colour=Plot, shape=Plot), size = 3)+
  #geom_text(data = rda_ssp, aes(x = RDA1_nor, y = RDA2_nor, label = species_names), check_overlap = F, colour = cor_ssp, size = 2.6, fontface="italic", alpha=.6)+ #Species name
  geom_point(data = rda_ssp, aes(x = RDA1_nor+0.09, y = RDA2_nor-0.38), size = 1, shape=3, colour = cor_ssp, alpha=.6)+ #Species point
  geom_segment(data = rda_arrows, aes(x = 0, xend = mult*RDA1, y = 0, yend = mult*RDA2), arrow = arrow(length = unit(0.25, "cm")), colour = 1, size = 1)+
  geom_text(data = rda_arrows, aes(x = (mult + mult/10)*RDA1, y = (mult + mult/10)*RDA2, label = vari), size = 4.5, hjust = 0.7, colour = cor_name_seta)+
  coord_cartesian(x = c(-1.1, 1.1), y = c(-1.3, 0.9))+
  labs(x = "RDA 1 (6.02%)", y = "RDA 2 (3.70%)")+
  scale_x_continuous(sec.axis = sec_axis(~./5.867226))+
  scale_y_continuous(sec.axis = sec_axis(~./5.759028))+
  scale_fill_manual(values = cor_ESECAFLOR)+
  scale_color_manual(values = cor_ESECAFLOR)+
  scale_shape_manual(values = c(16, 17))+
  theme_light()+
  theme(legend.position = c(0.1, 0.9),
        legend.background = element_rect(color = "black", size = .5, linetype = "solid"),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        axis.title = element_text(size = 12),
        axis.line.x.bottom = element_line(color = "black", size = .8),
        axis.text.x.bottom = element_text(color = "black"),
        axis.line.y.left = element_line(color = "black", size = .8),
        axis.text.y.left = element_text(color = "black"),
        axis.line.x.top = element_line(color = cor_ssp, size = .8),
        axis.text.x.top = element_text(color = cor_ssp),
        axis.line.y.right = element_line(color = cor_ssp, size = .8),
        axis.text.y.right = element_text(color = cor_ssp))


#####################
##3 - Save Figures----
library(gridExtra)

##Figure 1
imp_m.1 <- imp_m + ggtitle("(A)") + theme(plot.title=element_text(hjust=0))
imp_s.1 <- imp_s + ggtitle("(C)") + theme(plot.title=element_text(hjust=0))
lit_m.1 <- lit_m + ggtitle("(B)") + theme(plot.title=element_text(hjust=0))
lit_s.1 <- lit_s + ggtitle("(D)") + theme(plot.title=element_text(hjust=0))
library(grid)
matriz <- rbind(c(rep(1,4),NA, rep(2,10)),
                c(rep(3,4),NA, rep(4,10)))
png(file = "Fig1_22-07-2020.png", width = 3500, height = 2500, units = "px", res = 300)
grid.arrange(imp_m.1, lit_m.1, imp_s.1, lit_s.1, layout_matrix = matriz)
grid.text("Ants richnes", x = unit(0.33, "npc"),y = unit(.5, "npc"),gp = gpar(fontsize=20, fontfamily="sans"), rot=90)
dev.off()


#Figure 2
png(file = "Fig2_22-07-2020.png", width = 3900, height = 2200, units = "px", res = 300)
grid.arrange(rda_m +coord_cartesian(x = c(-1.1, .8), y = c(-.95, .8)) +ggtitle("(A)"),
             rda_s +coord_cartesian(x = c(-1.1, .7), y = c(-1, .7)) +ggtitle("(B)") +theme(legend.position="none"),
             ncol=2)
dev.off()


#Figure 3
png(file = "Fig3_22-07-2020.png", width = 2400, height = 1700, units = "px", res = 300)
graphClam
dev.off()


#FigureS2
png(file = "FigS2_22-07-2020.png", width = 2000, height = 1500, units = "px", res = 300)
cor_m
dev.off()


#FigureS3
png(file = "FigS3_22-07-2020.png", width = 2500, height = 1800, units = "px", res = 300)
cor_s
dev.off()


#Figure S4
png(file = "FigS4_22-07-2020.png", width = 2200, height = 1300, units = "px", res = 300)
tre_s
dev.off()


##4 - CITATION----
#Statistical analysis:
citation()
citation("nlme")
citation("MuMIn")
citation("vegan")
#Figures:
citation("ggplot2")
citation("ggrepel")
citation("GGally")
citation("gridExtra")

dev.off()