# Rare Community Index Supplemental Code


### Simulated data for example

# Simulated data consists of 600 sampled sites x1-x600
# Sites contain between 1-20 species detected at each
# There are 60 possible species, labelled sp1-sp60, which are observable in a site

set.seed(123) # Set seed for reproducibility
num_sites <- 600 # Number of sites sampled
num_species <- 60 # Number of species observable
species_list <- paste("sp", 1:num_species, sep = "") # Generate a list of species
site_ids <- paste("x", 1:num_sites, sep = "") # Generate site IDs
data <- data.frame(SiteID = character(), # Initialize an empty data frame
                   Species = character(),
                   stringsAsFactors = FALSE)
# Loop through each site and assign 1-20 random species
for (site in site_ids) {
  num_species <- sample(1:20, 1)  # Random number of species for this site
  species <- sample(species_list, num_species, replace = F)  # Random species
  site_data <- data.frame(SiteID = rep(site, num_species), Species = species, stringsAsFactors = FALSE)
  data <- rbind(data, site_data)
}
rm(site_data, num_sites, num_species, site, site_ids, species, species_list) # Cleanup environment

# Result
head(data)


### Function

# This function quantifies RCI scores for unique Sites/Samples based on the observed communities
# One could simply copy and paste this function for use in their own work
# calculate_RCI requires a dataset with unique Site IDs (samples) and species found at each site
# `SiteID` is a character vector of sample names
# `Species` is a character vector of species found at each site
# Data needs to be in long format as in the example dataset `data`

# The function: calculate_RCI()
# Code within the function is notated for clarity 
calculate_RCI <- function(SiteID, Species) { 
  
  Species_Rank <- as.data.frame(table(data$Species)) # Counting incidence per species
  colnames(Species_Rank) <- c("Species", "Count") # Renaming variables
  Species_Rank$Reciprocal = 1 / Species_Rank$Count # Calculate reciprocal of incidence
  unique_vals <- unique(Species_Rank$Reciprocal) # vector of unique reciprocal values used for species rankings
  Species_Rank$Rank <- match(Species_Rank$Reciprocal, sort(unique_vals)) # makes rank values consecutive while ensuring identical values have identical ranks
  Species_Rank <- subset(Species_Rank, select = c(Species, Rank)) # Remove unneeded columns
  Species_Rank <- Species_Rank[order(-Species_Rank$Rank),] # Order dataset (for use later)
  data_RCI <-  merge(data, Species_Rank, by = "Species", all.x = TRUE) # Appending species' rank to the dataset
  RankSum <- aggregate(Rank ~ SiteID, data = data_RCI, FUN = sum) # Calculate the sum of ranks per site
  Richness <- aggregate(Species ~ SiteID, data = data_RCI, FUN = length) # Calculate species richness per site
  colnames(RankSum) <- c("SiteID", "RankSum") # Renaming variables
  colnames(Richness) <- c("SiteID", "Richness") # Renaming variables
  data_RCI <- merge(data_RCI, RankSum, by = "SiteID", all.x = TRUE) # Appending sites' sum of ranks to the dataset
  data_RCI <- merge(data_RCI, Richness, by = "SiteID", all.x = TRUE) # Appending sites' species richness to the dataset
  data_RCI$MaxScore <- sapply(data_RCI$Richness, FUN=function(x){sum(Species_Rank$Rank[1:x])}) # Calculates the maximum sum of ranks for each site given species richness
  data_RCI$RCI <- data_RCI$RankSum / data_RCI$MaxScore # Calculate RCI value per site
  data_RCI <- subset(data_RCI, select = c(SiteID, Richness, RCI)) # Remove unneeded columns (species data)
  data_RCI <- unique(data_RCI) # Remove duplicate rows created by removing species data
  return(data_RCI)
}

# Example for calculate_RCI
data_RCI <- calculate_RCI(data$SiteID, data$Species) # Example for calculate_RCI

# Final Product
head(data_RCI)

# Demonstrating increasing RCI value as site species richness increases
# variation in RCI within species richnesses is the result of varying community rarity
plot(data_RCI$Richness, data_RCI$RCI)