# R script for comparing species diversity between periods
# Use the raw data to make various data sets of species abundances

# Load libraries
library(iNEXT)
library(ggplot2)
library(codyn)
library(permute)
library(lattice)
library(permute)
library(readxl)

# Convert data set 1 to data frame 1 for rarefaction analyses
df1 <- data.frame(Abundance_BioDivCL_1_)

# Convert data set 2 to data frame 2 for diversity analyses
df2 <- data.frame(Diversity_BioDivCL_1_)

# Convert data set 3 to data frame 3 for rank abundance analyses
df3 <- data.frame(RankAbunCurve_BioDivCL_1_)

# Calculate and compare the Shannon Diversity Index by site
community_diversity(df2, abundance.var = "abundance", replicate.var = "site", metric = c("Shannon"))

# Calculate evenness by site
community_structure(df2, time.var = NULL, abundance.var = "abundance", replicate.var = "site", metric = "SimpsonEvenness")

# Create rank abundance curves by site
rad <- radfit(df3)
plot(rad)

# Calculate the area difference in rank abundance curves by site
curve_difference(df2, species.var = "species", abundance.var = "abundance", replicate.var = "site")

# Calculate differences between rank abundance curves by site for four measures: richness, evenness, rank, and species composition
RAC_difference(df2, species.var = "species", abundance.var = "abundance", replicate.var = "site")

# Simulate species accumulation by number of individuals using bootstrap methods
out <- iNEXT(df1, datatype="abundance", endpoint = 2000, conf = 0.95, nboot = 1000, q=0)

# Generate rarefaction curves with extrapolation for species richness with 95% CI interval
ggiNEXT(out) + theme_classic() + theme(legend.position = "top") + scale_y_continuous(name="Species diversity")

# Estimate species richness via the Chao method
ChaoRichness(df1, datatype="abundance", conf = 0.95)

# Make subsets of data for arthropods, birds, herpetofuana, and plants for both time periods
# Import new data sets and convert them to data frames 4, 5, 6, & 7
# Convert data sets to data frames for rarefaction analyses
df4 <- data.frame(Abundance_BioDivCL_arthro)
df5 <- data.frame(Abundance_BioDivCL_birds)
df6 <- data.frame(Abundance_BioDivCL_herps)
df7 <- data.frame(Abundance_BioDivCL_plants)

# Simulate species accumulation by number of individuals by taxonomic group using bootstrap methods
outarthro <- iNEXT(df4, datatype="abundance", endpoint = 2000, conf = 0.95, nboot = 1000, q=0)
outbirds <- iNEXT(df5, datatype="abundance", endpoint = 2000, conf = 0.95, nboot = 1000, q=0)
outherps <- iNEXT(df6, datatype="abundance", endpoint = 2000, conf = 0.95, nboot = 1000, q=0)
outplants <- iNEXT(df7, datatype="abundance", endpoint = 2000, conf = 0.95, nboot = 1000, q=0)

# Generate rarefaction curves by taxonomic group with extrapolation for species richness with 95% CI interval
ggiNEXT(outarthro) + theme_classic() + theme(legend.position = "top") + scale_y_continuous(name="Species diversity")
ggiNEXT(outbirds) + theme_classic() + theme(legend.position = "top") + scale_y_continuous(name="Species diversity")
ggiNEXT(outherps) + theme_classic() + theme(legend.position = "top") + scale_y_continuous(name="Species diversity")
ggiNEXT(outplants) + theme_classic() + theme(legend.position = "top") + scale_y_continuous(name="Species diversity")

# Estimate species richness via the Chao method by taxonomic group
ChaoRichness(df4, datatype="abundance", conf = 0.95)
ChaoRichness(df5, datatype="abundance", conf = 0.95)
ChaoRichness(df6, datatype="abundance", conf = 0.95)
ChaoRichness(df7, datatype="abundance", conf = 0.95)