### Source: SIA-BRA
### Authors: Thaís Diniz-Reis and Luís A. Martinelli
### Journal: Global Ecology and Biogeography
### Year: 2022
### DOI: https://doi.org/10.1111/geb.13449

rm(list=ls())
graphics.off()
set.seed(1)

install.packages(c('readxl', 'dplyr', 'modeest','ggplot2', 'tidyverse', 'tidytext', 'ggtext', 'ggforce', 'ggdist', 'patchwork', 'ggpubr', 'lme4', 'rgdal','ggridges', 'huxtable', 'flextable', 'officer', 'multcomp', 'MuMIn', 'lmerTest', 'evaluate', 'effectsize', 'broom.mixed', 'viridis', 'hrbrthemes'))

library(readxl)     
library(dplyr)       
library(modeest)     
library(ggplot2)     
library(tidyverse)   
library(tidytext)    
library(ggtext)     
library(ggforce)     
library(ggdist)      
library(patchwork)   
library(ggpubr)      
library(lme4)        
library(rgdal)       
library(ggridges)    
library(huxtable)    
library(flextable)   
library(officer)    
library(multcomp)    
library(MuMIn)       
library(lmerTest)    
library(evaluate)   
library(effectsize)  
library(broom.mixed) 
library(viridis)   
library(hrbrthemes)  

### Entering data ###

# Data directory
choose.dir()
getwd()

# Data upload
data<-read.csv("SIA_BRA.csv", header=T, stringsAsFactors = T)

options(digits=3)

# Reclassification of original variables into a new variables
data$area<- ifelse(data$habitat == 'Terrestrial', 'Terrestrial',
                   ifelse(data$habitat == 'Freshwater', 'Freshwater', 'Marine'))
data$area <- as.factor (data$area)

data$paper <- paste(data$author,data$year)
data$paper <- as.factor(data$paper)

summary(data)

# Filter data into groups
ver <- data %>% filter (group == 'Vertebrates')
inv <- data %>% filter (group == 'Invertebrates')


#### Basic stats ####

id<-data %>% 
  group_by(group,area)

stat.c<-summarise(id,
                  avg=mean(d13C, na.rm=TRUE),
                  sd=sd(d13C,na.rm=TRUE),
                  median=median(d13C, na.rm=TRUE),
                  mode=mlv(d13C, na.rm=TRUE),
                  IQR=IQR(d13C, na.rm=TRUE),
                  min=min(d13C, na.rm=TRUE),
                  max=max(d13C, na.rm=TRUE),
                  count=n())
stat.c

obs13<-tibble(stat.c)
obs13
obs13<-mutate(obs13, conf.avg=(sd/sqrt(count))*1.96)
obs13

obs13<-mutate(obs13, "ci+95%" = avg+conf.avg)
obs13<-mutate(obs13, "ci-95%" = avg-conf.avg)
obs13_ft <- flextable(
  data = head(obs13), 
  col_keys = c('group',"area", "avg", "ci-95%", "ci+95%", "count"))
obs13_ft
save_as_docx("Table.carbon" = obs13_ft, path="C:/Users/thais/OneDrive/Área de Trabalho/trabalho/Table.carbon.docx")
print(obs13_ft, preview = "docx")

stat.n<-summarise(id,
                  avg=mean(d15N, na.rm=TRUE),
                  sd=sd(d15N,na.rm=TRUE),
                  median=median(d15N, na.rm=TRUE),
                  mode=mlv(d15N, na.rm=TRUE),
                  IQR=IQR(d15N, na.rm=TRUE),
                  min=min(d15N, na.rm=TRUE),
                  max=max(d15N, na.rm=TRUE),
                  count=n())
stat.n

obs15<-tibble(stat.n)
obs15
obs15<-mutate(obs15, conf.avg=(sd/sqrt(count))*1.96)
obs15

obs15<-mutate(obs15, "ci+95%" = avg+conf.avg)
obs15<-mutate(obs15, "ci-95%" = avg-conf.avg)
obs15_ft <- flextable(
  data = head(obs15), 
  col_keys = c('group',"area", "avg", "ci-95%", "ci+95%", "count"))
obs15_ft
save_as_docx("Table.nirogen" = obs15_ft, path="C:/Users/thais/OneDrive/Área de Trabalho/trabalho/Table.nitrogen.docx")
print(obs15_ft, preview = "docx")


#### Figures ####

## Customized Plot

