########################################
### Chapter 16 - Community Structure ###
### Species diversity and food webs ###
########################################

# Load libraries
library("tidyverse")
library("vegan")

# Load data set
data(BCI)

# prepare abundnace data set
Abundance <- BCI %>% 
  # add site names
  mutate(site = 1:50) %>% 
  gather(key = species, value = abundance, -site) %>% 
  filter(!abundance == 0) %>% 
  # reduce the data set to 2 sites (18 and 35)
  filter(site %in% c(35, 18)) %>% 
  # arrange the species according to their abundance, within each site
  arrange(site, -abundance) %>% 
  # Tweak the data a bit to make the them more uneven
  mutate(abundance = ifelse(site == 18 & species == "Faramea.occidentalis", 18, abundance),
         abundance = ifelse(site == 18 & species == "Oenocarpus.mapora", 18, abundance),
         abundance = ifelse(site == 18 & species == "Trichilia.tuberculata", 17, abundance))

# 1) Rank-abundance plot
RelAbundance <- Abundance %>% 
  # calculate the total nr of species, relative abundnace and rank the species according to their abundance.
  group_by(site) %>% 
  mutate(totalN = sum(abundance),
         relAbunance = abundance / totalN,
         rankAbundance = 1:n())

# Make a rank-abundance plot
ggplot(RelAbundance, aes(x = rankAbundance, y = log10(relAbunance), colour = factor(site)))  +
  geom_point() +
  geom_line()



# 2) Species richness and evenness
RelAbundance %>% 
  group_by(site) %>% 
  summarise(# species richness
            richnessS = n(),
            
            # Simpson indexes
            simpsonD = sum(relAbunance^2),
            simpsonIndexofD = 1 - simpsonD,
            simpsonRecipIndex = 1 / simpsonD,
            EvennessD = simpsonRecipIndex / richnessS,
            
            # Shannon index
            H = -sum(relAbunance * log(relAbunance)),
            Hmax = log(richnessS),
            EvennessH = H / Hmax,
            
            # Check values with vegan package
            diversityD = diversity(abundance, index = "simpson"),
            diversityH = diversity(abundance, index = "shannon"))



# 4) Foodweb, connectance and linkage density
foodweb <- tibble(FoodWeb = c("A", "B"),
       SpeciesRichness = c(11, 11),
       NrLinkage = c(14, 18))

# Calculate connectance and linkage density
foodweb %>% 
  mutate(MaxLinkage = SpeciesRichness^2,
         Connectance = NrLinkage / SpeciesRichness^2,
         LinkageDensity = NrLinkage / SpeciesRichness)