theme_set(theme_minimal(base_size = 12, base_family = "Open Sans"))
minhas_cores <-c('skyblue', 'darkblue','#75533d')

# Map: sites
# Figure 1

install.packages(c('plyr','sf','geobr','ggspatial'))

library(sf) # vetores
library(geobr) # Limites IBGE
library(ggplot2) # Graficos e mapas
library(ggspatial) # Elementos espaciais no ggplot2

# IBGE shape files upload
polygons_biomes <- sf::st_read("lm_bioma_250.shp")
polygons_brazil <- sf::st_read("gadm36_BRA_0.shp")
polygons_coast <- sf::st_read("Sistema_Costeiro_Marinho.shp")
polygons_southamerica <- sf::st_read("AMERICA_SUL.shp")
polygons_hidro <- sf::st_read("hid_trecho_massa_dagua_a.shp")
polygons_biomes$Bioma[polygons_biomes$Bioma=="Mata Atlântica"] <- "Atlantic"
polygons_biomes$Bioma[polygons_biomes$Bioma=="Amazônia"] <- "Amazon"
biomes_legend <- factor(c('Amazon','Atlantic','Caatinga', 'Cerrado', 'Pampa', 'Pantanal','Coastal-marine'))

# Map of the sampling sites
mapa <- ggplot() +
  geom_sf(data = polygons_southamerica, color = "darkgray", lwd = 0.35, fill = "gray88", show.legend = FALSE) +
  geom_sf(data = polygons_biomes, color = "transparent", aes(fill= Bioma), lwd = 0.2, alpha=0.7, show.legend = TRUE) +
  geom_sf(data = polygons_coast, color = "darkgray", lwd = 0.1, fill = "#16168C", alpha=0.5,show.legend = FALSE) +
  geom_sf(data = polygons_hidro, color = "transparent", lwd = 0.2, fill = 'skyblue', alpha=0.5,show.legend = FALSE) +
  geom_sf(data = polygons_brazil, color = "transparent", lwd = 0.1, fill = "transparent", show.legend = FALSE) +
  scale_fill_manual(values = c('#322318', "#50382a", "#75533d", "#a97a5b", "#c5a591", "#e1d1c6"))+
  geom_point(data = data, aes(x = long, y = lat), col= 'black', fill='white', size = 1, shape=21, stroke = 0.5, alpha = 1, show.legend = FALSE) +
  ylab("Latitude") +
  xlab("Longitude") +
  labs(fill='Biome')+
  coord_sf(xlim = c(-76, -28), ylim = c(-36, 8), expand = FALSE) +
  theme(axis.text.y = element_text(size = 8, color = "black"),
        axis.text.x = element_text(size = 8, color = "black"),
        axis.title.y = element_text(size = 10, color = "black"),
        axis.title.x = element_text(size = 10, color = "black"),
        legend.title=element_text(size=8, color = "black"),
        legend.text =element_text(size=8, color = "black"),
        legend.position = 'bottom',
        panel.background = element_rect(fill = "white"),
        panel.grid.major = element_line(color = "gray96", size = 0.50),
        panel.grid.minor = element_line(color = "gray96", size = 0.30),
        axis.line = element_line(color = "black", size = 0.5),
        panel.border = element_rect(color = "black", fill = NA, size = 0.5))+ 
  annotation_scale(location = "br") +
  annotation_north_arrow(location = "br", which_north = "true", 
                         pad_x = unit(0, "cm"), pad_y = unit(.8, "cm"),
                         style = north_arrow_fancy_orienteering)
plot(mapa)
ggsave("Figure 1.tiff", plot = mapa, units="px", width=1500, height=1500, dpi=300)


# Histograms
# Figure 2

#vertebrates
ver.c<- ggplot(ver,  aes(x=d13C, colour=area, fill=area, scales='free_y')) +
  geom_histogram(alpha=0.5, binwidth = 0.8, position="identity") + 
  scale_color_manual(values = minhas_cores) + 
  scale_fill_manual(values = minhas_cores)+
  scale_x_continuous(name=expression(paste(delta^{13}, "C (\u2030)")), 
                     breaks = seq(-45,0, 5),
                     limits=c(-45, 0)) + 
  scale_y_continuous(name= "Frequency of values",
                     breaks = seq(0,1000,250),
                     limits=c(0, 1000)) +
  annotate("text", x=-40, y=1000, label= "Vertebrates", size = 2.3)+
  theme(axis.line.x = element_line(size = 0.3, colour = "black"),
        axis.line.y = element_line(size = 0.3, colour = "black"),
        axis.line = element_line(size=1, colour = "black"),
        axis.ticks= element_line(colour="black"),
        axis.ticks.length = unit(.1, "cm"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(fill=NA, linetype=1),
        panel.background = element_blank(),
        axis.text.x=element_text(size=5.7, colour="black"),
        axis.text.y=element_text(size=5.7, colour="black"),
        axis.title.y=element_text(size=5.7, colour="black"),
        axis.title.x = element_text (size = 5.7, colour="black"), 
        legend.position = "bottom",
        legend.text = element_text(size=6, colour="black"),
        legend.title = element_blank())
ver.c

ver.n<- ggplot(ver,  aes(x=d15N, colour=area, fill=area, scales='free_y')) +
  geom_histogram(alpha=0.5, binwidth = 0.8, position="identity") + 
  scale_color_manual(values = minhas_cores) + 
  scale_fill_manual(values = minhas_cores)+
  scale_x_continuous(name=expression(paste(delta^{15}, "N (\u2030)")), 
                     breaks = seq(-10,30,5),
                     limits=c(-10, 30)) + 
  scale_y_continuous(name= "Frequency of values",
                     breaks = seq(0,600,100),
                     limits=c(0, 600)) +
  annotate("text", x=-6, y=600, label= "Vertebrates", size = 2.3)+
  theme(axis.line.x = element_line(size = 0.3, colour = "black"),
        axis.line.y = element_line(size = 0.3, colour = "black"),
        axis.line = element_line(size=1, colour = "black"),
        axis.ticks= element_line(colour="black"),
        axis.ticks.length = unit(.1, "cm"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(fill=NA, linetype=1),
        panel.background = element_blank(),
        axis.text.x=element_text(size=5.7, colour="black"),
        axis.text.y=element_text(size=5.7, colour="black"),
        axis.title.y=element_text(size=5.7, colour="black"),
        axis.title.x = element_text (size = 5.7, colour="black"), 
        legend.position = "bottom",
        legend.text = element_text(size=6, 
                                   
                                   colour="black"),
        legend.title = element_blank())
ver.n


#invertebrates
inv.c<- ggplot(inv,  aes(x=d13C, colour=area, fill=area, scales='free_y')) +
  geom_histogram(alpha=0.5, binwidth = 0.8, position="identity") + 
  scale_color_manual(values = minhas_cores) + 
  scale_fill_manual(values = minhas_cores)+
  scale_x_continuous(name=expression(paste(delta^{13}, "C (\u2030)")), 
                     breaks = seq(-45,0, 5),
                     limits=c(-45, 0)) + 
  scale_y_continuous(name= "Frequency of values",
                     breaks = seq(0,1000,250),
                     limits=c(0, 1000)) +
  annotate("text", x=-40, y=1000, label= "Invertebrates", size = 2.3)+
  theme(axis.line.x = element_line(size = 0.3, colour = "black"),
        axis.line.y = element_line(size = 0.3, colour = "black"),
        axis.line = element_line(size=1, colour = "black"),
        axis.ticks= element_line(colour="black"),
        axis.ticks.length = unit(.1, "cm"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(fill=NA, linetype=1),
        panel.background = element_blank(),
        axis.text.x=element_text(size=5.7, colour="black"),
        axis.text.y=element_text(size=5.7, colour="black"),
        axis.title.y=element_text(size=5.7, colour="black"),
        axis.title.x = element_text (size = 5.7, colour="black"), 
        legend.position = "bottom",
        legend.text = element_text(size=6, colour="black"),
        legend.title = element_blank())
inv.c


inv.n<- ggplot(inv, aes(x=d15N, colour=area, fill=area, scales='free_y')) +
  geom_histogram(alpha=0.5, binwidth = 0.8, position="identity") + 
  scale_color_manual(values = minhas_cores) + 
  scale_fill_manual(values = minhas_cores)+
  scale_x_continuous(name=expression(paste(delta^{15}, "N (\u2030)")), 
                     breaks = seq(-10,30,5),
                     limits=c(-10, 30)) + 
  scale_y_continuous(name= "Frequency of values",
                     breaks = seq(0,1300,250),
                     limits=c(0, 1300)) +
  annotate("text", x=-5, y=1300, label= "Invertebrates", size = 2.3)+
  theme(axis.line.x = element_line(size = 0.3, colour = "black"),
        axis.line.y = element_line(size = 0.3, colour = "black"),
        axis.line = element_line(size=1, colour = "black"),
        axis.ticks= element_line(colour="black"),
        axis.ticks.length = unit(.1, "cm"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(fill=NA, linetype=1),
        panel.background = element_blank(),
        axis.text.x=element_text(size=5.7, colour="black"),
        axis.text.y=element_text(size=5.7, colour="black"),
        axis.title.y=element_text(size=5.7, colour="black"),
        axis.title.x = element_text (size = 5.7, colour="black"), 
        legend.position = "bottom",
        legend.text = element_text(size=6, colour="black"),
        legend.title = element_blank())
inv.n

histogram<- ggarrange(ver.c,ver.n, inv.c, inv.n, ncol=2, nrow=2, common.legend = T, legend = 'bottom')
histogram
ggsave("Figure 2.tiff", plot = histogram, units="px", width=1500, height=1500, dpi=300)


## Modeling

data[, 'biome'] <- as.factor(data[, 'biome'])
str(data)

data[, 'habitat'] <- as.factor(data[, 'habitat'])
str(data)

#Vertebrates
data<- ver

# lmer model
md13C<-lmer(d13C~(1|biome/habitat/paper) + (1|tissue) + (1|diet), data=data)
md15N<-lmer(d15N~(1|biome/habitat/paper) + (1|tissue) + (1|diet), data=data) 

summary(md13C)
summary(md15N)

r.squaredGLMM(md13C)
r.squaredGLMM(md15N)

resmd13C<-residuals(md13C)
resmd15N<-residuals(md15N)

histmd13C<-hist(as.numeric(resmd13C),prob=TRUE, col="green", breaks = 50,  xlab="Residues", main="")
curvemd13C<-curve(dnorm(x, mean=mean(resmd13C), sd=sd(resmd13C)), add=TRUE, col="darkblue", lwd=2)

histmd15N<-hist(as.numeric(resmd15N),prob=TRUE, col="green", breaks = 50,  xlab="Residues", main="")
curvemd15N<-curve(dnorm(x, mean=mean(resmd15N), sd=sd(resmd15N)), add=TRUE, col="red", lwd=2)

#Invertebrates
data <- inv

# lmer model
md13C<-lmer(d13C~(1|biome/habitat/paper) + (1|tissue) + (1|diet), data=data)
md15N<-lmer(d15N~(1|biome/habitat/paper) + (1|tissue) + (1|diet), data=data) 

summary(md13C)
summary(md15N)

r.squaredGLMM(md13C)
r.squaredGLMM(md15N)

resmd13C<-residuals(md13C)
resmd15N<-residuals(md15N)

histmd13C<-hist(as.numeric(resmd13C),prob=TRUE, col="green", breaks = 50,  xlab="Residues", main="")
curvemd13C<-curve(dnorm(x, mean=mean(resmd13C), sd=sd(resmd13C)), add=TRUE, col="darkblue", lwd=2)

histmd15N<-hist(as.numeric(resmd15N),prob=TRUE, col="green", breaks = 50,  xlab="Residues", main="")
curvemd15N<-curve(dnorm(x, mean=mean(resmd15N), sd=sd(resmd15N)), add=TRUE, col="red", lwd=2)


data<-read.csv("SIA_BRA_model.csv", header=T, stringsAsFactors = T)
summary(data)

mod.plot<- data %>%  ggplot(aes(variable))+
  geom_bar(alpha=0.5, binwidth = 0.8)
ggplot(aes(x=variance_d13C, colour=variable, fill=variable, scales='free_y')) +
  
  mod.plot

## Baseline 
## Table 
data<-read.csv("SIA_BRA_baseline.csv", header=T)

id<-group_by(data,habitat) 
statd15N<-summarise(id,
                    avg=mean(d15N, na.rm=TRUE),
                    Sd=sd(d15N,na.rm=TRUE),
                    med=median(d15N, na.rm=TRUE),
                    IQR=IQR(d15N, na.rm=TRUE),
                    min=min(d15N, na.rm=TRUE), 
                    max=max(d15N, na.rm=TRUE),
                    count=n())
statd15N


modeld15N<-lm(d15N ~ habitat, data=data)
summary(modeld15N)
r.squaredGLMM(modeld15N)
plot(modeld15N)

histmd15N<-hist(as.numeric(resmd15N),prob=TRUE, col="green", breaks=30, xlab="Resídues - d15N", main="")
curvemd15N<-curve(dnorm(x, mean=mean(resmd15N), sd=sd(resmd15N)), add=TRUE, col="red", lwd=2)