---
Directions and Scripts for Calculating the Chessie BIBI
Last updated: December 8, 2025
Authors:
  Rikke D. Jepsen (rjepsen@icprb.org): Interstate Commission on the Potomac River Basin
  Zachary M. Smith (zachary.smith2@dec.ny.gov): New York State Department of Environmental Conservation
---
**Introduction**  
The following R-scripts and directions are intended as a supplement to the Smith et al. (2017) final report and appendices. These scripts will allow you to calculate the family-level, bioregion-specific Chesapeake Basin-wide Index of Biotic Integrity (Chessie BIBI) for streams and wadeable rivers in the Chesapeake Bay watershed.      

**1. Data Preparation**   
Most macroinvertebrate assemblage data are formatted in a similar manner, but data sets from separate entities often cannot be joined without some manipulation. The data preparation steps below manipulate the data to meet formatting requirements necessary for the BIBI calculation steps.       

A CSV file will be used as the input file to these scripts. Please refer to the "ChessieBIBI CSV setup guide and data_updatedJan2021.xlsx" for information about setting up your CSV file. Data may be entered in uppercase or lowercase, but all text will be converted to lowercase during the data preparation process. The document also contains several keys that explain codes and acronyms that were used. Some of the information that the codes pertain to can be found in domain lists on the bottom of the Water Quality Portal User Guide (https://www.waterqualitydata.us/portal_userguide/). It is important to title all the columns in the same way as they appear in the document, or the scripts will not run as written.               

**1.1** 
Load the packages below into the environment. The packages must be installed before they can be opened. To check if a package is already installed, go to the "Packages"" tab (RStudio's lower right quadrant) and either scroll to find it or type in the name in the search bar. To install a package, either use `install.packages("name")` or the Install button in the Packages tab. To open an installed library, type `library(name)`.    
```{r}
library(tidyverse)
library(stringr)
library(renv)
```
**1.2**
Save this RMarkdown file in a folder where you would like the outputs of this program to be saved and have saved your input CSV file. Use "setwd" to connect RStudio to this folder.

Import the input CSV by first inserting the name of your input CSV where indicated (for example, "Chessie BIBI biology.csv"). 

The next chunk of code will perform some cleanup on the data, including trimmimg unwanted spaces and setting all text to lowercase. Also, the *event_id* and *sample_number* columns will be concatenated in a new column called *unique_id*. This will be the sample identifier in the remainder of the code. 
```{r}
setwd("COPY YOUR WORKING DIRECTORY PATH HERE. CHANGE ALL \ to /") 
data <- read.csv("NAME OF YOUR INPUT CSV.csv") 

prep.df <- data %>%
  rename_all(tolower) %>%
  rename_all(stringr::str_trim) %>%
  mutate_if(is.character, tolower) %>%
  mutate_if(is.character, stringr::str_trim) %>%
  mutate_all(. , .funs=tolower) %>%
  mutate(tsn_final = str_replace_all(tsn_final, "(?<![0-9])0+", ""))
```
**1.3** 
Perform method standardizations as described on page 5 of Smith et al. (2017), with the exception of taxa roll-ups/downs, which have already been completed in the master taxa list that was used to create the input CSV. The following code chunks will perform taxa exclusions, remove samples with macroinvertebrate counts less than 70, exclude samples collected from December through February, include only samples from 1st - 5th order streams, restrict the sampling method to kick-net or similar procedures, and remove samples collected from pools (code BE124P). Last, the Functional Feeding Group (FFG) and habit codes are converted to full names.

2025 Update: In previous analyses, sites identified as Strahler Stream Order (SSO) 5 or higher at a 1:100,000 scale were excluded from the analysis dataset because the index was developed for streams and wadeable rivers identified as 4th order or less at this scale. Map resolution has improved and has led to the delineation of more headwater streams, which has significantly changed SSO designations. We now classify wadeable streams and rivers as SSO 5 and under. If you are calculating Chessie BIBI scores on your own dataset, please use best professional judgement to determine the cutoff for wadeable streams and rivers.  

```{r}
prep.df <- prep.df %>%
  filter(phylum %in% c("annelida", "arthropoda", "mollusca", "platyhelminthes")) %>%
  filter(subphylum %in% c("clitellata", "crustacea", "hexapoda", "rhabditophora") | phylum %in% c("annelida", "arthropoda", "mollusca", "platyhelminthes")) %>%
  filter(!class %in% c("branchiopoda", "maxillopoda", "ostracoda")) %>%
  filter(!order %in% c("hymenoptera")) %>%
  filter(!family %in% c("gerridae", "hebridae", "veliidae", "hydrometridae", "saldidae"))

prep.df <- prep.df %>%
  mutate(reporting_value = as.numeric(reporting_value)) %>%
  mutate(month = as.numeric(month)) %>%
  mutate(aspt = as.numeric(aspt)) %>%
  mutate(bibi_tv = as.numeric(bibi_tv)) %>%
  mutate(month = as.numeric(month)) %>%
  mutate(strahler_stream_order = as.numeric(strahler_stream_order)) %>%
  mutate(g_method = as.numeric(g_method)) %>%
  group_by(unique_id) %>%
  filter(sum(reporting_value) > 70) %>%
  ungroup()

prep.df <- prep.df %>%
  filter(strahler_stream_order <= 5,
         g_method %in% c(49, 57, 58, 59, 86, 87, 89, 92, 93, 94, 101, 103, 104, 106, 107, 108, 110),
         !month %in% c(12, 1, 2),
         !sample_method %in% c("be124p"))

prep.df <- prep.df %>%
  mutate(ffg = case_when(
    is.na(bibi_ffg) | bibi_ffg == "" ~ as.character(NA),
    bibi_ffg == "cg" ~ "gather",
    bibi_ffg == "pr" ~ "predator",
    bibi_ffg == "sc" ~ "scrape",
    bibi_ffg == "cf" ~ "filter",
    bibi_ffg == "sh" ~ "shred",
    bibi_ffg == "om" ~ "omnivore",
    bibi_ffg == "pa" ~ "parasite",
    bibi_ffg %in% c("pc", "pi") ~ "pierce",
    TRUE ~ "ERROR"),
  habit = case_when(
    is.na(bibi_habit) | bibi_habit == "" ~ as.character(NA),
    bibi_habit == "sp" ~ "sprawl",
    bibi_habit == "cn" ~ "cling",
    bibi_habit == "cb" ~ "climb",
    bibi_habit == "bu" ~ "burrow",
    bibi_habit == "sw" ~ "swim",
    bibi_habit == "sk" ~ "skate",
    TRUE ~ "ERROR"))
```
**1.4**
Use `select` to retain only the indicated columns and then sum the counts of repeated taxa in samples. The resulting data frame is **taxa.df**.
```{r}
taxa.df <- prep.df %>% 
  select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, reporting_value, final_id, phylum, subphylum, class, subclass, 
         order, suborder, family, subfamily, tribe, genus, species, habit, ffg, bibi_tv, aspt) %>%
  group_by_at(vars(-reporting_value)) %>%
  summarize(reporting_value = sum(reporting_value)) %>%
  ungroup()
```
**1.5** 
Perform probabilistic rarefaction as described in Smith et al. (2017) page 7 and Appendix D. The Chessie BIBI uses rarefied counts of 100 because that was the lowest common denominator among all the agencies and programs represented in the Chessie BIBI database.  

WARNING: every time you rerun the rarefaction scripts below, the resulting data frame will have different rarefied counts for rare species (usually the singletons or doubletons in your dataset).
```{r}
taxa2.df <- taxa.df %>%                                                                      
  group_by(unique_id) %>%                                                                               
  mutate(sample_sum = sum(reporting_value)) %>%                                                         
  mutate(ratio = 100*(reporting_value/sample_sum)) %>%                                                  
  mutate(ratio_int = ifelse(ratio >= 1, round(ratio), as.integer(ratio))) %>%                           
  mutate(cumulative = ifelse(ratio_int >= 1, sum(ratio_int), 0)) %>%                                    
  mutate(rare_needed = as.integer(ifelse(ratio_int >= 1, 100-cumulative, 100-(max(cumulative))))) %>%   
  mutate(count= n()) %>%
  mutate(random = sample(100, size = n(), replace = FALSE)) %>%                                         
  arrange(ratio_int, random) %>%                                                                        
  ungroup()                                                                                             
 
 taxa3.df <- taxa2.df %>%
   group_by(unique_id) %>%
   mutate(rare_num = case_when(ratio_int == 0 ~ TRUE, ratio_int != 0 ~ FALSE)) %>%                      
   summarise(total_rare = sum(rare_num)) %>%                                                            
   ungroup()
 
 taxa_join <- left_join(taxa2.df, taxa3.df, by = "unique_id") %>%                            
   group_by(unique_id) %>%
   mutate(rare_needed = ifelse(rare_needed > total_rare, total_rare, rare_needed)) %>%                   
   mutate(rare_needed = ifelse(rare_needed < 1, 0, rare_needed)) %>%                                   
   mutate(rare_2b_removed = (total_rare - rare_needed)+1) %>%                                          
   slice(median(rare_2b_removed):n()) %>%                                                                       
   mutate(rare_count =(ifelse(ratio_int == 0, 1, ratio_int))) %>%                                       
   mutate(rare_count_total = sum(rare_count)) %>%                                                        
   ungroup()
   
 taxa.df4 <- taxa_join %>%                                                                               
  select(unique_id, icprb_bioregion_id, rare_count, phylum, subphylum, class, subclass, order, suborder, family, subfamily, tribe, genus, habit, species, ffg, bibi_tv)
```
**2. Calculating the Biological Metrics**  
Now that the data are in the proper format and the reporting values have been rarefied, the biological metrics can be calculated. The metric formulas are located in Smith et al. (2017) Appendix E and are described in the "Metrics used in Chessie BIBI.pdf" on the FTP site. Please note that the appendix contains descriptions of all the metrics tested during the development of the Chessie BIBI, so includes many that did not eventually end up as part of the BIBI. Also, please note that if you downloaded Appendix E before 1/23/2019, you will have to download it again as several typos have been corrected.  

The scripts below compute for each unique sampling event (unique_id) all of the possible family-level metrics needed for the Chessie BIBI. Metrics specific to each bioregion are selected and used to calculate the bioregion-specific BIBI indices in a later step.      

The first step in calculating the metrics is to generate a new data frame for the non-rarefied and rarefied data, both of which are needed because not all metric calculations require rarefied data. Then, a data frame, **metrics.key**, will be generated and metric calculations appended. Note that if your data, and therefore your input CSV, does not include a sample time column, remove all instances of "sample_time" from the scripts or you will get a warning about the field being missing.
```{r}
# Create a data frame of the non-rarefied and rarefied data.
non_rare <- taxa.df %>%
   arrange(unique_id) %>%
   ungroup()

rare <- taxa.df4 %>%
  arrange(unique_id) %>%
  ungroup()

# Create the metrics.key data frame.
metrics.key <- taxa.df %>%              
 select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time) %>%
 unique() %>%
 arrange(unique_id) 

# ESTABLISH FUNCTIONS USED TO CALCULATE SOME METRICS
###--------------------------------------------------------------------------
original_order <- function(mod.df, original.df, unique.id.col) {
  unique.id.col <- rlang::enquo(unique.id.col)

  final.df <- original.df %>%
    dplyr::select(!!unique.id.col) %>%
    dplyr::distinct() %>%
    dplyr::left_join(mod.df, by = rlang::quo_name(unique.id.col))

  return(final.df)
}
###--------------------------------------------------------------------------
## CALCULATE THE METRICS ##

sample_sum <- non_rare %>%
 group_by(unique_id) %>%
 summarise(sample_sum = sum(reporting_value)) 

#
 aspt_mod <- function(long.df, unique.id.col, count.col, tol.col) {
  unique.id.col <- rlang::enquo(unique.id.col)
  tol.col <- rlang::enquo(tol.col )
  count.col <- rlang::enquo(count.col )
  #----------------------------------------------------------------------------
  long.df <- long.df %>%
    dplyr::filter((!!count.col) > 0)
  #----------------------------------------------------------------------------
 score.vec <- long.df %>%
    filter(!is.na(!!!tol.col),
           !rlang::UQ(tol.col) %in% c("")) %>%
    select(!!unique.id.col, !!tol.col, !!count.col) %>%
    mutate(score = rlang::UQ(count.col) * rlang::UQ(tol.col)) %>%
    group_by(!!unique.id.col) %>%
    summarize(score = sum(score),
              taxa = sum(!!count.col),
              score = score / taxa) %>%
    original_order(long.df, !!unique.id.col) %>%
    pull(score)

  return(score.vec)
}

metrics.key$aspt_mod <- aspt_mod(non_rare, unique_id, reporting_value, aspt)

# pct_gastro_oligo 
pct_gastro_oligo <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(class)) %>% 
 summarise(sum = sum(reporting_value[class %in% c("gastropoda", "clitellata")]))  
 metrics.key$pct_gastro_oligo <- (pct_gastro_oligo$sum / sample_sum$sample_sum )*100

# pct_diptera
pct_diptera <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>% 
 summarise(sum = sum(reporting_value[order %in% c("diptera")]))  
 metrics.key$pct_diptera <- (pct_diptera$sum / sample_sum$sample_sum)*100

# GOLD  
gold_GOL <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(class)) %>% 
 summarise(sum_GOL = sum(reporting_value[class %in% c("gastropoda", "clitellata")]))

gold_Dip <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>% 
 summarise(sum_Dip = sum(reporting_value[order %in% c("diptera")]))  

metrics.key$gold <- (1-(gold_GOL$sum_GOL + gold_Dip$sum_Dip)/sample_sum$sample_sum)*100

# HBI    
hbi_multiplier <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>% 
 summarise(multiplier = sum(reporting_value*bibi_tv)) 

hbi_divisor <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>% 
 summarise(divisor = sum(reporting_value))

metrics.key$hbi <- hbi_multiplier$multiplier / hbi_divisor$divisor

# pct_arthropoda      
pct_arthropoda <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(phylum)) %>% 
 summarise(sum = sum(reporting_value[phylum %in% c("arthropoda")]))  
 metrics.key$pct_arthropoda <- (pct_arthropoda$sum / sample_sum$sample_sum )*100

# pct_burrow    
pct_burrow <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(habit)) %>% 
 summarise(sum = sum(reporting_value[habit %in% c("burrow")]))  
 metrics.key$pct_burrow <- (pct_burrow$sum / sample_sum$sample_sum )*100
 
# pct_chironomidae   
pct_chironomidae <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(family)) %>% 
 summarise(sum = sum(reporting_value[family %in% c("chironomidae")]))  
 metrics.key$pct_chironomidae <- (pct_chironomidae$sum / sample_sum$sample_sum )*100

# pct_cling  
pct_cling <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(habit)) %>% 
 summarise(sum = sum(reporting_value[habit %in% c("cling")]))  
 metrics.key$pct_cling <- (pct_cling$sum / sample_sum$sample_sum )*100

#pct_collect 
pct_collect <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>% 
 summarise(sum = sum(reporting_value[ffg %in% c("filter", "gather")]))  
 metrics.key$pct_collect <- (pct_collect$sum / sample_sum$sample_sum )*100 

# pct_cote    
pct_cote <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>% 
 summarise(sum = sum(reporting_value[order %in% c("coleoptera", "trichoptera", "ephemeroptera", "odonata")]))  
 metrics.key$pct_cote <- (pct_cote$sum / sample_sum$sample_sum )*100

# pct_ephemeroptera    
pct_ephemeroptera <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>% 
 summarise(sum = sum(reporting_value[order %in% c("ephemeroptera")]))  
 metrics.key$pct_ephemeroptera <- (pct_ephemeroptera$sum / sample_sum$sample_sum )*100

# pct_baetidae    
pct_baetidae <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(family)) %>% 
 summarise(sum = sum(reporting_value[family %in% c("baetidae")]))  
 metrics.key$pct_baetidae <- (pct_baetidae$sum / sample_sum$sample_sum )*100
 
# pct_ephemeroptera_no_baetid       
pct_ephemeroptera_no_baetid <- pct_baetidae
pct_ephemeroptera_no_baetid$sum_ephemeroptera <- pct_ephemeroptera$sum 
pct_ephemeroptera_no_baetid$ephem_no_baetidae <- pct_ephemeroptera_no_baetid$sum_ephemeroptera - pct_ephemeroptera_no_baetid$sum
metrics.key$pct_ephemeroptera_no_baetid <- (pct_ephemeroptera_no_baetid$ephem_no_baetidae/sample_sum$sample_sum)*100

# pct_ept  
pct_ept <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>% 
 summarise(sum = sum(reporting_value[order %in% c("ephemeroptera", "plecoptera", "trichoptera")]))  
 metrics.key$pct_ept <- (pct_ept$sum / sample_sum$sample_sum )*100
  
# pct_hydropsychidae    
pct_hydropsychidae <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(family)) %>% 
 summarise(sum = sum(reporting_value[family %in% c("hydropsychidae")]))  
 metrics.key$pct_hydropsychidae <- (pct_hydropsychidae$sum / sample_sum$sample_sum )*100

# pct_ept_no_hydro    
pct_ept_hydro <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>% 
 summarise(sum = sum(reporting_value[family %in% c("hydropsychidae")]))
 pct_ept_hydro$subtract <- pct_ept$sum - pct_ept_hydro$sum
 metrics.key$pct_ept_no_hydro <- (pct_ept_hydro$subtract/pct_ept$sum)*100

# pct_euholognatha  
pct_euholognatha <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(suborder)) %>% 
 summarise(sum = sum(reporting_value[suborder %in% c("euholognatha")]))
 metrics.key$pct_euholognatha <- (pct_euholognatha$sum / sample_sum$sample_sum )*100
  
# pct_filter    
pct_filter <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>% 
 summarise(sum = sum(reporting_value[ffg %in% c("filter")]))  
 metrics.key$pct_filter <- (pct_filter$sum / sample_sum$sample_sum )*100

# pct_heptageniidae    
pct_heptageniidae <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(family)) %>% 
 summarise(sum = sum(reporting_value[family %in% c("heptageniidae")]))  
 metrics.key$pct_heptageniidae <- (pct_heptageniidae$sum / sample_sum$sample_sum )*100

# pct_hexapoda    
pct_hexapoda <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(subphylum)) %>% 
 summarise(sum = sum(reporting_value[subphylum %in% c("hexapoda")]))  
 metrics.key$pct_hexapoda <- (pct_hexapoda$sum / sample_sum$sample_sum )*100

# pct_hydro_ept    
metrics.key$pct_hydro_ept <- ifelse(metrics.key$pct_ept == 0, 0, (metrics.key$pct_hydropsychidae / metrics.key$pct_ept) * 100)

# pct_intol_0_3    
pct_intol_0_3 <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>% 
 summarise(sum = sum(reporting_value[bibi_tv %in% c("0", "1", "2", "3")]))  
 metrics.key$pct_intol_0_3 <- (pct_intol_0_3$sum / sample_sum$sample_sum)*100
  
# pct_intol_0_4   
pct_intol_0_4 <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>% 
 summarise(sum = sum(reporting_value[bibi_tv %in% c("0", "1", "2", "3", "4")]))  
 metrics.key$pct_intol_0_4 <- (pct_intol_0_4$sum / sample_sum$sample_sum)*100

# pct_amph_iso  
pct_amph_iso <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>% 
 summarise(sum = sum(reporting_value[order %in% c("amphipoda", "isopoda")]))  
 metrics.key$pct_amph_iso <- (pct_amph_iso$sum / sample_sum$sample_sum )*100  

# pct_ephemerellidae  
pct_ephemerellidae <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(family)) %>% 
 summarise(sum = sum(reporting_value[family %in% c("ephemerellidae")]))  
 metrics.key$pct_ephemerellidae <- (pct_ephemerellidae$sum / sample_sum$sample_sum )*100 

# pct_limestone    
metrics.key$pct_limestone <- metrics.key$pct_amph_iso + metrics.key$pct_ephemerellidae

# pct_mod_tol_4_6 
pct_mod_tol_4_6 <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>% 
 summarise(sum = sum(reporting_value[bibi_tv %in% c("4", "5", "6")]))  
 metrics.key$pct_mod_tol_4_6 <- (pct_mod_tol_4_6$sum / sample_sum$sample_sum)*100

# pct_trichoptera    
pct_trichoptera <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>% 
 summarise(sum = sum(reporting_value[order %in% c("trichoptera")]))  
 metrics.key$pct_trichoptera <- (pct_trichoptera$sum / sample_sum$sample_sum )*100

# pct_non_hydrop_trichoptera  
pct_non_hydrop_trichoptera <- pct_trichoptera
pct_non_hydrop_trichoptera$hydropsychidae <- pct_hydropsychidae$sum
pct_non_hydrop_trichoptera$subtract <- pct_non_hydrop_trichoptera$sum - pct_non_hydrop_trichoptera$hydropsychidae
metrics.key$pct_non_hydrop_trichoptera <- (pct_non_hydrop_trichoptera$subtract/pct_non_hydrop_trichoptera$sum)*100

# pct_odonata  
pct_odonata <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>% 
 summarise(sum = sum(reporting_value[order %in% c("odonata")]))  
 metrics.key$pct_odonata <- (pct_odonata$sum / sample_sum$sample_sum )*100
  
# pct_pisciforma  
pct_pisciforma <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(suborder)) %>% 
 summarise(sum = sum(reporting_value[suborder %in% c("pisciforma")]))  
 metrics.key$pct_pisciforma <- (pct_pisciforma$sum / sample_sum$sample_sum )*100

# pct_predator  
# In next interation of Chessie BIBI calculations, "piercer" will be added to %predator calculations
pct_predator <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>% 
 summarise(sum = sum(reporting_value[ffg %in% c("predator")]))  
 metrics.key$pct_predator <- (pct_predator$sum / sample_sum$sample_sum )*100 
  
# pct_pterygota  
pct_pterygota <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(subclass)) %>% 
 summarise(sum = sum(reporting_value[subclass %in% c("pterygota")]))  
 metrics.key$pct_pterygota <- (pct_pterygota$sum / sample_sum$sample_sum )*100

# pct_retreat_caddisfly  
pct_retreat_caddisfly <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(suborder)) %>% 
 summarise(sum = sum(reporting_value[suborder %in% c("annulipalpia")]))  
 metrics.key$pct_retreat_caddisfly <- (pct_retreat_caddisfly$sum / sample_sum$sample_sum )*100

# pct_scrape   
pct_scrape <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>% 
 summarise(sum = sum(reporting_value[ffg %in% c("scrape")]))  
 metrics.key$pct_scrape <- (pct_scrape$sum / sample_sum$sample_sum )*100
   
# pct_sprawl   
pct_sprawl <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(habit)) %>% 
 summarise(sum = sum(reporting_value[habit %in% c("sprawl")]))  
 metrics.key$pct_sprawl <- (pct_sprawl$sum / sample_sum$sample_sum )*100

# pct_systellognatha   
pct_systellognatha <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(suborder)) %>% 
 summarise(sum = sum(reporting_value[suborder %in% c("systellognatha")]))  
 metrics.key$pct_systellognatha <- (pct_systellognatha$sum / sample_sum$sample_sum )*100

# pct_tolerant_7_10   
pct_tolerant_7_10 <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>% 
 summarise(sum = sum(reporting_value[bibi_tv %in% c("7", "8", "9", "10")]))  
 metrics.key$pct_tolerant_7_10 <- (pct_tolerant_7_10$sum / sample_sum$sample_sum)*100

# pct_urban_intol   
pct_urban_intol <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(family)) %>% 
 summarise(sum = sum(reporting_value[family %in% c(
    "aeshnidae",
    "ameletidae",
    "asellidae",
    "athericidae",
    "baetidae",
    "brachycentridae",
    "caenidae",
    "calamoceratidae",
    "cambaridae",
    "capniidae",
    "ceratopogonidae",
    "chironomidae",
    "chloroperlidae",
    "cordulegastridae",
    "corduliidae",
    "corydalidae",
    "crangonyctidae",
    "elmidae",
    "ephemerellidae",
    "ephemeridae",
    "glossosomatidae",
    "gomphidae",
    "heptageniidae",
    "hydropsychidae",
    "hydroptilidae",
    "isonychiidae",
    "lepidostomatidae",
    "leptophlebiidae",
    "leuctridae",
    "libellulidae",
    "metretopodidae",
    "nemouridae",
    "odontoceridae",
    "peltoperlidae",
    "perlidae",
    "perlodidae",
    "philopotamidae",
    "phryganeidae",
    "polycentropodidae",
    "polymitarcyidae",
    "potamanthidae",
    "psephenidae",
    "pteronarcyidae",
    "rhyacophilidae",
    "sericostomatidae",
    "sialidae",
    "simuliidae",
    "stratiomyidae",
    "tabanidae",
    "taeniopterygidae",
    "tipulidae",
    "uenoidae",
    "viviparidae")]))  
  metrics.key$pct_urban_intol <- (pct_urban_intol$sum / sample_sum$sample_sum )*100

# Richness and diversity metrics require rarefied data
# Calculate family-level taxa richness for each unique_id
family_richness <- rare %>%
 group_by(unique_id) %>%
 summarise(family_richness = n_distinct(family)) 

# Margalefs
margalefs_numerator <- family_richness %>%
 mutate(family_richness = family_richness - 1)

margalefs_denominator <- rare %>%
 group_by(unique_id) %>%
 summarise(sample_sum = log(sum(rare_count)))

metrics.key$margalefs <- margalefs_numerator$family_richness / margalefs_denominator$sample_sum

# pct_ept_rich
pct_ept_rich <- rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>%
 filter(order %in% c("ephemeroptera", "plecoptera", "trichoptera")) %>%
 summarise(rich_ept = n_distinct(family))
 pct_ept_rich2 <- merge(pct_ept_rich, family_richness, by= "unique_id", all=TRUE)
 pct_ept_rich2[is.na(pct_ept_rich2)] <- 0
 metrics.key$pct_ept_rich <- (pct_ept_rich2$rich_ept / pct_ept_rich2$family_richness)*100
    
# pct_ept_rich_no_tol    
pct_ept_rich_no_tol <- rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>%
 filter(order %in% c("ephemeroptera", "plecoptera", "trichoptera")) %>%
 filter(bibi_tv %in% c("0","1","2","3","4","5","6")) %>%
 summarise(rich_ept_no_tol = n_distinct(family))
 pct_ept_rich_no_tol2 <- merge(pct_ept_rich_no_tol, family_richness, by= "unique_id", all=TRUE)
 pct_ept_rich_no_tol2[is.na(pct_ept_rich_no_tol2)] <- 0
 metrics.key$pct_ept_rich_no_tol <- (pct_ept_rich_no_tol2$rich_ept_no_tol / pct_ept_rich_no_tol2$family_richness)*100

# Pielou
family_richness_pielou <- rare %>%
  group_by(unique_id) %>%
  summarise(family_richness = n_distinct(family)) %>%
  mutate(ln_richness = log(family_richness)) # in R, log is natural log

sample_sum_rare <- rare %>%
  group_by(unique_id) %>%
  summarise(sample_sum = sum(rare_count)) 
sample_sum_join <- sample_sum_rare %>%
  left_join(rare, sample_sum_rare, by = "unique_id") %>%
  select(unique_id, sample_sum)

pielou_shannon <- rare %>%
  group_by(unique_id, family) %>%
  summarise(taxa_sum = sum(rare_count)) %>%
  left_join(sample_sum_join, pielou_shannon, by = "unique_id") %>%
  distinct() %>%
  mutate(Pi = taxa_sum / sample_sum) %>%
  mutate(Pi_lnPi = -(Pi * log(Pi))) %>%
  summarise(Shannon = sum(Pi_lnPi))

metrics.key$pielou <- pielou_shannon$Shannon / family_richness_pielou$ln_richness

# rich_burrow  
rich_burrow <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(habit)) %>%
 filter(habit %in% c("burrow")) %>%
 summarise(rich_burrow = n_distinct(family)) %>%
 select(unique_id, rich_burrow) 
 metrics.key <- merge(metrics.key, rich_burrow, by="unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0 
  
# rich_climb
rich_climb <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(habit)) %>%
 filter(habit %in% c("climb")) %>%
 summarise(rich_climb = n_distinct(family)) %>%
 select(unique_id, rich_climb) 
 metrics.key <- merge(metrics.key, rich_climb, by="unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0 
  
# rich_cling 
rich_cling <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(habit)) %>%
 filter(habit %in% c("cling")) %>%
 summarise(rich_cling = n_distinct(family)) %>%
 select(unique_id, rich_cling) 
 metrics.key <- merge(metrics.key, rich_cling, by="unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0 
 
# rich_collect   
rich_collect_filter <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>%
 filter(ffg %in% c("filter")) %>%
 summarise(rich_collect_filter = n_distinct(family)) %>%
 select(unique_id, rich_collect_filter)
 
rich_collect_gather <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>%
 filter(ffg %in% c("gather")) %>%
 summarise(rich_collect_gather = n_distinct(family)) %>%
 select(unique_id, rich_collect_gather)
 rich_collect <- merge(rich_collect_filter, rich_collect_gather, by = "unique_id", all=TRUE)
 rich_collect[is.na(rich_collect)] <- 0  
 rich_collect$rich_collect <- rich_collect$rich_collect_filter + rich_collect$rich_collect_gather
 rich_collect$rich_collect_gather <- NULL
 rich_collect$rich_collect_filter <- NULL
 metrics.key <- merge(metrics.key, rich_collect, by = "unique_id", all=TRUE)
   
# rich_ephemeroptera     
rich_ephemeroptera <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>%
 filter(order %in% c("ephemeroptera")) %>%
 summarise(rich_ephemeroptera = n_distinct(family)) %>%
 select(unique_id, rich_ephemeroptera) 
 metrics.key <- merge(metrics.key, rich_ephemeroptera, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0
     
# rich_ept    
rich_ept <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>%
 filter(order %in% c("ephemeroptera", "plecoptera", "trichoptera")) %>%
 summarise(rich_ept = n_distinct(family)) %>%
 select(unique_id, rich_ept) 
 metrics.key <- merge(metrics.key, rich_ept, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0   
 
# rich_filter    
rich_filter <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>%
 filter(ffg %in% c("filter")) %>%
 summarise(rich_filter = n_distinct(family)) %>%
 select(unique_id, rich_filter) 
 metrics.key <- merge(metrics.key, rich_filter, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0   

# rich_gather    
rich_gather <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>%
 filter(ffg %in% c("gather")) %>%
 summarise(rich_gather = n_distinct(family)) %>%
 select(unique_id, rich_gather) 
 metrics.key <- merge(metrics.key, rich_gather, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0   

# rich_intol 
rich_intol0 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("0")) %>%
 summarise(rich_intol0 = n_distinct(family)) %>%
 select(unique_id, rich_intol0)

rich_intol1 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("1")) %>%
 summarise(rich_intol1 = n_distinct(family)) %>%
 select(unique_id, rich_intol1)

rich_intol2 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("2")) %>%
 summarise(rich_intol2 = n_distinct(family)) %>%
 select(unique_id, rich_intol2)

rich_intol3 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("3")) %>%
 summarise(rich_intol3 = n_distinct(family)) %>%
 select(unique_id, rich_intol3)
 rich_intol01 <- merge(rich_intol0, rich_intol1, by= "unique_id", all=TRUE)
 rich_intol01[is.na(rich_intol01)] <- 0
 rich_intol012 <- merge(rich_intol01, rich_intol2, by= "unique_id", all=TRUE)
 rich_intol012[is.na(rich_intol012)] <- 0
 rich_intol0123 <- merge(rich_intol012, rich_intol3, by= "unique_id", all=TRUE)
 rich_intol0123[is.na(rich_intol0123)] <- 0
 rich_intol0123$rich_intol <- rich_intol0123$rich_intol0 + rich_intol0123$rich_intol1 + rich_intol0123$rich_intol2 + rich_intol0123$rich_intol3
 
rich_intol0123 <- rich_intol0123 %>%
 select(unique_id, rich_intol)
 metrics.key <- merge(metrics.key, rich_intol0123, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0   

# rich_modtol  
rich_modtol4 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("4")) %>%
 summarise(rich_modtol4 = n_distinct(family)) %>%
 select(unique_id, rich_modtol4)

rich_modtol5 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("5")) %>%
 summarise(rich_modtol5 = n_distinct(family)) %>%
 select(unique_id, rich_modtol5)

rich_modtol6 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("6")) %>%
 summarise(rich_modtol6 = n_distinct(family)) %>%
 select(unique_id, rich_modtol6)
 rich_modtol45 <- merge(rich_modtol4, rich_modtol5, by= "unique_id", all=TRUE)
 rich_modtol45[is.na(rich_modtol45)] <- 0
 rich_modtol456 <- merge(rich_modtol45, rich_modtol6, by= "unique_id", all=TRUE)
 rich_modtol456[is.na(rich_modtol456)] <- 0
 rich_modtol456$rich_modtol <- rich_modtol456$rich_modtol4 + rich_modtol456$rich_modtol5 + rich_modtol456$rich_modtol6
 
 rich_modtol456 <- rich_modtol456 %>%
 select(unique_id, rich_modtol)
 metrics.key <- merge(metrics.key, rich_modtol456, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0  
  
# rich_plecoptera  
rich_plecoptera <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>%
 filter(order %in% c("plecoptera")) %>%
 summarise(rich_plecoptera = n_distinct(family)) %>%
 select(unique_id, rich_plecoptera) 
 metrics.key <- merge(metrics.key, rich_plecoptera, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0  

# rich_predator    
rich_predator <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>%
 filter(ffg %in% c("predator")) %>%
 summarise(rich_predator = n_distinct(family)) %>%
 select(unique_id, rich_predator) 
 metrics.key <- merge(metrics.key, rich_predator, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0 
 
# rich_shred    
rich_shred <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>%
 filter(ffg %in% c("shred")) %>%
 summarise(rich_shred = n_distinct(family)) %>%
 select(unique_id, rich_shred) 
 metrics.key <- merge(metrics.key, rich_shred, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0 

# rich_sprawl    
rich_sprawl <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(habit)) %>%
 filter(habit %in% c("sprawl")) %>%
 summarise(rich_sprawl = n_distinct(family)) %>%
 select(unique_id, rich_sprawl) 
 metrics.key <- merge(metrics.key, rich_sprawl, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0 
    
# rich_tol    
rich_tol7 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("7")) %>%
 summarise(rich_tol7 = n_distinct(family)) %>%
 select(unique_id, rich_tol7)

rich_tol8 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("8")) %>%
 summarise(rich_tol8 = n_distinct(family)) %>%
 select(unique_id, rich_tol8)

rich_tol9 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("9")) %>%
 summarise(rich_tol9 = n_distinct(family)) %>%
 select(unique_id, rich_tol9)

rich_tol10 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("10")) %>%
 summarise(rich_tol10 = n_distinct(family)) %>%
 select(unique_id, rich_tol10)
 rich_tol78 <- merge(rich_tol7, rich_tol8, by= "unique_id", all=TRUE)
 rich_tol78[is.na(rich_tol78)] <- 0
 rich_tol789 <- merge(rich_tol78, rich_tol9, by= "unique_id", all=TRUE)
 rich_tol789[is.na(rich_tol789)] <- 0
 rich_tol78910 <- merge(rich_tol789, rich_tol10, by= "unique_id", all=TRUE)
 rich_tol78910[is.na(rich_tol78910)] <- 0
 rich_tol78910$rich_tol <- rich_tol78910$rich_tol7 + rich_tol78910$rich_tol8 + rich_tol78910$rich_tol9 + rich_tol78910$rich_tol10
 
rich_tol78910 <- rich_tol78910 %>%
 select(unique_id, rich_tol)
 metrics.key <- merge(metrics.key, rich_tol78910, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0   
    
# rich_trichoptera  
rich_trichoptera <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>%
 filter(order %in% c("trichoptera")) %>%
 summarise(rich_trichoptera = n_distinct(family)) %>%
 select(unique_id, rich_trichoptera) 
 metrics.key <- merge(metrics.key, rich_trichoptera, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0 
```
If you would like to export **metrics.key** as a CSV, run the following script to save the file to your working directory. Don't forget to replace the text with your own where indicated.
```{r}
setwd("WORKING DIRECTORY")
write.csv(metrics.key, file= "FILE NAME.csv")
```
**3. Metric Scoring**  
**3.1**
Subset **metrics.key** into twelve dataframes, one for each bioregion.
```{r}
NAPU <- metrics.key %>% filter(icprb_bioregion_id == "napu")
NCA <- metrics.key %>% filter(icprb_bioregion_id == "nca")
CA <- metrics.key %>% filter(icprb_bioregion_id == "ca")
MAC<- metrics.key %>% filter(icprb_bioregion_id == "mac")
SEP <- metrics.key %>% filter(icprb_bioregion_id == "sep")
BLUE <- metrics.key %>% filter(icprb_bioregion_id == "blue")
NRV <- metrics.key %>% filter(icprb_bioregion_id == "nrv")
SRV <- metrics.key %>% filter(icprb_bioregion_id == "srv")
UNP <- metrics.key %>% filter(icprb_bioregion_id == "unp")
SGV <- metrics.key %>% filter(icprb_bioregion_id == "sgv")
LNP <- metrics.key %>% filter(icprb_bioregion_id == "lnp")
PIED <- metrics.key %>% filter(icprb_bioregion_id == "pied")
```
**3.2**
Subset each of the twelve bioregion-specific dataframes so that they contain only the metrics specific to that bioregion. A list of the metrics used for each bioregion can be found in Table 15, Page 44 of Smith et al. (2017). One change has been made from the original list of metrics: pct_predator was removed from NAPU, as it was later deemed unnecessary. 
```{r}
#Northern Appalachian Plateau and Uplands (NAPU)
NAPU_score <- NAPU %>% 
  select(
  "unique_id",
  "icprb_bioregion_id",
  "station_id",
  "sample_date",
  "sample_time",
  "margalefs", 
  "pct_cling",
  "pct_collect",
  "pct_ephemeroptera",
  "pct_ept",
  "pct_ept_no_hydro",
  "pct_ept_rich_no_tol",
  "pct_filter",
  "pct_intol_0_3",
  "pct_non_hydrop_trichoptera",
  "pct_scrape",
  "pielou",
  "rich_burrow",
  "rich_cling",
  "rich_gather",
  "rich_modtol",
  "rich_plecoptera",
  "rich_shred",
  "rich_sprawl")
       
# North Central Appalachians (NCA)
NCA_score <- NCA %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "aspt_mod",
    "pct_cote",
    "pct_ephemeroptera",
    "pct_heptageniidae",
    "pct_limestone",
    "rich_cling",
    "rich_ephemeroptera",
    "rich_ept",
    "rich_plecoptera")

# Central Appalachians (CA)
CA_score <- CA %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "pct_ept_rich",
    "pct_euholognatha",
    "pct_heptageniidae",
    "pct_hydro_ept",
    "pct_mod_tol_4_6",
    "pct_non_hydrop_trichoptera",
    "pct_pisciforma",
    "pct_systellognatha",
    "rich_burrow",
    "rich_collect",
    "rich_filter",
    "rich_plecoptera",
    "rich_predator")

# Middle Atlantic Coastal Plain (MAC)
MAC_score <- MAC %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "aspt_mod",
    "gold",
    "pct_arthropoda",
    "pct_chironomidae",
    "pct_collect",
    "pct_mod_tol_4_6",
    "pct_odonata",
    "pct_predator",
    "pct_scrape",
    "pct_urban_intol",
    "rich_climb",
    "rich_tol")

# Southeastern Plains (SEP)
SEP_score <- SEP %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "aspt_mod",
    "gold",
    "hbi",
    "pct_cling",
    "pct_cote",
    "pct_ept_rich_no_tol",
    "rich_cling",
    "rich_collect",
    "rich_ephemeroptera",
    "rich_filter",
    "rich_trichoptera")

# Blue Ridge (BLUE)
BLUE_score <- BLUE %>% 
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "pct_burrow",
    "pct_ept_rich_no_tol",
    "pct_scrape",
    "pct_systellognatha",
    "rich_ephemeroptera")

# Northern Ridge and Valley (NRV)
NRV_score <- NRV %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "aspt_mod",
    "pct_ephemeroptera_no_baetid",
    "pct_ept_rich",
    "rich_collect",
    "rich_ephemeroptera",
    "rich_predator",
    "rich_shred",
    "rich_trichoptera")

# Southern Ridge and Valley (SRV)
SRV_score <- SRV %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "aspt_mod",
    "pct_burrow",
    "pct_cling",
    "pct_cote",
    "pct_ephemeroptera_no_baetid",
    "pct_ept_rich_no_tol",
    "pct_predator",
    "pct_scrape",
    "rich_cling",
    "rich_collect",
    "rich_ephemeroptera",
    "rich_intol",
    "rich_plecoptera",
    "rich_predator",
    "rich_trichoptera")

# Upper Northern Piedmont (UNP)
UNP_score <- UNP %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "hbi",
    "pct_collect",
    "pct_ept_rich_no_tol",
    "pct_filter",
    "pct_hexapoda",
    "pct_pisciforma",
    "pct_retreat_caddisfly",
    "pct_scrape",
    "pct_urban_intol",
    "rich_ephemeroptera",
    "rich_plecoptera",
    "rich_trichoptera")

# Southern Great Valley (SGV)
SGV_score <- SGV %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "aspt_mod",
    "gold",
    "pct_cling",
    "pct_collect",
    "pct_ephemeroptera_no_baetid",
    "pct_ept_no_hydro",
    "pct_ept_rich_no_tol",
    "pct_filter",
    "pct_heptageniidae",
    "pct_mod_tol_4_6",
    "pct_non_hydrop_trichoptera",
    "pct_pisciforma",
    "pct_predator",
    "pct_pterygota",
    "pct_scrape",
    "pct_sprawl",
    "pct_tolerant_7_10",
    "rich_collect",
    "rich_ephemeroptera",
    "rich_filter",
    "rich_trichoptera")

# Lower Northern Piedmont (LNP)
LNP_score <- LNP %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "aspt_mod",
    "pct_cling",
    "pct_collect",
    "pct_ept_rich",
    "pct_mod_tol_4_6",
    "pct_predator",
    "pct_pterygota",
    "rich_cling",
    "rich_ephemeroptera")

# Piedmont (PIED)
PIED_score <- PIED %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "pct_cling",
    "pct_ept_rich",
    "pct_intol_0_4",
    "pct_tolerant_7_10",
    "pct_urban_intol",
    "rich_cling",
    "rich_ephemeroptera",
    "rich_trichoptera")
```
####**3.3**
Score the metrics using the approach described on pages 20-21 in Smith et al. (2017).
```{r}
# NAPU Metric Scoring
NAPU_score$score_margalefs <- 
  case_when(
    NAPU_score$margalefs > 2.83 ~ 100,
    NAPU_score$margalefs < 2.74 ~ 0,
    TRUE ~ ((NAPU_score$margalefs - 2.74)/(2.83-2.74))*100)

NAPU_score$score_pct_cling <- 
  case_when(
    NAPU_score$pct_cling > 59.19 ~ 100,
    NAPU_score$pct_cling < 50.98 ~ 0,
    TRUE ~ ((NAPU_score$pct_cling-50.98)/(59.19-50.98))*100)

NAPU_score$score_pct_collect <- 
   case_when(
     NAPU_score$pct_collect < 71.84 ~ 100,
     NAPU_score$pct_collect > 82.8 ~ 0,
     TRUE ~ ((82.8-NAPU_score$pct_collect)/(82.8-71.84))*100)

NAPU_score$score_pct_ephemeroptera <- 
  case_when(
    NAPU_score$pct_ephemeroptera > 31.47 ~ 100,
    NAPU_score$pct_ephemeroptera < 19.62 ~ 0,
    TRUE ~ ((NAPU_score$pct_ephemeroptera-19.62)/(31.47-19.62))*100)

NAPU_score$score_pct_ept <- 
  case_when(
    NAPU_score$pct_ept > 62.26 ~ 100,
    NAPU_score$pct_ept < 54.08 ~ 0,
    TRUE ~ ((NAPU_score$pct_ept-54.08)/(62.26-54.08))*100)

NAPU_score$score_pct_ept_no_hydro <- 
  case_when(
    NAPU_score$pct_ept_no_hydro > 87.96 ~ 100,
    NAPU_score$pct_ept_no_hydro < 70.51 ~ 0,
    TRUE ~ ((NAPU_score$pct_ept_no_hydro-70.51)/(87.96-70.51))*100)

NAPU_score$score_pct_ept_rich_no_tol <- 
  case_when(
    NAPU_score$pct_ept_rich_no_tol > 45.32 ~ 100,
    NAPU_score$pct_ept_rich_no_tol < 32.85 ~ 0,
    TRUE ~ ((NAPU_score$pct_ept_rich_no_tol-32.85)/(45.32-32.85))*100)

NAPU_score$score_pct_filter <- 
  case_when(
    NAPU_score$pct_filter < 15.51 ~ 100,
    NAPU_score$pct_filter > 22.62 ~ 0,
    TRUE ~ ((22.62-NAPU_score$pct_filter)/(22.62-15.51))*100)

NAPU_score$score_pct_intol_0_3 <- 
  case_when(
    NAPU_score$pct_intol_0_3 > 44.7 ~ 100,
    NAPU_score$pct_intol_0_3 < 24.22 ~ 0,
    TRUE ~ ((NAPU_score$pct_intol_0_3-24.22)/(44.7-24.22))*100)

NAPU_score$score_pct_non_hydrop_trichoptera <-
  case_when(
    NAPU_score$pct_non_hydrop_trichoptera > 45.08 ~ 100,
    NAPU_score$pct_non_hydrop_trichoptera < 20.6 ~ 0,
    TRUE ~ ((NAPU_score$pct_non_hydrop_trichoptera-20.6)/(45.08-20.6))*100)

NAPU_score$score_pct_scrape <- 
  case_when(
    NAPU_score$pct_scrape > 9.36 ~ 100,
    NAPU_score$pct_scrape < 6.89 ~ 0,
    TRUE ~ ((NAPU_score$pct_scrape-6.89)/(9.36-6.89))*100)

NAPU_score$score_pielou <- 
  case_when(
    NAPU_score$pielou < 0.79 ~ 100,
    NAPU_score$pielou > 0.80 ~ 0,
    TRUE ~ ((0.80-NAPU_score$pielou)/(0.80-0.79))*100)

NAPU_score$score_rich_burrow <- 
  case_when(
    NAPU_score$rich_burrow > 2.95 ~ 100,
    NAPU_score$rich_burrow < 1.93 ~ 0,
    TRUE ~ ((NAPU_score$rich_burrow-1.93)/(2.95-1.93))*100)

NAPU_score$score_rich_cling <- 
  case_when(
    NAPU_score$rich_cling > 8.0 ~ 100,
    NAPU_score$rich_cling < 6.77 ~ 0,
    TRUE ~ ((NAPU_score$rich_cling-6.77)/(8.0-6.77))*100)

NAPU_score$score_rich_gather <- 
  case_when(
    NAPU_score$rich_gather < 5.0 ~ 100,
    NAPU_score$rich_gather > 6.01 ~ 0,
    TRUE ~ ((6.01-NAPU_score$rich_gather)/(6.01-5.0))*100)

NAPU_score$score_rich_modtol <- 
  case_when(
     NAPU_score$rich_modtol < 6.0 ~ 100,
     NAPU_score$rich_modtol > 7.03 ~ 0,
     TRUE ~ ((7.03-NAPU_score$rich_modtol)/(7.03-6.0))*100)

NAPU_score$score_rich_plecoptera <- 
  case_when(
    NAPU_score$rich_plecoptera > 2.13 ~ 100,
    NAPU_score$rich_plecoptera < 0.70 ~ 0,
    TRUE ~ ((NAPU_score$rich_plecoptera-0.70)/(2.13-0.70))*100)

NAPU_score$score_rich_shred <- 
  case_when(
    NAPU_score$rich_shred > 2.0 ~ 100,
    NAPU_score$rich_shred < 0.84 ~ 0,
    TRUE ~ ((NAPU_score$rich_shred-0.84)/(2.0-0.84))*100)

NAPU_score$score_rich_sprawl <- 
  case_when(
    NAPU_score$rich_sprawl < 1.0 ~ 100,
    NAPU_score$rich_sprawl > 2.3 ~ 0,
    TRUE ~ ((2.3-NAPU_score$rich_sprawl)/(2.3-1.0))*100)

# NCA Metric Scoring
NCA_score$score_aspt_mod <-
  case_when(
    NCA_score$aspt_mod > 7.54 ~ 100,
    NCA_score$aspt_mod < 6.67 ~ 0,
    TRUE ~ ((NCA_score$aspt_mod - 6.67)/(7.54-6.67))*100)
  
NCA_score$score_pct_cote <-
  case_when(
    NCA_score$pct_cote > 63.95 ~ 100,
    NCA_score$pct_cote < 47.71 ~ 0,
    TRUE ~ ((NCA_score$pct_cote - 47.71)/(63.95-47.71))*100)

NCA_score$score_pct_ephemeroptera <-
  case_when(
    NCA_score$pct_ephemeroptera > 40.73 ~ 100,
    NCA_score$pct_ephemeroptera < 9.02 ~ 0,
    TRUE ~ ((NCA_score$pct_ephemeroptera - 9.02)/(40.73-9.02))*100)

NCA_score$score_pct_heptageniidae <-
  case_when(
   NCA_score$pct_heptageniidae > 13.72 ~ 100,
   NCA_score$pct_heptageniidae  < 0 ~ 0,
   TRUE ~ ((NCA_score$pct_heptageniidae - 0)/(13.72-0))*100)

NCA_score$score_pct_limestone <-
  case_when(
    NCA_score$pct_limestone > 8.93 ~ 100,
    NCA_score$pct_limestone < 0 ~ 0,
    TRUE ~ ((NCA_score$pct_limestone - 0)/(8.93-0))*100)

NCA_score$score_rich_cling <-
  case_when(
    NCA_score$rich_cling > 10.0 ~ 100,
    NCA_score$rich_cling < 7.75 ~ 0,
    TRUE ~ ((NCA_score$rich_cling-7.75)/(10.0-7.75))*100)

NCA_score$score_rich_ephemeroptera <- 
  case_when(
    NCA_score$rich_ephemeroptera > 4.0 ~ 100,
    NCA_score$rich_ephemeroptera < 3.83 ~ 0,
    TRUE ~ ((NCA_score$rich_ephemeroptera - 3.83)/(4.0-3.83))*100)

NCA_score$score_rich_ept <-
  case_when(
    NCA_score$rich_ept > 12.0 ~ 100,
    NCA_score$rich_ept < 8.02 ~ 0,
    TRUE ~ ((NCA_score$rich_ept-8.02)/(12.0-8.02))*100)

NCA_score$score_rich_plecoptera <-
  case_when(
    NCA_score$rich_plecoptera > 4.04 ~ 100,
    NCA_score$rich_plecoptera < 3.45 ~ 0,
    TRUE ~ ((NCA_score$rich_plecoptera - 3.45)/(4.04-3.45))*100)

# CA Metric Scoring
CA_score$score_pct_ept_rich <-
  case_when(
    CA_score$pct_ept_rich > 66.82 ~ 100,
    CA_score$pct_ept_rich < 56.82 ~ 0,
    TRUE ~ ((CA_score$pct_ept_rich-56.82)/(66.82-56.82))*100)

CA_score$score_pct_euholognatha <-
  case_when(
    CA_score$pct_euholognatha > 17.98 ~ 100,
    CA_score$pct_euholognatha < 0.23 ~ 0,
    TRUE ~ ((CA_score$pct_euholognatha-0.23)/(17.98-0.23))*100)

CA_score$score_pct_heptageniidae <-
  case_when(
    CA_score$pct_heptageniidae > 6.69 ~ 100,
    CA_score$pct_heptageniidae < 0 ~ 0,
    TRUE ~ ((CA_score$pct_heptageniidae-0)/(6.69-0))*100)

CA_score$score_pct_hydro_ept <-
  case_when(
    CA_score$pct_hydro_ept < 12.5 ~ 100,
    CA_score$pct_hydro_ept > 21.73 ~ 0,
    TRUE ~ ((21.73-CA_score$pct_hydro_ept)/(21.73-12.5))*100)

CA_score$score_pct_mod_tol_4_6 <- 
  case_when(
    CA_score$pct_mod_tol_4_6 < 43.04 ~ 100,
    CA_score$pct_mod_tol_4_6 > 57.26 ~ 0,
    TRUE ~ ((57.26-CA_score$pct_mod_tol_4_6)/(57.26-43.04))*100)

CA_score$score_pct_non_hydrop_trichoptera <-
  case_when(
    CA_score$pct_non_hydrop_trichoptera > 40.4 ~ 100,
    CA_score$pct_non_hydrop_trichoptera < 25.02 ~ 0,
    TRUE ~ ((CA_score$pct_non_hydrop_trichoptera - 25.02)/(40.4-25.02))*100)

CA_score$score_pct_pisciforma <-
  case_when(
    CA_score$pct_pisciforma > 12.33 ~ 100,
    CA_score$pct_pisciforma < 5.49 ~ 0,
    TRUE ~ ((CA_score$pct_pisciforma-5.49)/(12.33-5.49))*100)

CA_score$score_pct_systellognatha <-
  case_when(
    CA_score$pct_systellognatha > 5.44 ~ 100,
    CA_score$pct_systellognatha < 0.08 ~ 0,
    TRUE ~ ((CA_score$pct_systellognatha-0.08)/(5.44-0.08))*100)

CA_score$score_rich_burrow <- 
  case_when(
    CA_score$rich_burrow > 3.0 ~ 100,
    CA_score$rich_burrow < 1.4 ~ 0,
    TRUE ~ ((CA_score$rich_burrow-1.4)/(3.0-1.4))*100)

CA_score$score_rich_collect <- 
  case_when(
    CA_score$rich_collect > 7.0 ~ 100,
    CA_score$rich_collect < 5.13 ~ 0,
    TRUE ~ ((CA_score$rich_collect-5.13)/(7.0-5.13))*100)

CA_score$score_rich_filter <-
  case_when(
    CA_score$rich_filter > 2.91 ~ 100,
    CA_score$rich_filter < 1.38 ~ 0,
    TRUE ~ ((CA_score$rich_filter-1.38)/(2.91-1.38))*100)

CA_score$score_rich_plecoptera <-
  case_when(
    CA_score$rich_plecoptera > 4.0 ~ 100,
    CA_score$rich_plecoptera < 0.47 ~ 0,
    TRUE ~ ((CA_score$rich_plecoptera-0.47)/(4.0-0.47))*100)

CA_score$score_rich_predator <- 
  case_when(
    CA_score$rich_predator > 3.0 ~ 100,
    CA_score$rich_predator < 1.44 ~ 0,
    TRUE ~ ((CA_score$rich_predator-1.44)/(3.0-1.44))*100)

# MAC Metric Scoring
MAC_score$score_aspt_mod <-
  case_when(
    MAC_score$aspt_mod > 3.27 ~ 100,
    MAC_score$aspt_mod < 2.9 ~ 0,
    TRUE ~ ((MAC_score$aspt_mod-2.9)/(3.27-2.9))*100)

MAC_score$score_gold <-
  case_when(
    MAC_score$gold > 58.96 ~ 100,
    MAC_score$gold < 33.91 ~ 0,
    TRUE ~ ((MAC_score$gold-33.91)/(58.96-33.91))*100)

MAC_score$score_pct_arthropoda <-
  case_when(
    MAC_score$pct_arthropoda < 74.01 ~ 100,
    MAC_score$pct_arthropoda > 95.38 ~ 0,
    TRUE ~ ((95.38-MAC_score$pct_arthropoda)/(95.38-74.01))*100)

MAC_score$score_pct_chironomidae <-
  case_when(
    MAC_score$pct_chironomidae < 23.01 ~ 100,
    MAC_score$pct_chironomidae > 42.68 ~ 0,
    TRUE ~ ((42.68-MAC_score$pct_chironomidae)/(42.68-23.01))*100)
  
MAC_score$score_pct_collect <-
  case_when(
    MAC_score$pct_collect < 73.25 ~ 100,
    MAC_score$pct_collect > 96.19 ~ 0,
    TRUE ~ ((96.19 - MAC_score$pct_collect)/(96.19-73.25))*100)

MAC_score$score_pct_mod_tol_4_6 <-
  case_when(
    MAC_score$pct_mod_tol_4_6 < 52.49 ~ 100,
    MAC_score$pct_mod_tol_4_6 > 75.48 ~ 0,
    TRUE ~ ((75.48-MAC_score$pct_mod_tol_4_6)/(75.48-52.49))*100)

MAC_score$score_pct_odonata <-
  case_when(
    MAC_score$pct_odonata > 1.82 ~ 100,
    MAC_score$pct_odonata < 0.43 ~ 0,
    TRUE ~ ((MAC_score$pct_odonata - 0.43)/(1.82-0.43))*100)

MAC_score$score_pct_predator <-
  case_when(
    MAC_score$pct_predator > 5.16 ~ 100,
    MAC_score$pct_predator < 4.01 ~ 0,
    TRUE ~ ((MAC_score$pct_predator - 4.01)/(5.16-4.01))*100)

MAC_score$score_pct_scrape <-
  case_when(
    MAC_score$pct_scrape > 12.07 ~ 100,
    MAC_score$pct_scrape < 0 ~ 0,
    TRUE ~ ((MAC_score$pct_scrape-0)/(12.07-0))*100)

MAC_score$score_pct_urban_intol <-
  case_when(
    MAC_score$pct_urban_intol < 65.28 ~ 100,
    MAC_score$pct_urban_intol > 71.78 ~ 0,
    TRUE ~ ((71.78-MAC_score$pct_urban_intol)/(71.78-65.28))*100)

MAC_score$score_rich_climb <-
  case_when(
    MAC_score$rich_climb > 2.45 ~ 100,
    MAC_score$rich_climb < 0 ~ 0,
    TRUE ~ ((MAC_score$rich_climb-0)/(2.45-0))*100)

MAC_score$score_rich_tol <-
  case_when(
    MAC_score$rich_tol > 4.0 ~ 100,
    MAC_score$rich_tol < 2.16 ~ 0,
    TRUE ~ ((MAC_score$rich_tol-2.16)/(4.0-2.16))*100)

# SEP Metric Scoring
SEP_score$score_aspt_mod <-
  case_when(
    SEP_score$aspt_mod > 4.65 ~ 100,
    SEP_score$aspt_mod < 3.18 ~ 0,
    TRUE ~ ((SEP_score$aspt_mod-3.18)/(4.65-3.18))*100)

SEP_score$score_gold <-
  case_when(
    SEP_score$gold > 50.3 ~ 100,
    SEP_score$gold < 32.16 ~ 0,
    TRUE ~ ((SEP_score$gold-32.16)/(50.3-32.16))*100)

SEP_score$score_hbi <-
  case_when(
    SEP_score$hbi < 5.12 ~ 100,
    SEP_score$hbi > 5.86 ~ 0,
    TRUE ~ ((5.86-SEP_score$hbi)/(5.86-5.12))*100)

SEP_score$pct_cling <-
  case_when(
    SEP_score$pct_cling > 43.05 ~ 100,
    SEP_score$pct_cling < 0 ~ 0,
    TRUE ~ ((SEP_score$pct_cling-0)/(43.05-0))*100)

SEP_score$score_pct_cote <-
  case_when(
    SEP_score$pct_cote > 26.37 ~ 100,
    SEP_score$pct_cote < 1.62 ~ 0,
    TRUE ~ ((SEP_score$pct_cote-1.62)/(26.37-1.62))*100)

SEP_score$score_pct_ept_rich_no_tol <-
  case_when(
    SEP_score$pct_ept_rich_no_tol > 24.73 ~ 100,
    SEP_score$pct_ept_rich_no_tol < 16.49 ~ 0,
    TRUE ~ ((SEP_score$pct_ept_rich_no_tol-16.49)/(24.73-16.49))*100)

SEP_score$score_rich_cling <-
  case_when(
    SEP_score$rich_cling > 4.93 ~ 100,
    SEP_score$rich_cling < 1.03 ~ 0,
    TRUE ~ ((SEP_score$rich_cling-1.03)/(4.93-1.03))*100)

SEP_score$score_rich_collect <-
  case_when(
    SEP_score$rich_collect > 8.0 ~ 100,
    SEP_score$rich_collect < 4.62 ~ 0,
    TRUE ~ ((SEP_score$rich_collect-4.62)/(8.0-4.62))*100)

SEP_score$score_rich_ephemeroptera <-
  case_when(
    SEP_score$rich_ephemeroptera > 2.0 ~ 100,
    SEP_score$rich_ephemeroptera < 0 ~ 0,
    TRUE ~ ((SEP_score$rich_ephemeroptera-0)/(2.0-0))*100)

SEP_score$score_rich_filter <-
  case_when(
    SEP_score$rich_filter  > 2.31 ~ 100,
    SEP_score$rich_filter < 0.04 ~ 0,
    TRUE ~ ((SEP_score$rich_filter-0.04)/(2.31-0.04))*100)

SEP_score$score_rich_trichoptera <-
  case_when(
    SEP_score$rich_trichoptera > 2.0 ~ 100,
    SEP_score$rich_trichoptera < 0.41 ~ 0,
    TRUE ~ ((SEP_score$rich_trichoptera-0.41)/(2.0-0.41))*100)

# BLUE Metric Scoring
BLUE_score$score_pct_burrow <-
  case_when(
    BLUE_score$pct_burrow < 12.65 ~ 100,
    BLUE_score$pct_burrow > 40.71 ~ 0,
    TRUE ~ ((40.71-BLUE_score$pct_burrow)/(40.71-12.65))*100)

BLUE_score$score_pct_ept_rich_no_tol <-
  case_when(
    BLUE_score$pct_ept_rich_no_tol > 52.78 ~ 100,
    BLUE_score$pct_ept_rich_no_tol < 34.74 ~ 0,
    TRUE ~ ((BLUE_score$pct_ept_rich_no_tol-34.74)/(52.78-34.74))*100)

BLUE_score$score_pct_scrape <-
  case_when(
    BLUE_score$pct_scrape > 15.23 ~ 100,
    BLUE_score$pct_scrape < 0 ~ 0,
    TRUE ~ ((BLUE_score$pct_scrape-0)/(15.23-0))*100)

BLUE_score$score_pct_systellognatha <-
  case_when(
    BLUE_score$pct_systellognatha > 6.06 ~ 100,
    BLUE_score$pct_systellognatha < 0 ~ 0,
    TRUE ~ ((BLUE_score$pct_systellognatha-0)/(6.06-0))*100)

BLUE_score$score_rich_ephemeroptera <-
  case_when(
    BLUE_score$rich_ephemeroptera > 4.0 ~ 100,
    BLUE_score$rich_ephemeroptera < 0.33 ~ 0,
    TRUE ~ ((BLUE_score$rich_ephemeroptera-0.33)/(4.0-0.33))*100)

# NRV Metric Scoring
NRV_score$score_aspt_mod <-
  case_when(
    NRV_score$aspt_mod > 6.27 ~ 100,
    NRV_score$aspt_mod < 4.59 ~ 0,
    TRUE ~ ((NRV_score$aspt_mod-4.59)/(6.27-4.59))*100)

NRV_score$score_pct_ephemeroptera_no_baetid <-
  case_when(
    NRV_score$pct_ephemeroptera_no_baetid > 23.69 ~ 100,
    NRV_score$pct_ephemeroptera_no_baetid < 6.01 ~ 0,
    TRUE ~ ((NRV_score$pct_ephemeroptera_no_baetid-6.01)/(23.69-6.01))*100)

NRV_score$score_pct_ept_rich <-
  case_when(
    NRV_score$pct_ept_rich > 64.7 ~ 100,
    NRV_score$pct_ept_rich < 47.53 ~ 0,
    TRUE ~ ((NRV_score$pct_ept_rich-47.53)/(64.7-47.53))*100)

NRV_score$score_rich_collect <-
  case_when(
    NRV_score$rich_collect > 8.0 ~ 100,
    NRV_score$rich_collect < 7.14 ~ 0,
    TRUE ~ ((NRV_score$rich_collect-7.14)/(8.0-7.14))*100)

NRV_score$score_rich_ephemeroptera <-
  case_when(
    NRV_score$rich_ephemeroptera > 4.0 ~ 100,
    NRV_score$rich_ephemeroptera < 1.8 ~ 0,
    TRUE ~ ((NRV_score$rich_ephemeroptera-1.8)/(4.0-1.8))*100)

NRV_score$score_rich_predator <-
  case_when(
    NRV_score$rich_predator > 3.0 ~ 100,
    NRV_score$rich_predator < 0.32 ~ 0,
    TRUE ~ ((NRV_score$rich_predator-0.32)/(3.0-0.32))*100)

NRV_score$score_rich_shred <-
  case_when(
    NRV_score$rich_shred > 2.0 ~ 100,
    NRV_score$rich_shred < 1.27 ~ 0,
    TRUE ~ ((NRV_score$rich_shred-1.27)/(2.0-1.27))*100)

NRV_score$score_rich_trichoptera <-
  case_when(
    NRV_score$rich_trichoptera > 3.0 ~ 100,
    NRV_score$rich_trichoptera < 2.53 ~ 0,
    TRUE ~ ((NRV_score$rich_trichoptera-2.53)/(3.0-2.53))*100)

# UNP Metric Scoring
UNP_score$score_hbi <-
  case_when(
    UNP_score$hbi < 4.01 ~ 100,
    UNP_score$hbi > 5.74 ~ 0,
    TRUE ~ ((5.74-UNP_score$hbi)/(5.74-4.01))*100)

UNP_score$score_pct_collect <-
  case_when(
    UNP_score$pct_collect < 80.32 ~ 100,
    UNP_score$pct_collect > 90.97 ~ 0,
    TRUE ~ ((90.97-UNP_score$pct_collect)/(90.97-80.32))*100)

UNP_score$score_pct_ept_rich_no_tol <-
  case_when(
    UNP_score$pct_ept_rich_no_tol > 37.09 ~ 100,
    UNP_score$pct_ept_rich_no_tol < 22.1 ~ 0,
    TRUE ~ ((UNP_score$pct_ept_rich_no_tol-22.1)/(37.09-22.1))*100)

UNP_score$score_pct_filter <-
  case_when(
    UNP_score$pct_filter > 18.14 ~ 100,
    UNP_score$pct_filter < 1.91 ~ 0,
    TRUE ~ ((UNP_score$pct_filter - 1.91)/(18.14-1.91))*100)

UNP_score$score_pct_hexapoda <-
  case_when(
    UNP_score$pct_hexapoda > 99.5 ~ 100,
    UNP_score$pct_hexapoda < 96.59 ~ 0,
    TRUE ~ ((UNP_score$pct_hexapoda-96.59)/(99.5-96.59))*100)

UNP_score$score_pct_pisciforma <-
  case_when(
    UNP_score$pct_pisciforma > 9.43 ~ 100,
    UNP_score$pct_pisciforma < 0 ~ 0,
    TRUE ~ ((UNP_score$pct_pisciforma-0)/(9.43-0))*100)

UNP_score$score_pct_retreat_caddisfly <-
  case_when(
    UNP_score$pct_retreat_caddisfly > 10.97 ~ 100,
    UNP_score$pct_retreat_caddisfly < 2.89 ~ 0,
    TRUE ~ ((UNP_score$pct_retreat_caddisfly-2.89)/(10.97-2.89))*100)

UNP_score$score_pct_scrape <-
  case_when(
    UNP_score$pct_scrape > 5.72 ~ 100,
    UNP_score$pct_scrape < 0.22 ~ 0,
    TRUE ~ ((UNP_score$pct_scrape-0.22)/(5.72-0.22))*100)

UNP_score$score_pct_urban_intol <- 
  case_when(
    UNP_score$pct_urban_intol > 98.46 ~ 100,
    UNP_score$pct_urban_intol < 91.71 ~ 0,
    TRUE ~ ((UNP_score$pct_urban_intol-91.71)/(98.46-91.71))*100)

UNP_score$score_rich_ephemeroptera <-
  case_when(
    UNP_score$rich_ephemeroptera > 3.0 ~ 100,
    UNP_score$rich_ephemeroptera < 1.44 ~ 0,
    TRUE ~ ((UNP_score$rich_ephemeroptera-1.44)/(3.0-1.44))*100)

UNP_score$score_rich_plecoptera <-
  case_when(
    UNP_score$rich_plecoptera > 2.23 ~ 100,
    UNP_score$rich_plecoptera < 0 ~ 0,
    TRUE ~ ((UNP_score$rich_plecoptera-0)/(2.23-0))*100)

UNP_score$score_rich_trichoptera <-
  case_when(
    UNP_score$rich_trichoptera > 3.0 ~ 100,
    UNP_score$rich_trichoptera < 1.01 ~ 0,
    TRUE ~ ((UNP_score$rich_trichoptera-1.01)/(3.0-1.01))*100)

# SRV Metric Scoring
SRV_score$score_aspt_mod <-
  case_when(
    SRV_score$aspt_mod > 6.44 ~ 100,
    SRV_score$aspt_mod < 4.95 ~ 0,
    TRUE ~ ((SRV_score$aspt_mod-4.95)/(6.44-4.95))*100)

SRV_score$score_pct_burrow <-
  case_when(
    SRV_score$pct_burrow < 17.13 ~ 100,
    SRV_score$pct_burrow > 31.28 ~ 0,
    TRUE ~ ((31.28-SRV_score$pct_burrow)/(31.28-17.13))*100)

SRV_score$score_pct_cling <-
  case_when(
    SRV_score$pct_cling > 58.45 ~ 100,
    SRV_score$pct_cling < 40.16 ~ 100,
    TRUE ~ ((SRV_score$pct_cling-40.16)/(58.45-40.16))*100)

SRV_score$score_pct_cote <-
  case_when(
    SRV_score$pct_cote > 57.55 ~ 100,
    SRV_score$pct_cote < 43.78 ~ 0,
    TRUE ~ ((SRV_score$pct_cote-43.78)/(57.55-43.78))*100)

SRV_score$score_pct_ephemeroptera_no_baetid <-
  case_when(
    SRV_score$pct_ephemeroptera_no_baetid > 21.0 ~ 100,
    SRV_score$pct_ephemeroptera_no_baetid < 2.8 ~ 0,
    TRUE ~ ((SRV_score$pct_ephemeroptera_no_baetid-2.8)/(21.0-2.8))*100)

SRV_score$score_pct_ept_rich_no_tol <-
  case_when(
    SRV_score$pct_ept_rich_no_tol > 47.16 ~ 100,
    SRV_score$pct_ept_rich_no_tol < 34.94 ~ 0,
    TRUE ~ ((SRV_score$pct_ept_rich_no_tol-34.94)/(47.16-34.94))*100)

SRV_score$score_pct_predator <-
  case_when(
    SRV_score$pct_predator > 7.12 ~ 100,
    SRV_score$pct_predator < 2.24 ~ 0,
    TRUE ~ ((SRV_score$pct_predator-2.24)/(7.12-2.24))*100)

SRV_score$score_pct_scrape <-
  case_when(
    SRV_score$pct_scrape > 10.88 ~ 100,
    SRV_score$pct_scrape < 0.47 ~ 0,
    TRUE ~ ((SRV_score$pct_scrape-0.47)/(10.88-0.47))*100)

SRV_score$score_rich_cling <- 
  case_when(
    SRV_score$rich_cling > 8.43 ~ 100,
    SRV_score$rich_cling < 6.16 ~ 0,
    TRUE ~ ((SRV_score$rich_cling-6.16)/(8.43-6.16))*100)

SRV_score$score_rich_collect <-
  case_when(
    SRV_score$rich_collect > 7.97 ~ 100,
    SRV_score$rich_collect < 6.02 ~ 0,
    TRUE ~ ((SRV_score$rich_collect-6.02)/(7.97-6.02))*100)

SRV_score$score_rich_ephemeroptera <-
  case_when(
    SRV_score$rich_ephemeroptera > 4.0 ~ 100,
    SRV_score$rich_ephemeroptera < 3.66 ~ 0,
    TRUE ~ ((SRV_score$rich_ephemeroptera-3.66)/(4.0-3.66))*100)

SRV_score$score_rich_intol <-
  case_when(
    SRV_score$rich_intol > 8.0 ~ 100,
    SRV_score$rich_intol < 4.15 ~ 0,
    TRUE ~ ((SRV_score$rich_intol-4.15)/(8.0-4.15))*100)

SRV_score$score_rich_plecoptera <-
  case_when(
    SRV_score$rich_plecoptera > 3.0 ~ 100,
    SRV_score$rich_plecoptera < 0.65 ~ 0,
    TRUE ~ ((SRV_score$rich_plecoptera-0.65)/(3.0-0.65))*100)

SRV_score$score_rich_predator <-
  case_when(
    SRV_score$rich_predator > 3.0 ~ 100,
    SRV_score$rich_predator < 2.72 ~ 0,
    TRUE ~ ((SRV_score$rich_predator-2.72)/(3.0-2.72))*100)

SRV_score$score_rich_trichoptera <-
  case_when(
    SRV_score$rich_trichoptera > 3.0 ~ 100,
    SRV_score$rich_trichoptera < 2.79 ~ 0,
    TRUE ~ ((SRV_score$rich_trichoptera-2.79)/(3.0-2.79))*100)

# SGV Metric Scoring
SGV_score$score_aspt_mod <-
  case_when(
    SGV_score$aspt_mod > 6.44 ~ 100,
    SGV_score$aspt_mod < 3.94 ~ 0,
    TRUE ~ ((SGV_score$aspt_mod-3.94)/(6.44-3.94))*100)

SGV_score$score_gold <-
  case_when(
    SGV_score$gold > 80.89 ~ 100,
    SGV_score$gold < 64.95 ~ 0,
    TRUE ~ ((SGV_score$gold-64.95)/(80.89-64.95))*100)

SGV_score$score_pct_cling <-
  case_when(
    SGV_score$pct_cling > 62.23 ~ 100,
    SGV_score$pct_cling < 39.7 ~ 0,
    TRUE ~ ((SGV_score$pct_cling-39.7)/(62.23-39.7))*100)

SGV_score$score_pct_collect <-
  case_when(
    SGV_score$pct_collect < 71.17 ~ 100,
    SGV_score$pct_collect > 94.57 ~ 0,
    TRUE ~ ((94.57-SGV_score$pct_collect)/(94.57-71.17))*100)

SGV_score$score_pct_ephemeroptera_no_baetid <-
  case_when(
    SGV_score$pct_ephemeroptera_no_baetid > 23.94 ~ 100,
    SGV_score$pct_ephemeroptera_no_baetid < 5.23 ~ 0,
    TRUE ~ ((SGV_score$pct_ephemeroptera_no_baetid-5.23)/(23.94-5.23))*100)

SGV_score$score_pct_ept_no_hydro <-
  case_when(
    SGV_score$pct_ept_no_hydro > 83.21 ~ 100,
    SGV_score$pct_ept_no_hydro < 73.53 ~ 0,
    TRUE ~ ((SGV_score$pct_ept_no_hydro-73.53)/(83.21-73.53))*100)

SGV_score$score_pct_ept_rich_no_tol <-
  case_when(
    SGV_score$pct_ept_rich_no_tol > 32.97 ~ 100,
    SGV_score$pct_ept_rich_no_tol < 19.39 ~ 0,
    TRUE ~ ((SGV_score$pct_ept_rich_no_tol-19.39)/(32.97-19.39))*100)

SGV_score$score_pct_filter <-
  case_when(
    SGV_score$pct_filter > 22.79 ~ 100,
    SGV_score$pct_filter < 8.6 ~ 0,
    TRUE ~ ((SGV_score$pct_filter-8.6)/(22.79-8.6))*100)

SGV_score$score_pct_heptageniidae <-
  case_when(
    SGV_score$pct_heptageniidae > 7.89 ~ 100,
    SGV_score$pct_heptageniidae < 0 ~ 0,
    TRUE ~ ((SGV_score$pct_heptageniidae-0)/(7.89-0))*100)

SGV_score$score_pct_mod_tol_4_6 <-
  case_when(
    SGV_score$pct_mod_tol_4_6 < 52.49 ~ 100,
    SGV_score$pct_mod_tol_4_6 > 76.71 ~ 0,
    TRUE ~ ((76.71-SGV_score$pct_mod_tol_4_6)/(76.71-52.49))*100)

SGV_score$score_pct_non_hydrop_trichoptera <-
  case_when(
    SGV_score$pct_non_hydrop_trichoptera > 30.54 ~ 100,
    SGV_score$pct_non_hydrop_trichoptera < 10.57 ~ 0,
    TRUE ~ ((SGV_score$pct_non_hydrop_trichoptera-10.57)/(30.54-10.57))*100)

SGV_score$score_pct_pisciforma <-
  case_when(
    SGV_score$pct_pisciforma > 17.07 ~ 100,
    SGV_score$pct_pisciforma < 0.1 ~ 0,
    TRUE ~ ((SGV_score$pct_pisciforma-0.1)/(17.07-0.1))*100)

SGV_score$score_pct_predator <-
  case_when(
    SGV_score$pct_predator > 5.28 ~ 100,
    SGV_score$pct_predator < 0.02 ~ 0,
    TRUE ~ ((SGV_score$pct_predator-0.02)/(5.28-0.02))*100)

SGV_score$score_pct_pterygota <-
  case_when(
    SGV_score$pct_pterygota > 96.88 ~ 100,
    SGV_score$pct_pterygota < 90.2 ~ 0,
    TRUE ~ ((SGV_score$pct_pterygota-90.2)/(96.88-90.2))*100)

SGV_score$score_pct_scrape <-
  case_when(
    SGV_score$pct_scrape > 15.43 ~ 100,
    SGV_score$pct_scrape < 0.32 ~ 0,
    TRUE ~ ((SGV_score$pct_scrape-0.32)/(15.43-0.32))*100)

SGV_score$score_pct_sprawl <-
  case_when(
    SGV_score$pct_sprawl < 1.79 ~ 100,
    SGV_score$pct_sprawl > 4.54 ~ 0,
    TRUE ~ ((4.54-SGV_score$pct_sprawl)/(4.54-1.79))*100)

SGV_score$score_pct_tolerant_7_10 <-
  case_when(
    SGV_score$pct_tolerant_7_10 < 1.73 ~ 100,
    SGV_score$pct_tolerant_7_10 > 4.3 ~ 0,
    TRUE ~ ((4.3-SGV_score$pct_tolerant_7_10)/(4.3-1.73))*100)

SGV_score$score_rich_collect <-
  case_when(
    SGV_score$rich_collect > 8.0 ~ 100,
    SGV_score$rich_collect < 6.32 ~ 0,
    TRUE ~ ((SGV_score$rich_collect-6.32)/(8.0-6.32))*100)

SGV_score$score_rich_ephemeroptera <-
  case_when(
    SGV_score$rich_ephemeroptera > 3.94 ~ 100,
    SGV_score$rich_ephemeroptera < 0.61 ~ 0,
    TRUE ~ ((SGV_score$rich_ephemeroptera-0.61)/(3.94-0.61))*100)

SGV_score$score_rich_filter <-
  case_when(
    SGV_score$rich_filter > 3.0 ~ 100,
    SGV_score$rich_filter < 1.36 ~ 0,
    TRUE ~ ((SGV_score$rich_filter-1.36)/(3.0-1.36))*100)

SGV_score$score_rich_trichoptera <-
  case_when(
    SGV_score$rich_trichoptera > 2.54 ~ 100,
    SGV_score$rich_trichoptera < 0.08 ~ 0,
    TRUE ~ ((SGV_score$rich_trichoptera-0.08)/(2.54-0.08))*100)

# LNP Metric Scoring
LNP_score$score_aspt_mod <-
  case_when(
    LNP_score$aspt_mod > 6.62 ~ 100,
    LNP_score$aspt_mod < 2.7 ~ 0,
    TRUE ~ ((LNP_score$aspt_mod-2.7)/(6.62-2.7))*100)

LNP_score$score_pct_cling <-
  case_when(
    LNP_score$pct_cling > 62.77 ~ 100,
    LNP_score$pct_cling < 4.54 ~ 0,
    TRUE ~ ((LNP_score$pct_cling-4.54)/(62.77-4.54))*100)

LNP_score$score_pct_collect <-
  case_when(
    LNP_score$pct_collect < 72.82 ~ 100,
    LNP_score$pct_collect > 92.49 ~ 0,
    TRUE ~ ((92.49-LNP_score$pct_collect)/(92.49-72.82))*100)

LNP_score$score_pct_ept_rich <-
  case_when(
    LNP_score$pct_ept_rich > 52.32 ~ 100,
    LNP_score$pct_ept_rich < 38.07 ~ 0,
    TRUE ~ ((LNP_score$pct_ept_rich-38.07)/(52.32-38.07))*100)

LNP_score$score_pct_mod_tol_4_6 <-
  case_when(
    LNP_score$pct_mod_tol_4_6 < 46.54 ~ 100,
    LNP_score$pct_mod_tol_4_6 > 77.78 ~ 0,
    TRUE ~ ((77.78-LNP_score$pct_mod_tol_4_6)/(77.78-46.54))*100)

LNP_score$score_pct_predator <-
  case_when(
    LNP_score$pct_predator > 7.78 ~ 100,
    LNP_score$pct_predator < 1.43 ~ 0,
    TRUE ~ ((LNP_score$pct_predator-1.43)/(7.78-1.43))*100)

LNP_score$score_pct_pterygota <-
  case_when(
    LNP_score$pct_pterygota > 98.36 ~ 100,
    LNP_score$pct_pterygota < 92.09 ~ 0,
    TRUE ~ ((LNP_score$pct_pterygota-92.09)/(98.36-92.09))*100)

LNP_score$score_rich_cling <-
  case_when(
    LNP_score$rich_cling > 7.0 ~ 100,
    LNP_score$rich_cling < 3.69 ~ 0,
    TRUE ~ ((LNP_score$rich_cling-3.69)/(7.0-3.69))*100)

LNP_score$score_rich_ephemeroptera <-
  case_when(
    LNP_score$rich_ephemeroptera > 3.0 ~ 100,
    LNP_score$rich_ephemeroptera < 1.7 ~ 0,
    TRUE ~ ((LNP_score$rich_ephemeroptera-1.7)/(3.0-1.7))*100)

# PIED Metric Scoring
PIED_score$score_pct_cling <-
  case_when(
    PIED_score$pct_cling > 62.62 ~ 100,
    PIED_score$pct_cling < 32.33 ~ 0,
    TRUE ~ ((PIED_score$pct_cling-32.33)/(62.62-32.33))*100)

PIED_score$score_pct_ept_rich <-
  case_when(
    PIED_score$pct_ept_rich > 53.58 ~ 100,
    PIED_score$pct_ept_rich < 38.51 ~ 0,
    TRUE ~ ((PIED_score$pct_ept_rich-38.51)/(53.58-38.51))*100)

PIED_score$score_pct_intol_0_4 <-
  case_when(
    PIED_score$pct_intol_0_4 > 58.15 ~ 100,
    PIED_score$pct_intol_0_4 < 27.97 ~ 0,
    TRUE ~ ((PIED_score$pct_intol_0_4-27.97)/(58.15-27.97))*100)

PIED_score$score_pct_tolerant_7_10 <- 
  case_when(
    PIED_score$pct_tolerant_7_10 < 1.6 ~ 100,
    PIED_score$pct_tolerant_7_10 > 4.84 ~ 0,
    TRUE ~ ((4.84-PIED_score$pct_tolerant_7_10)/(4.84-1.6))*100)

PIED_score$score_pct_urban_intol <-
  case_when(
    PIED_score$pct_urban_intol > 95.41 ~ 100,
    PIED_score$pct_urban_intol < 86.93 ~ 0,
    TRUE ~ ((PIED_score$pct_urban_intol-86.93)/(95.41-86.93))*100)

PIED_score$score_rich_cling <-
  case_when(
    PIED_score$rich_cling > 8.0 ~ 100,
    PIED_score$rich_cling < 2.96 ~ 0,
    TRUE ~ ((PIED_score$rich_cling-2.96)/(8.0-2.96))*100)
  
PIED_score$score_rich_ephemeroptera <-  
  case_when(
    PIED_score$rich_ephemeroptera > 3.6 ~ 100,
    PIED_score$rich_ephemeroptera < 0.86 ~ 0,
    TRUE ~ ((PIED_score$rich_ephemeroptera-0.86)/(3.6-0.86))*100)

PIED_score$score_rich_trichoptera <-
  case_when(
    PIED_score$rich_trichoptera > 2.03 ~ 100,
    PIED_score$rich_trichoptera < 0.45 ~ 0,
    TRUE ~ ((PIED_score$rich_trichoptera-0.45)/(2.03-0.45))*100)
```
**4. Calculating the Chessie BIBI Scores**  
Average the metric score columns for each unique_id. **If you added or removed columns, you will need to adjust the column ranges!**

The resulting mean value is the Chessie BIBI for each unique_id, and it is located in a new column called *chessie_bibi*.
```{r}
NAPU_score$chessie_bibi <- rowMeans(NAPU_score[ ,25:43], na.rm = T)
NCA_score$chessie_bibi  <- rowMeans(NCA_score[ ,15:23], na.rm = T)
CA_score$chessie_bibi   <- rowMeans(CA_score[ ,19:31], na.rm = T)
MAC_score$chessie_bibi  <- rowMeans(MAC_score[ ,18:29], na.rm = T)
SEP_score$chessie_bibi  <- rowMeans(SEP_score[ ,17:26], na.rm = T)
BLUE_score$chessie_bibi <- rowMeans(BLUE_score[ ,11:15], na.rm = T) 
NRV_score$chessie_bibi  <- rowMeans(NRV_score[ ,14:21], na.rm = T)
SRV_score$chessie_bibi  <- rowMeans(SRV_score[ ,21:35], na.rm = T)
UNP_score$chessie_bibi  <- rowMeans(UNP_score[ ,18:29], na.rm = T)
SGV_score$chessie_bibi  <- rowMeans(SGV_score[ ,27:47], na.rm = T) 
LNP_score$chessie_bibi  <- rowMeans(LNP_score[ ,15:23], na.rm = T)
PIED_score$chessie_bibi <- rowMeans(PIED_score[ ,14:21], na.rm = T) 
```
###**5. Assigning Narrative Ratings** 
Narrative ratings are assigned as Excellent, Good, Fair, Poor, or Very Poor based on the Chessie BIBI scores. The assignments are based on Table 5, Page 27 in Smith et al. (2017).
```{r}
NAPU_score$rating <-
  case_when(
    NAPU_score$chessie_bibi < 17.93 ~ "very_poor",
    NAPU_score$chessie_bibi >= 17.93 & NAPU_score$chessie_bibi < 35.85 ~ "poor",
    NAPU_score$chessie_bibi >= 35.85 & NAPU_score$chessie_bibi < 47.42 ~ "fair",
    NAPU_score$chessie_bibi >= 47.42 & NAPU_score$chessie_bibi < 61.95 ~ "good",
    NAPU_score$chessie_bibi >= 61.95 ~ "excellent",
    TRUE ~ "ERROR")

NCA_score$rating <-
  case_when(
    NCA_score$chessie_bibi < 15.73 ~ "very_poor",
    NCA_score$chessie_bibi >= 15.73 & NCA_score$chessie_bibi < 31.46 ~ "poor",
    NCA_score$chessie_bibi >= 31.46 & NCA_score$chessie_bibi < 56.25 ~ "fair",
    NCA_score$chessie_bibi >= 56.25 & NCA_score$chessie_bibi < 78.66 ~ "good",
    NCA_score$chessie_bibi >= 78.66 ~ "excellent",
    TRUE ~ "ERROR")
  
CA_score$rating <-
  case_when(
    CA_score$chessie_bibi < 27.88 ~ "very_poor",
    CA_score$chessie_bibi >= 27.88 & CA_score$chessie_bibi < 55.77 ~ "poor",
    CA_score$chessie_bibi >= 55.77 & CA_score$chessie_bibi < 69.55 ~ "fair",
    CA_score$chessie_bibi >= 69.55 & CA_score$chessie_bibi < 78.56 ~ "good",
    CA_score$chessie_bibi >= 78.56 ~ "excellent",
    TRUE ~ "ERROR")

MAC_score$rating <-
  case_when(
    MAC_score$chessie_bibi < 22.72 ~ "very_poor",
    MAC_score$chessie_bibi >= 22.72 & MAC_score$chessie_bibi < 45.43 ~ "poor",
    MAC_score$chessie_bibi >= 45.43 & MAC_score$chessie_bibi < 63.18 ~ "fair",
    MAC_score$chessie_bibi >= 63.18 & MAC_score$chessie_bibi < 76.85 ~ "good",
    MAC_score$chessie_bibi >= 76.85 ~ "excellent",
    TRUE ~ "ERROR")
  
SEP_score$rating <-
  case_when(
    SEP_score$chessie_bibi < 16.28 ~ "very_poor",
    SEP_score$chessie_bibi >= 16.28 & SEP_score$chessie_bibi < 32.55 ~ "poor",
    SEP_score$chessie_bibi >= 32.55 & SEP_score$chessie_bibi < 56.66 ~ "fair",
    SEP_score$chessie_bibi >= 56.66 & SEP_score$chessie_bibi < 83.90 ~ "good",
    SEP_score$chessie_bibi >= 83.90 ~ "excellent",
    TRUE ~ "ERROR")
  
BLUE_score$rating <-
  case_when(
    BLUE_score$chessie_bibi < 30.87 ~ "very_poor",
    BLUE_score$chessie_bibi >= 30.87 & BLUE_score$chessie_bibi < 61.74 ~ "poor",
    BLUE_score$chessie_bibi >= 61.74 & BLUE_score$chessie_bibi < 82.21 ~ "fair",
    BLUE_score$chessie_bibi >= 82.21 & BLUE_score$chessie_bibi < 92.07 ~ "good",
    BLUE_score$chessie_bibi >= 92.07 ~ "excellent",
    TRUE ~ "ERROR")
  
NRV_score$rating <-
  case_when(
    NRV_score$chessie_bibi < 19.10 ~ "very_poor",
    NRV_score$chessie_bibi >= 19.10 & NRV_score$chessie_bibi < 38.21 ~ "poor",
    NRV_score$chessie_bibi >= 38.21 & NRV_score$chessie_bibi < 50.81 ~ "fair",
    NRV_score$chessie_bibi >= 50.81 & NRV_score$chessie_bibi < 75.00 ~ "good",
    NRV_score$chessie_bibi >= 75.00 ~ "excellent",
    TRUE ~ "ERROR")
  
SRV_score$rating <-
  case_when(
    SRV_score$chessie_bibi < 21.58 ~ "very_poor",
    SRV_score$chessie_bibi >= 21.58 & SRV_score$chessie_bibi < 43.16 ~ "poor",
    SRV_score$chessie_bibi >= 43.16 & SRV_score$chessie_bibi < 57.97 ~ "fair",
    SRV_score$chessie_bibi >= 57.97 & SRV_score$chessie_bibi < 71.93 ~ "good",
    SRV_score$chessie_bibi >= 71.93 ~ "excellent",
    TRUE ~ "ERROR")

UNP_score$rating <-
  case_when(
    UNP_score$chessie_bibi < 31.30 ~ "very_poor",
    UNP_score$chessie_bibi >= 31.30 & UNP_score$chessie_bibi < 62.59 ~ "poor",
    UNP_score$chessie_bibi >= 62.59 & UNP_score$chessie_bibi < 69.85 ~ "fair",
    UNP_score$chessie_bibi >= 69.85 & UNP_score$chessie_bibi < 80.52 ~ "good",
    UNP_score$chessie_bibi >= 80.52 ~ "excellent",
    TRUE ~ "ERROR")

SGV_score$rating <-
  case_when(
    SGV_score$chessie_bibi < 29.93 ~ "very_poor",
    SGV_score$chessie_bibi >= 29.93 & SGV_score$chessie_bibi < 59.86 ~ "poor",
    SGV_score$chessie_bibi >= 59.86 & SGV_score$chessie_bibi < 66.02 ~ "fair",
    SGV_score$chessie_bibi >= 66.02 & SGV_score$chessie_bibi < 76.48 ~ "good",
    SGV_score$chessie_bibi >= 76.48 ~ "excellent",
    TRUE ~ "ERROR")

LNP_score$rating <-
  case_when(
    LNP_score$chessie_bibi < 35.16 ~ "very_poor",
    LNP_score$chessie_bibi >= 35.16 & LNP_score$chessie_bibi < 70.32 ~ "poor",
    LNP_score$chessie_bibi >= 70.32 & LNP_score$chessie_bibi < 80.15 ~ "fair",
    LNP_score$chessie_bibi >= 80.15 & LNP_score$chessie_bibi < 91.30 ~ "good",
    LNP_score$chessie_bibi >= 91.30 ~ "excellent",
    TRUE ~ "ERROR")

PIED_score$rating <-
  case_when(
    PIED_score$chessie_bibi < 30.31 ~ "very_poor",
    PIED_score$chessie_bibi >= 30.31 & PIED_score$chessie_bibi < 60.62 ~ "poor",
    PIED_score$chessie_bibi >= 60.62 & PIED_score$chessie_bibi < 73.18 ~ "fair",
    PIED_score$chessie_bibi >= 73.18 & PIED_score$chessie_bibi < 81.83 ~ "good",
    PIED_score$chessie_bibi >= 81.83 ~ "excellent",
    TRUE ~ "ERROR")
```
**6. Final Dataframe**  
Select the indicated columns from each bioregion's dataframe. **Remember to delete "sample_time" from the scripts if you did not include this column.**Merge these data into one large data frame, **bibi_ratings_combined**, using a series of `rbind()` calls. `rbind()` only accepts two data frames, which is why it is repeated several times. 
```{r}
NAPU_sub <- NAPU_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
NCA_sub <- NCA_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
CA_sub <- CA_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
MAC_sub <- MAC_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
SEP_sub <- SEP_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
BLUE_sub <- BLUE_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
NRV_sub <- NRV_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
SRV_sub <- SRV_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
UNP_sub <- UNP_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
SGV_sub <- SGV_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
LNP_sub <- LNP_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
PIED_sub <- PIED_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)

bind1 <- rbind(NAPU_sub, NCA_sub)
bind2 <- rbind(bind1, CA_sub)
bind3 <- rbind(bind2, MAC_sub)
bind4 <- rbind(bind3, SEP_sub)
bind5 <- rbind(bind4, BLUE_sub)
bind6 <- rbind(bind5, NRV_sub)
bind7 <- rbind(bind6, SRV_sub)
bind8 <- rbind(bind7, UNP_sub)
bind9 <- rbind(bind8, SGV_sub)
bind10 <- rbind(bind9, LNP_sub)
bibi_ratings_combined <- rbind(bind10, PIED_sub)
```
To export this final dataframe as a CSV, use the following script. Don't forget to replace the text with your own.
```{r}
setwd("WORKING DIRECTORY")
write.csv(bibi_ratings_combined, file= "FILE NAME.csv")
```
**Optional Next Step**  
You now have a dataframe that contains the Chessie BIBI score, narrative rating, and information about the sampling event for each unique_id in your dataset. It is advised that you proceed to the next step, which is verifying those scores and ratings via repetitions of the Chessie BIBI program.

The reason that we recommend completing the verification step is that the dataset and, ultimately, the scores produced by rarefaction vary slightly with each repetition. Taxa richness and diversity metrics for a stream macroinvertebrate sample are influenced by the total number of organisms in the sample. Higher total counts make it more likely to find rarer species in the sample. In order to compare richness and diversity across all samples, the R-program for calculating the Chessie BIBI index performs a probabilistic rarefaction and normalizes all sample counts to 100 organisms. The probabilistic rarefaction steps are described in Smith et al. (2017) page 7 and Appendix D and implemented in Step 1.5 of the Chessie BIBI program. As warned in the report and program documentation, the normalized sample created by the probabilistic rarefaction step is slightly different each time the program is run. This is because it randomly selects the rarer taxa to include in the BIBI calculation. If the program randomly selects a rare taxon that belongs to a commonly found family, taxa richness is unaffected because all genera are rolled-up to family level for the family-level, bioregion-specific Chessie BIBI index. If the program randomly picks a rare taxon that belongs to a rarely found family, the calculated richness and diversity metrics are higher. Richness and diversity are sensitive to disturbance and one or more of them are included in all of the twelve bioregion-specific indexes. A difference of one or two families in the count of a sample's families can change the final BIBI score and sometimes its rating.    

The following scripts will remove some of the uncertainty by determining which index rating (Excellent, Good, Fair, Poor, or Very Poor) has the highest probability of being the correct rating from repeated runs of the Chessie BIBI program.

**1.Function Creation**
Create a function that will let you repeat the Chessie BIBI program a set number of times. This function places all of the above scripts to calculate the metrics, scores, and ratings inside of a wrapper that can then be reapeatedly run in the next step.
```{r}
 ChessieBIBI_score <- function(){
      # Rarefaction
      taxa2.df <- taxa.df %>%                                                                      
        group_by(unique_id) %>%                                                                               
        mutate(sample_sum = sum(reporting_value)) %>%                                                         
        mutate(ratio = 100*(reporting_value/sample_sum)) %>%                                                  
        mutate(ratio_int = ifelse(ratio >= 1, round(ratio), as.integer(ratio))) %>%                           
        mutate(cumulative = ifelse(ratio_int >= 1, sum(ratio_int), 0)) %>%                                    
        mutate(rare_needed = as.integer(ifelse(ratio_int >= 1, 100-cumulative, 100-(max(cumulative))))) %>%   
        mutate(count= n()) %>%
        mutate(random = sample(100, size = n(), replace = FALSE)) %>%                                         
        arrange(ratio_int, random) %>%                                                                        
        ungroup()                                                                                             
       
       taxa3.df <- taxa2.df %>%
         group_by(unique_id) %>%
         mutate(rare_num = case_when(ratio_int == 0 ~ TRUE, ratio_int != 0 ~ FALSE)) %>%                      
         summarise(total_rare = sum(rare_num)) %>%                                                            
         ungroup()
       
       taxa_join <- left_join(taxa2.df, taxa3.df, by = "unique_id") %>%                            
         group_by(unique_id) %>%
         mutate(rare_needed = ifelse(rare_needed > total_rare, total_rare, rare_needed)) %>%                   
         mutate(rare_needed = ifelse(rare_needed < 1, 0, rare_needed)) %>%                                   
         mutate(rare_2b_removed = (total_rare - rare_needed)+1) %>%                                          
         slice(median(rare_2b_removed):n()) %>%                                                                       
         mutate(rare_count =(ifelse(ratio_int == 0, 1, ratio_int))) %>%                                       
         mutate(rare_count_total = sum(rare_count)) %>%                                                        
         ungroup()
         
       taxa.df4 <- taxa_join %>%                                                                               
        select(unique_id, icprb_bioregion_id, rare_count, phylum, subphylum, class, subclass, order, suborder, family, subfamily, tribe, genus, habit, species, ffg, bibi_tv)
      
       # Create a data frame of the non-rarefied and rarefied data.
      non_rare <- taxa.df %>%
      arrange(unique_id) %>%
      ungroup()

      rare <- taxa.df4 %>%
      arrange(unique_id) %>%
      ungroup()

# Create the metrics.key data frame.
metrics.key <- taxa.df %>%              
 select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time) %>%
 unique() %>%
 arrange(unique_id) 

###--------------------------------------------------------------------------
# ESTABLISH FUNCTIONS USED TO CALCULATE SOME METRICS
###--------------------------------------------------------------------------
original_order <- function(mod.df, original.df, unique.id.col) {
  unique.id.col <- rlang::enquo(unique.id.col)

  final.df <- original.df %>%
    dplyr::select(!!unique.id.col) %>%
    dplyr::distinct() %>%
    dplyr::left_join(mod.df, by = rlang::quo_name(unique.id.col))

  return(final.df)
}
###--------------------------------------------------------------------------
## CALCULATE THE METRICS ##
# Calculate how many individuals are in each unique_id (unique sampling event)
sample_sum <- non_rare %>%
 group_by(unique_id) %>%
 summarise(sample_sum = sum(reporting_value)) 

#
 aspt_mod <- function(long.df, unique.id.col, count.col, tol.col) {
  unique.id.col <- rlang::enquo(unique.id.col)
  tol.col <- rlang::enquo(tol.col )
  count.col <- rlang::enquo(count.col )
  #----------------------------------------------------------------------------
  long.df <- long.df %>%
    dplyr::filter((!!count.col) > 0)
  #----------------------------------------------------------------------------
 score.vec <- long.df %>%
    filter(!is.na(!!!tol.col),
           !rlang::UQ(tol.col) %in% c("")) %>%
    select(!!unique.id.col, !!tol.col, !!count.col) %>%
    mutate(score = rlang::UQ(count.col) * rlang::UQ(tol.col)) %>%
    group_by(!!unique.id.col) %>%
    summarize(score = sum(score),
              taxa = sum(!!count.col),
              score = score / taxa) %>%
    original_order(long.df, !!unique.id.col) %>%
    pull(score)

  return(score.vec)
}

metrics.key$aspt_mod <- aspt_mod(non_rare, unique_id, reporting_value, aspt)

# pct_gastro_oligo 
pct_gastro_oligo <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(class)) %>% 
 summarise(sum = sum(reporting_value[class %in% c("gastropoda", "clitellata")]))  
 metrics.key$pct_gastro_oligo <- (pct_gastro_oligo$sum / sample_sum$sample_sum )*100

# pct_diptera
pct_diptera <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>% 
 summarise(sum = sum(reporting_value[order %in% c("diptera")]))  
 metrics.key$pct_diptera <- (pct_diptera$sum / sample_sum$sample_sum)*100

# GOLD  
metrics.key$gold <- 100 - (metrics.key$pct_gastro_oligo + metrics.key$pct_diptera)

# HBI    
hbi_multiplier <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>% 
 summarise(multiplier = sum(reporting_value*bibi_tv)) 

hbi_divisor <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>% 
 summarise(divisor = sum(reporting_value))

metrics.key$hbi <- hbi_multiplier$multiplier / hbi_divisor$divisor

# pct_arthropoda      
pct_arthropoda <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(phylum)) %>% 
 summarise(sum = sum(reporting_value[phylum %in% c("arthropoda")]))  
 metrics.key$pct_arthropoda <- (pct_arthropoda$sum / sample_sum$sample_sum )*100

# pct_burrow    
pct_burrow <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(habit)) %>% 
 summarise(sum = sum(reporting_value[habit %in% c("burrow")]))  
 metrics.key$pct_burrow <- (pct_burrow$sum / sample_sum$sample_sum )*100
 
# pct_chironomidae   
pct_chironomidae <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(family)) %>% 
 summarise(sum = sum(reporting_value[family %in% c("chironomidae")]))  
 metrics.key$pct_chironomidae <- (pct_chironomidae$sum / sample_sum$sample_sum )*100

# pct_cling  
pct_cling <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(habit)) %>% 
 summarise(sum = sum(reporting_value[habit %in% c("cling")]))  
 metrics.key$pct_cling <- (pct_cling$sum / sample_sum$sample_sum )*100

#pct_collect 
pct_collect <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>% 
 summarise(sum = sum(reporting_value[ffg %in% c("filter", "gather")]))  
 metrics.key$pct_collect <- (pct_collect$sum / sample_sum$sample_sum )*100 

# pct_cote    
pct_cote <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>% 
 summarise(sum = sum(reporting_value[order %in% c("coleoptera", "trichoptera", "ephemeroptera", "odonata")]))  
 metrics.key$pct_cote <- (pct_cote$sum / sample_sum$sample_sum )*100

# pct_ephemeroptera    
pct_ephemeroptera <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>% 
 summarise(sum = sum(reporting_value[order %in% c("ephemeroptera")]))  
 metrics.key$pct_ephemeroptera <- (pct_ephemeroptera$sum / sample_sum$sample_sum )*100

# pct_baetidae    
pct_baetidae <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(family)) %>% 
 summarise(sum = sum(reporting_value[family %in% c("baetidae")]))  
 metrics.key$pct_baetidae <- (pct_baetidae$sum / sample_sum$sample_sum )*100
 
# pct_ephemeroptera_no_baetid       
pct_ephemeroptera_no_baetid <- pct_baetidae
pct_ephemeroptera_no_baetid$sum_ephemeroptera <- pct_ephemeroptera$sum 
pct_ephemeroptera_no_baetid$ephem_no_baetidae <- pct_ephemeroptera_no_baetid$sum_ephemeroptera - pct_ephemeroptera_no_baetid$sum
metrics.key$pct_ephemeroptera_no_baetid <- (pct_ephemeroptera_no_baetid$ephem_no_baetidae/sample_sum$sample_sum)*100

# pct_ept  
pct_ept <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>% 
 summarise(sum = sum(reporting_value[order %in% c("ephemeroptera", "plecoptera", "trichoptera")]))  
 metrics.key$pct_ept <- (pct_ept$sum / sample_sum$sample_sum )*100
  
# pct_hydropsychidae    
pct_hydropsychidae <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(family)) %>% 
 summarise(sum = sum(reporting_value[family %in% c("hydropsychidae")]))  
 metrics.key$pct_hydropsychidae <- (pct_hydropsychidae$sum / sample_sum$sample_sum )*100

# pct_ept_no_hydro    
pct_ept_hydro <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>% 
 summarise(sum = sum(reporting_value[family %in% c("hydropsychidae")]))
 pct_ept_hydro$subtract <- pct_ept$sum - pct_ept_hydro$sum
 metrics.key$pct_ept_no_hydro <- (pct_ept_hydro$subtract/pct_ept$sum)*100

# pct_euholognatha  
pct_euholognatha <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(suborder)) %>% 
 summarise(sum = sum(reporting_value[suborder %in% c("euholognatha")]))
 metrics.key$pct_euholognatha <- (pct_euholognatha$sum / sample_sum$sample_sum )*100
  
# pct_filter    
pct_filter <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>% 
 summarise(sum = sum(reporting_value[ffg %in% c("filter")]))  
 metrics.key$pct_filter <- (pct_filter$sum / sample_sum$sample_sum )*100

# pct_heptageniidae    
pct_heptageniidae <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(family)) %>% 
 summarise(sum = sum(reporting_value[family %in% c("heptageniidae")]))  
 metrics.key$pct_heptageniidae <- (pct_heptageniidae$sum / sample_sum$sample_sum )*100

# pct_hexapoda    
pct_hexapoda <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(subphylum)) %>% 
 summarise(sum = sum(reporting_value[subphylum %in% c("hexapoda")]))  
 metrics.key$pct_hexapoda <- (pct_hexapoda$sum / sample_sum$sample_sum )*100

# pct_hydro_ept    
metrics.key$pct_hydro_ept <- ifelse(metrics.key$pct_ept == 0, 0, (metrics.key$pct_hydropsychidae / metrics.key$pct_ept) * 100)

# pct_intol_0_3    
pct_intol_0_3 <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>% 
 summarise(sum = sum(reporting_value[bibi_tv %in% c("0", "1", "2", "3")]))  
 metrics.key$pct_intol_0_3 <- (pct_intol_0_3$sum / sample_sum$sample_sum)*100
  
# pct_intol_0_4   
pct_intol_0_4 <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>% 
 summarise(sum = sum(reporting_value[bibi_tv %in% c("0", "1", "2", "3", "4")]))  
 metrics.key$pct_intol_0_4 <- (pct_intol_0_4$sum / sample_sum$sample_sum)*100

# pct_amph_iso  
pct_amph_iso <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>% 
 summarise(sum = sum(reporting_value[order %in% c("amphipoda", "isopoda")]))  
 metrics.key$pct_amph_iso <- (pct_amph_iso$sum / sample_sum$sample_sum )*100  

# pct_ephemerellidae  
pct_ephemerellidae <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(family)) %>% 
 summarise(sum = sum(reporting_value[family %in% c("ephemerellidae")]))  
 metrics.key$pct_ephemerellidae <- (pct_ephemerellidae$sum / sample_sum$sample_sum )*100 

# pct_limestone    
metrics.key$pct_limestone <- metrics.key$pct_amph_iso + metrics.key$pct_ephemerellidae

# pct_mod_tol_4_6 
pct_mod_tol_4_6 <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>% 
 summarise(sum = sum(reporting_value[bibi_tv %in% c("4", "5", "6")]))  
 metrics.key$pct_mod_tol_4_6 <- (pct_mod_tol_4_6$sum / sample_sum$sample_sum)*100

# pct_trichoptera    
pct_trichoptera <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>% 
 summarise(sum = sum(reporting_value[order %in% c("trichoptera")]))  
 metrics.key$pct_trichoptera <- (pct_trichoptera$sum / sample_sum$sample_sum )*100

# pct_non_hydrop_trichoptera  
pct_non_hydrop_trichoptera <- pct_trichoptera
pct_non_hydrop_trichoptera$hydropsychidae <- pct_hydropsychidae$sum
pct_non_hydrop_trichoptera$subtract <- pct_non_hydrop_trichoptera$sum - pct_non_hydrop_trichoptera$hydropsychidae
metrics.key$pct_non_hydrop_trichoptera <- (pct_non_hydrop_trichoptera$subtract/pct_non_hydrop_trichoptera$sum)*100

# pct_odonata  
pct_odonata <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>% 
 summarise(sum = sum(reporting_value[order %in% c("odonata")]))  
 metrics.key$pct_odonata <- (pct_odonata$sum / sample_sum$sample_sum )*100
  
# pct_pisciforma  
pct_pisciforma <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(suborder)) %>% 
 summarise(sum = sum(reporting_value[suborder %in% c("pisciforma")]))  
 metrics.key$pct_pisciforma <- (pct_pisciforma$sum / sample_sum$sample_sum )*100

# pct_predator  
pct_predator <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>% 
 summarise(sum = sum(reporting_value[ffg %in% c("predator")]))  
 metrics.key$pct_predator <- (pct_predator$sum / sample_sum$sample_sum )*100 
  
# pct_pterygota  
pct_pterygota <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(subclass)) %>% 
 summarise(sum = sum(reporting_value[subclass %in% c("pterygota")]))  
 metrics.key$pct_pterygota <- (pct_pterygota$sum / sample_sum$sample_sum )*100

# pct_retreat_caddisfly  
pct_retreat_caddisfly <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(suborder)) %>% 
 summarise(sum = sum(reporting_value[suborder %in% c("annulipalpia")]))  
 metrics.key$pct_retreat_caddisfly <- (pct_retreat_caddisfly$sum / sample_sum$sample_sum )*100

# pct_scrape   
pct_scrape <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>% 
 summarise(sum = sum(reporting_value[ffg %in% c("scrape")]))  
 metrics.key$pct_scrape <- (pct_scrape$sum / sample_sum$sample_sum )*100
   
# pct_sprawl   
pct_sprawl <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(habit)) %>% 
 summarise(sum = sum(reporting_value[habit %in% c("sprawl")]))  
 metrics.key$pct_sprawl <- (pct_sprawl$sum / sample_sum$sample_sum )*100

# pct_systellognatha   
pct_systellognatha <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(suborder)) %>% 
 summarise(sum = sum(reporting_value[suborder %in% c("systellognatha")]))  
 metrics.key$pct_systellognatha <- (pct_systellognatha$sum / sample_sum$sample_sum )*100

# pct_tolerant_7_10   
pct_tolerant_7_10 <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>% 
 summarise(sum = sum(reporting_value[bibi_tv %in% c("7", "8", "9", "10")]))  
 metrics.key$pct_tolerant_7_10 <- (pct_tolerant_7_10$sum / sample_sum$sample_sum)*100

# pct_urban_intol   
pct_urban_intol <- non_rare %>%
 group_by(unique_id) %>%
 filter(!is.na(family)) %>% 
 summarise(sum = sum(reporting_value[family %in% c(
    "aeshnidae",
    "ameletidae",
    "asellidae",
    "athericidae",
    "baetidae",
    "brachycentridae",
    "caenidae",
    "calamoceratidae",
    "cambaridae",
    "capniidae",
    "ceratopogonidae",
    "chironomidae",
    "chloroperlidae",
    "cordulegastridae",
    "corduliidae",
    "corydalidae",
    "crangonyctidae",
    "elmidae",
    "ephemerellidae",
    "ephemeridae",
    "glossosomatidae",
    "gomphidae",
    "heptageniidae",
    "hydropsychidae",
    "hydroptilidae",
    "isonychiidae",
    "lepidostomatidae",
    "leptophlebiidae",
    "leuctridae",
    "libellulidae",
    "metretopodidae",
    "nemouridae",
    "odontoceridae",
    "peltoperlidae",
    "perlidae",
    "perlodidae",
    "philopotamidae",
    "phryganeidae",
    "polycentropodidae",
    "polymitarcyidae",
    "potamanthidae",
    "psephenidae",
    "pteronarcyidae",
    "rhyacophilidae",
    "sericostomatidae",
    "sialidae",
    "simuliidae",
    "stratiomyidae",
    "tabanidae",
    "taeniopterygidae",
    "tipulidae",
    "uenoidae",
    "viviparidae")]))  
  metrics.key$pct_urban_intol <- (pct_urban_intol$sum / sample_sum$sample_sum )*100

# Richness and diversity metrics require rarefied data
# Calculate family-level taxa richness for each unique_id
family_richness <- rare %>%
 group_by(unique_id) %>%
 summarise(family_richness = n_distinct(family)) 

# Margalefs
margalefs_numerator <- family_richness %>%
 mutate(family_richness = family_richness - 1)

margalefs_denominator <- rare %>%
 group_by(unique_id) %>%
 summarise(sample_sum = log(sum(rare_count)))

metrics.key$margalefs <- margalefs_numerator$family_richness / margalefs_denominator$sample_sum

# pct_ept_rich
pct_ept_rich <- rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>%
 filter(order %in% c("ephemeroptera", "plecoptera", "trichoptera")) %>%
 summarise(rich_ept = n_distinct(family))
 pct_ept_rich2 <- merge(pct_ept_rich, family_richness, by= "unique_id", all=TRUE)
 pct_ept_rich2[is.na(pct_ept_rich2)] <- 0
 metrics.key$pct_ept_rich <- (pct_ept_rich2$rich_ept / pct_ept_rich2$family_richness)*100
    
# pct_ept_rich_no_tol    
pct_ept_rich_no_tol <- rare %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>%
 filter(order %in% c("ephemeroptera", "plecoptera", "trichoptera")) %>%
 filter(bibi_tv %in% c("0","1","2","3","4","5","6")) %>%
 summarise(rich_ept_no_tol = n_distinct(family))
 pct_ept_rich_no_tol2 <- merge(pct_ept_rich_no_tol, family_richness, by= "unique_id", all=TRUE)
 pct_ept_rich_no_tol2[is.na(pct_ept_rich_no_tol2)] <- 0
 metrics.key$pct_ept_rich_no_tol <- (pct_ept_rich_no_tol2$rich_ept_no_tol / pct_ept_rich_no_tol2$family_richness)*100

# Pielou
family_richness_pielou <- rare %>%
  group_by(unique_id) %>%
  summarise(family_richness = n_distinct(family)) %>%
  mutate(ln_richness = log(family_richness)) # in R, log is natural log

sample_sum_rare <- rare %>%
  group_by(unique_id) %>%
  summarise(sample_sum = sum(rare_count)) 
sample_sum_join <- sample_sum_rare %>%
  left_join(rare, sample_sum_rare, by = "unique_id") %>%
  select(unique_id, sample_sum)

pielou_shannon <- rare %>%
  group_by(unique_id, family) %>%
  summarise(taxa_sum = sum(rare_count)) %>%
  left_join(sample_sum_join, pielou_shannon, by = "unique_id") %>%
  distinct() %>%
  mutate(Pi = taxa_sum / sample_sum) %>%
  mutate(Pi_lnPi = -(Pi * log(Pi))) %>%
  summarise(Shannon = sum(Pi_lnPi))

metrics.key$pielou <- pielou_shannon$Shannon / family_richness_pielou$ln_richness

# rich_burrow  
rich_burrow <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(habit)) %>%
 filter(habit %in% c("burrow")) %>%
 summarise(rich_burrow = n_distinct(family)) %>%
 select(unique_id, rich_burrow) 
 metrics.key <- merge(metrics.key, rich_burrow, by="unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0 
  
# rich_climb
rich_climb <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(habit)) %>%
 filter(habit %in% c("climb")) %>%
 summarise(rich_climb = n_distinct(family)) %>%
 select(unique_id, rich_climb) 
 metrics.key <- merge(metrics.key, rich_climb, by="unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0 
  
# rich_cling 
rich_cling <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(habit)) %>%
 filter(habit %in% c("cling")) %>%
 summarise(rich_cling = n_distinct(family)) %>%
 select(unique_id, rich_cling) 
 metrics.key <- merge(metrics.key, rich_cling, by="unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0 
 
# rich_collect   
rich_collect_filter <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>%
 filter(ffg %in% c("filter")) %>%
 summarise(rich_collect_filter = n_distinct(family)) %>%
 select(unique_id, rich_collect_filter)
 
rich_collect_gather <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>%
 filter(ffg %in% c("gather")) %>%
 summarise(rich_collect_gather = n_distinct(family)) %>%
 select(unique_id, rich_collect_gather)
 rich_collect <- merge(rich_collect_filter, rich_collect_gather, by = "unique_id", all=TRUE)
 rich_collect[is.na(rich_collect)] <- 0  
 rich_collect$rich_collect <- rich_collect$rich_collect_filter + rich_collect$rich_collect_gather
 rich_collect$rich_collect_gather <- NULL
 rich_collect$rich_collect_filter <- NULL
 metrics.key <- merge(metrics.key, rich_collect, by = "unique_id", all=TRUE)
   
# rich_ephemeroptera     
rich_ephemeroptera <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>%
 filter(order %in% c("ephemeroptera")) %>%
 summarise(rich_ephemeroptera = n_distinct(family)) %>%
 select(unique_id, rich_ephemeroptera) 
 metrics.key <- merge(metrics.key, rich_ephemeroptera, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0
     
# rich_ept    
rich_ept <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>%
 filter(order %in% c("ephemeroptera", "plecoptera", "trichoptera")) %>%
 summarise(rich_ept = n_distinct(family)) %>%
 select(unique_id, rich_ept) 
 metrics.key <- merge(metrics.key, rich_ept, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0   
 
# rich_filter    
rich_filter <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>%
 filter(ffg %in% c("filter")) %>%
 summarise(rich_filter = n_distinct(family)) %>%
 select(unique_id, rich_filter) 
 metrics.key <- merge(metrics.key, rich_filter, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0   

# rich_gather    
rich_gather <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>%
 filter(ffg %in% c("gather")) %>%
 summarise(rich_gather = n_distinct(family)) %>%
 select(unique_id, rich_gather) 
 metrics.key <- merge(metrics.key, rich_gather, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0   

# rich_intol 
rich_intol0 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("0")) %>%
 summarise(rich_intol0 = n_distinct(family)) %>%
 select(unique_id, rich_intol0)

rich_intol1 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("1")) %>%
 summarise(rich_intol1 = n_distinct(family)) %>%
 select(unique_id, rich_intol1)

rich_intol2 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("2")) %>%
 summarise(rich_intol2 = n_distinct(family)) %>%
 select(unique_id, rich_intol2)

rich_intol3 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("3")) %>%
 summarise(rich_intol3 = n_distinct(family)) %>%
 select(unique_id, rich_intol3)
 rich_intol01 <- merge(rich_intol0, rich_intol1, by= "unique_id", all=TRUE)
 rich_intol01[is.na(rich_intol01)] <- 0
 rich_intol012 <- merge(rich_intol01, rich_intol2, by= "unique_id", all=TRUE)
 rich_intol012[is.na(rich_intol012)] <- 0
 rich_intol0123 <- merge(rich_intol012, rich_intol3, by= "unique_id", all=TRUE)
 rich_intol0123[is.na(rich_intol0123)] <- 0
 rich_intol0123$rich_intol <- rich_intol0123$rich_intol0 + rich_intol0123$rich_intol1 + rich_intol0123$rich_intol2 + rich_intol0123$rich_intol3
 
rich_intol0123 <- rich_intol0123 %>%
 select(unique_id, rich_intol)
 metrics.key <- merge(metrics.key, rich_intol0123, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0   

# rich_modtol  
rich_modtol4 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("4")) %>%
 summarise(rich_modtol4 = n_distinct(family)) %>%
 select(unique_id, rich_modtol4)

rich_modtol5 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("5")) %>%
 summarise(rich_modtol5 = n_distinct(family)) %>%
 select(unique_id, rich_modtol5)

rich_modtol6 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("6")) %>%
 summarise(rich_modtol6 = n_distinct(family)) %>%
 select(unique_id, rich_modtol6)
 rich_modtol45 <- merge(rich_modtol4, rich_modtol5, by= "unique_id", all=TRUE)
 rich_modtol45[is.na(rich_modtol45)] <- 0
 rich_modtol456 <- merge(rich_modtol45, rich_modtol6, by= "unique_id", all=TRUE)
 rich_modtol456[is.na(rich_modtol456)] <- 0
 rich_modtol456$rich_modtol <- rich_modtol456$rich_modtol4 + rich_modtol456$rich_modtol5 + rich_modtol456$rich_modtol6
 
 rich_modtol456 <- rich_modtol456 %>%
 select(unique_id, rich_modtol)
 metrics.key <- merge(metrics.key, rich_modtol456, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0  
  
# rich_plecoptera  
rich_plecoptera <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>%
 filter(order %in% c("plecoptera")) %>%
 summarise(rich_plecoptera = n_distinct(family)) %>%
 select(unique_id, rich_plecoptera) 
 metrics.key <- merge(metrics.key, rich_plecoptera, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0  

# rich_predator    
rich_predator <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>%
 filter(ffg %in% c("predator")) %>%
 summarise(rich_predator = n_distinct(family)) %>%
 select(unique_id, rich_predator) 
 metrics.key <- merge(metrics.key, rich_predator, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0 
 
# rich_shred    
rich_shred <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(ffg)) %>%
 filter(ffg %in% c("shred")) %>%
 summarise(rich_shred = n_distinct(family)) %>%
 select(unique_id, rich_shred) 
 metrics.key <- merge(metrics.key, rich_shred, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0 

# rich_sprawl    
rich_sprawl <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(habit)) %>%
 filter(habit %in% c("sprawl")) %>%
 summarise(rich_sprawl = n_distinct(family)) %>%
 select(unique_id, rich_sprawl) 
 metrics.key <- merge(metrics.key, rich_sprawl, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0 
    
# rich_tol    
rich_tol7 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("7")) %>%
 summarise(rich_tol7 = n_distinct(family)) %>%
 select(unique_id, rich_tol7)

rich_tol8 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("8")) %>%
 summarise(rich_tol8 = n_distinct(family)) %>%
 select(unique_id, rich_tol8)

rich_tol9 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("9")) %>%
 summarise(rich_tol9 = n_distinct(family)) %>%
 select(unique_id, rich_tol9)

rich_tol10 <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(bibi_tv)) %>%
 filter(bibi_tv %in% c("10")) %>%
 summarise(rich_tol10 = n_distinct(family)) %>%
 select(unique_id, rich_tol10)
 rich_tol78 <- merge(rich_tol7, rich_tol8, by= "unique_id", all=TRUE)
 rich_tol78[is.na(rich_tol78)] <- 0
 rich_tol789 <- merge(rich_tol78, rich_tol9, by= "unique_id", all=TRUE)
 rich_tol789[is.na(rich_tol789)] <- 0
 rich_tol78910 <- merge(rich_tol789, rich_tol10, by= "unique_id", all=TRUE)
 rich_tol78910[is.na(rich_tol78910)] <- 0
 rich_tol78910$rich_tol <- rich_tol78910$rich_tol7 + rich_tol78910$rich_tol8 + rich_tol78910$rich_tol9 + rich_tol78910$rich_tol10
 
rich_tol78910 <- rich_tol78910 %>%
 select(unique_id, rich_tol)
 metrics.key <- merge(metrics.key, rich_tol78910, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0   
    
# rich_trichoptera  
rich_trichoptera <- rare  %>%
 group_by(unique_id) %>%
 filter(!is.na(order)) %>%
 filter(order %in% c("trichoptera")) %>%
 summarise(rich_trichoptera = n_distinct(family)) %>%
 select(unique_id, rich_trichoptera) 
 metrics.key <- merge(metrics.key, rich_trichoptera, by = "unique_id", all=TRUE)
 metrics.key[is.na(metrics.key)] <- 0

# Subset data by bioregion
NAPU <- metrics.key %>% filter(icprb_bioregion_id == "napu")
NCA <- metrics.key %>% filter(icprb_bioregion_id == "nca")
CA <- metrics.key %>% filter(icprb_bioregion_id == "ca")
MAC<- metrics.key %>% filter(icprb_bioregion_id == "mac")
SEP <- metrics.key %>% filter(icprb_bioregion_id == "sep")
BLUE <- metrics.key %>% filter(icprb_bioregion_id == "blue")
NRV <- metrics.key %>% filter(icprb_bioregion_id == "nrv")
SRV <- metrics.key %>% filter(icprb_bioregion_id == "srv")
UNP <- metrics.key %>% filter(icprb_bioregion_id == "unp")
SGV <- metrics.key %>% filter(icprb_bioregion_id == "sgv")
LNP <- metrics.key %>% filter(icprb_bioregion_id == "lnp")
PIED <- metrics.key %>% filter(icprb_bioregion_id == "pied")

#Select metrics for each bioregion
#Northern Appalachian Plateau and Uplands (NAPU)
NAPU_score <- NAPU %>% 
  select(
  "unique_id",
  "icprb_bioregion_id",
  "station_id",
  "sample_date",
  "sample_time",
  "margalefs", 
  "pct_cling",
  "pct_collect",
  "pct_ephemeroptera",
  "pct_ept",
  "pct_ept_no_hydro",
  "pct_ept_rich_no_tol",
  "pct_filter",
  "pct_intol_0_3",
  "pct_non_hydrop_trichoptera",
  "pct_scrape",
  "pielou",
  "rich_burrow",
  "rich_cling",
  "rich_gather",
  "rich_modtol",
  "rich_plecoptera",
  "rich_shred",
  "rich_sprawl")
       
# North Central Appalachians (NCA)
NCA_score <- NCA %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "aspt_mod",
    "pct_cote",
    "pct_ephemeroptera",
    "pct_heptageniidae",
    "pct_limestone",
    "rich_cling",
    "rich_ephemeroptera",
    "rich_ept",
    "rich_plecoptera")

# Central Appalachians (CA)
CA_score <- CA %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "pct_ept_rich",
    "pct_euholognatha",
    "pct_heptageniidae",
    "pct_hydro_ept",
    "pct_mod_tol_4_6",
    "pct_non_hydrop_trichoptera",
    "pct_pisciforma",
    "pct_systellognatha",
    "rich_burrow",
    "rich_collect",
    "rich_filter",
    "rich_plecoptera",
    "rich_predator")

# Middle Atlantic Coastal Plain (MAC)
MAC_score <- MAC %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "aspt_mod",
    "gold",
    "pct_arthropoda",
    "pct_chironomidae",
    "pct_collect",
    "pct_mod_tol_4_6",
    "pct_odonata",
    "pct_predator",
    "pct_scrape",
    "pct_urban_intol",
    "rich_climb",
    "rich_tol")

# Southeastern Plains (SEP)
SEP_score <- SEP %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "aspt_mod",
    "gold",
    "hbi",
    "pct_cling",
    "pct_cote",
    "pct_ept_rich_no_tol",
    "rich_cling",
    "rich_collect",
    "rich_ephemeroptera",
    "rich_filter",
    "rich_trichoptera")

# Blue Ridge (BLUE)
BLUE_score <- BLUE %>% 
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "pct_burrow",
    "pct_ept_rich_no_tol",
    "pct_scrape",
    "pct_systellognatha",
    "rich_ephemeroptera")

# Northern Ridge and Valley (NRV)
NRV_score <- NRV %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "aspt_mod",
    "pct_ephemeroptera_no_baetid",
    "pct_ept_rich",
    "rich_collect",
    "rich_ephemeroptera",
    "rich_predator",
    "rich_shred",
    "rich_trichoptera")

# Southern Ridge and Valley (SRV)
SRV_score <- SRV %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "aspt_mod",
    "pct_burrow",
    "pct_cling",
    "pct_cote",
    "pct_ephemeroptera_no_baetid",
    "pct_ept_rich_no_tol",
    "pct_predator",
    "pct_scrape",
    "rich_cling",
    "rich_collect",
    "rich_ephemeroptera",
    "rich_intol",
    "rich_plecoptera",
    "rich_predator",
    "rich_trichoptera")

# Upper Northern Piedmont (UNP)
UNP_score <- UNP %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "hbi",
    "pct_collect",
    "pct_ept_rich_no_tol",
    "pct_filter",
    "pct_hexapoda",
    "pct_pisciforma",
    "pct_retreat_caddisfly",
    "pct_scrape",
    "pct_urban_intol",
    "rich_ephemeroptera",
    "rich_plecoptera",
    "rich_trichoptera")

# Southern Great Valley (SGV)
SGV_score <- SGV %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "aspt_mod",
    "gold",
    "pct_cling",
    "pct_collect",
    "pct_ephemeroptera_no_baetid",
    "pct_ept_no_hydro",
    "pct_ept_rich_no_tol",
    "pct_filter",
    "pct_heptageniidae",
    "pct_mod_tol_4_6",
    "pct_non_hydrop_trichoptera",
    "pct_pisciforma",
    "pct_predator",
    "pct_pterygota",
    "pct_scrape",
    "pct_sprawl",
    "pct_tolerant_7_10",
    "rich_collect",
    "rich_ephemeroptera",
    "rich_filter",
    "rich_trichoptera")

# Lower Northern Piedmont (LNP)
LNP_score <- LNP %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "aspt_mod",
    "pct_cling",
    "pct_collect",
    "pct_ept_rich",
    "pct_mod_tol_4_6",
    "pct_predator",
    "pct_pterygota",
    "rich_cling",
    "rich_ephemeroptera")

# Piedmont (PIED)
PIED_score <- PIED %>%
  select(
    "unique_id",
    "icprb_bioregion_id",
    "station_id",
    "sample_date",
    "sample_time",
    "pct_cling",
    "pct_ept_rich",
    "pct_intol_0_4",
    "pct_tolerant_7_10",
    "pct_urban_intol",
    "rich_cling",
    "rich_ephemeroptera",
    "rich_trichoptera")

# Score the metrics
# NAPU Metric Scoring
NAPU_score$score_margalefs <- 
  case_when(
    NAPU_score$margalefs > 2.83 ~ 100,
    NAPU_score$margalefs < 2.74 ~ 0,
    TRUE ~ ((NAPU_score$margalefs - 2.74)/(2.83-2.74))*100)

NAPU_score$score_pct_cling <- 
  case_when(
    NAPU_score$pct_cling > 59.19 ~ 100,
    NAPU_score$pct_cling < 50.98 ~ 0,
    TRUE ~ ((NAPU_score$pct_cling-50.98)/(59.19-50.98))*100)

NAPU_score$score_pct_collect <- 
   case_when(
     NAPU_score$pct_collect < 71.84 ~ 100,
     NAPU_score$pct_collect > 82.8 ~ 0,
     TRUE ~ ((82.8-NAPU_score$pct_collect)/(82.8-71.84))*100)

NAPU_score$score_pct_ephemeroptera <- 
  case_when(
    NAPU_score$pct_ephemeroptera > 31.47 ~ 100,
    NAPU_score$pct_ephemeroptera < 19.62 ~ 0,
    TRUE ~ ((NAPU_score$pct_ephemeroptera-19.62)/(31.47-19.62))*100)

NAPU_score$score_pct_ept <- 
  case_when(
    NAPU_score$pct_ept > 62.26 ~ 100,
    NAPU_score$pct_ept < 54.08 ~ 0,
    TRUE ~ ((NAPU_score$pct_ept-54.08)/(62.26-54.08))*100)

NAPU_score$score_pct_ept_no_hydro <- 
  case_when(
    NAPU_score$pct_ept_no_hydro > 87.96 ~ 100,
    NAPU_score$pct_ept_no_hydro < 70.51 ~ 0,
    TRUE ~ ((NAPU_score$pct_ept_no_hydro-70.51)/(87.96-70.51))*100)

NAPU_score$score_pct_ept_rich_no_tol <- 
  case_when(
    NAPU_score$pct_ept_rich_no_tol > 45.32 ~ 100,
    NAPU_score$pct_ept_rich_no_tol < 32.85 ~ 0,
    TRUE ~ ((NAPU_score$pct_ept_rich_no_tol-32.85)/(45.32-32.85))*100)

NAPU_score$score_pct_filter <- 
  case_when(
    NAPU_score$pct_filter < 15.51 ~ 100,
    NAPU_score$pct_filter > 22.62 ~ 0,
    TRUE ~ ((22.62-NAPU_score$pct_filter)/(22.62-15.51))*100)

NAPU_score$score_pct_intol_0_3 <- 
  case_when(
    NAPU_score$pct_intol_0_3 > 44.7 ~ 100,
    NAPU_score$pct_intol_0_3 < 24.22 ~ 0,
    TRUE ~ ((NAPU_score$pct_intol_0_3-24.22)/(44.7-24.22))*100)

NAPU_score$score_pct_non_hydrop_trichoptera <-
  case_when(
    NAPU_score$pct_non_hydrop_trichoptera > 45.08 ~ 100,
    NAPU_score$pct_non_hydrop_trichoptera < 20.6 ~ 0,
    TRUE ~ ((NAPU_score$pct_non_hydrop_trichoptera-20.6)/(45.08-20.6))*100)

NAPU_score$score_pct_scrape <- 
  case_when(
    NAPU_score$pct_scrape > 9.36 ~ 100,
    NAPU_score$pct_scrape < 6.89 ~ 0,
    TRUE ~ ((NAPU_score$pct_scrape-6.89)/(9.36-6.89))*100)

NAPU_score$score_pielou <- 
  case_when(
    NAPU_score$pielou < 0.79 ~ 100,
    NAPU_score$pielou > 0.80 ~ 0,
    TRUE ~ ((0.80-NAPU_score$pielou)/(0.80-0.79))*100)

NAPU_score$score_rich_burrow <- 
  case_when(
    NAPU_score$rich_burrow > 2.95 ~ 100,
    NAPU_score$rich_burrow < 1.93 ~ 0,
    TRUE ~ ((NAPU_score$rich_burrow-1.93)/(2.95-1.93))*100)

NAPU_score$score_rich_cling <- 
  case_when(
    NAPU_score$rich_cling > 8.0 ~ 100,
    NAPU_score$rich_cling < 6.77 ~ 0,
    TRUE ~ ((NAPU_score$rich_cling-6.77)/(8.0-6.77))*100)

NAPU_score$score_rich_gather <- 
  case_when(
    NAPU_score$rich_gather < 5.0 ~ 100,
    NAPU_score$rich_gather > 6.01 ~ 0,
    TRUE ~ ((6.01-NAPU_score$rich_gather)/(6.01-5.0))*100)

NAPU_score$score_rich_modtol <- 
  case_when(
     NAPU_score$rich_modtol < 6.0 ~ 100,
     NAPU_score$rich_modtol > 7.03 ~ 0,
     TRUE ~ ((7.03-NAPU_score$rich_modtol)/(7.03-6.0))*100)

NAPU_score$score_rich_plecoptera <- 
  case_when(
    NAPU_score$rich_plecoptera > 2.13 ~ 100,
    NAPU_score$rich_plecoptera < 0.70 ~ 0,
    TRUE ~ ((NAPU_score$rich_plecoptera-0.70)/(2.13-0.70))*100)

NAPU_score$score_rich_shred <- 
  case_when(
    NAPU_score$rich_shred > 2.0 ~ 100,
    NAPU_score$rich_shred < 0.84 ~ 0,
    TRUE ~ ((NAPU_score$rich_shred-0.84)/(2.0-0.84))*100)

NAPU_score$score_rich_sprawl <- 
  case_when(
    NAPU_score$rich_sprawl < 1.0 ~ 100,
    NAPU_score$rich_sprawl > 2.3 ~ 0,
    TRUE ~ ((2.3-NAPU_score$rich_sprawl)/(2.3-1.0))*100)

# NCA Metric Scoring
NCA_score$score_aspt_mod <-
  case_when(
    NCA_score$aspt_mod > 7.54 ~ 100,
    NCA_score$aspt_mod < 6.67 ~ 0,
    TRUE ~ ((NCA_score$aspt_mod - 6.67)/(7.54-6.67))*100)
  
NCA_score$score_pct_cote <-
  case_when(
    NCA_score$pct_cote > 63.95 ~ 100,
    NCA_score$pct_cote < 47.71 ~ 0,
    TRUE ~ ((NCA_score$pct_cote - 47.71)/(63.95-47.71))*100)

NCA_score$score_pct_ephemeroptera <-
  case_when(
    NCA_score$pct_ephemeroptera > 40.73 ~ 100,
    NCA_score$pct_ephemeroptera < 9.02 ~ 0,
    TRUE ~ ((NCA_score$pct_ephemeroptera - 9.02)/(40.73-9.02))*100)

NCA_score$score_pct_heptageniidae <-
  case_when(
   NCA_score$pct_heptageniidae > 13.72 ~ 100,
   NCA_score$pct_heptageniidae  < 0 ~ 0,
   TRUE ~ ((NCA_score$pct_heptageniidae - 0)/(13.72-0))*100)

NCA_score$score_pct_limestone <-
  case_when(
    NCA_score$pct_limestone > 8.93 ~ 100,
    NCA_score$pct_limestone < 0 ~ 0,
    TRUE ~ ((NCA_score$pct_limestone - 0)/(8.93-0))*100)

NCA_score$score_rich_cling <-
  case_when(
    NCA_score$rich_cling > 10.0 ~ 100,
    NCA_score$rich_cling < 7.75 ~ 0,
    TRUE ~ ((NCA_score$rich_cling-7.75)/(10.0-7.75))*100)

NCA_score$score_rich_ephemeroptera <- 
  case_when(
    NCA_score$rich_ephemeroptera > 4.0 ~ 100,
    NCA_score$rich_ephemeroptera < 3.83 ~ 0,
    TRUE ~ ((NCA_score$rich_ephemeroptera - 3.83)/(4.0-3.83))*100)

NCA_score$score_rich_ept <-
  case_when(
    NCA_score$rich_ept > 12.0 ~ 100,
    NCA_score$rich_ept < 8.02 ~ 0,
    TRUE ~ ((NCA_score$rich_ept-8.02)/(12.0-8.02))*100)

NCA_score$score_rich_plecoptera <-
  case_when(
    NCA_score$rich_plecoptera > 4.04 ~ 100,
    NCA_score$rich_plecoptera < 3.45 ~ 0,
    TRUE ~ ((NCA_score$rich_plecoptera - 3.45)/(4.04-3.45))*100)

# CA Metric Scoring
CA_score$score_pct_ept_rich <-
  case_when(
    CA_score$pct_ept_rich > 66.82 ~ 100,
    CA_score$pct_ept_rich < 56.82 ~ 0,
    TRUE ~ ((CA_score$pct_ept_rich-56.82)/(66.82-56.82))*100)

CA_score$score_pct_euholognatha <-
  case_when(
    CA_score$pct_euholognatha > 17.98 ~ 100,
    CA_score$pct_euholognatha < 0.23 ~ 0,
    TRUE ~ ((CA_score$pct_euholognatha-0.23)/(17.98-0.23))*100)

CA_score$score_pct_heptageniidae <-
  case_when(
    CA_score$pct_heptageniidae > 6.69 ~ 100,
    CA_score$pct_heptageniidae < 0 ~ 0,
    TRUE ~ ((CA_score$pct_heptageniidae-0)/(6.69-0))*100)

CA_score$score_pct_hydro_ept <-
  case_when(
    CA_score$pct_hydro_ept < 12.5 ~ 100,
    CA_score$pct_hydro_ept > 21.73 ~ 0,
    TRUE ~ ((21.73-CA_score$pct_hydro_ept)/(21.73-12.5))*100)

CA_score$score_pct_mod_tol_4_6 <- 
  case_when(
    CA_score$pct_mod_tol_4_6 < 43.04 ~ 100,
    CA_score$pct_mod_tol_4_6 > 57.26 ~ 0,
    TRUE ~ ((57.26-CA_score$pct_mod_tol_4_6)/(57.26-43.04))*100)

CA_score$score_pct_non_hydrop_trichoptera <-
  case_when(
    CA_score$pct_non_hydrop_trichoptera > 40.4 ~ 100,
    CA_score$pct_non_hydrop_trichoptera < 25.02 ~ 0,
    TRUE ~ ((CA_score$pct_non_hydrop_trichoptera - 25.02)/(40.4-25.02))*100)

CA_score$score_pct_pisciforma <-
  case_when(
    CA_score$pct_pisciforma > 12.33 ~ 100,
    CA_score$pct_pisciforma < 5.49 ~ 0,
    TRUE ~ ((CA_score$pct_pisciforma-5.49)/(12.33-5.49))*100)

CA_score$score_pct_systellognatha <-
  case_when(
    CA_score$pct_systellognatha > 5.44 ~ 100,
    CA_score$pct_systellognatha < 0.08 ~ 0,
    TRUE ~ ((CA_score$pct_systellognatha-0.08)/(5.44-0.08))*100)

CA_score$score_rich_burrow <- 
  case_when(
    CA_score$rich_burrow > 3.0 ~ 100,
    CA_score$rich_burrow < 1.4 ~ 0,
    TRUE ~ ((CA_score$rich_burrow-1.4)/(3.0-1.4))*100)

CA_score$score_rich_collect <- 
  case_when(
    CA_score$rich_collect > 7.0 ~ 100,
    CA_score$rich_collect < 5.13 ~ 0,
    TRUE ~ ((CA_score$rich_collect-5.13)/(7.0-5.13))*100)

CA_score$score_rich_filter <-
  case_when(
    CA_score$rich_filter > 2.91 ~ 100,
    CA_score$rich_filter < 1.38 ~ 0,
    TRUE ~ ((CA_score$rich_filter-1.38)/(2.91-1.38))*100)

CA_score$score_rich_plecoptera <-
  case_when(
    CA_score$rich_plecoptera > 4.0 ~ 100,
    CA_score$rich_plecoptera < 0.47 ~ 0,
    TRUE ~ ((CA_score$rich_plecoptera-0.47)/(4.0-0.47))*100)

CA_score$score_rich_predator <- 
  case_when(
    CA_score$rich_predator > 3.0 ~ 100,
    CA_score$rich_predator < 1.44 ~ 0,
    TRUE ~ ((CA_score$rich_predator-1.44)/(3.0-1.44))*100)

# MAC Metric Scoring
MAC_score$score_aspt_mod <-
  case_when(
    MAC_score$aspt_mod > 3.27 ~ 100,
    MAC_score$aspt_mod < 2.9 ~ 0,
    TRUE ~ ((MAC_score$aspt_mod-2.9)/(3.27-2.9))*100)

MAC_score$score_gold <-
  case_when(
    MAC_score$gold > 58.96 ~ 100,
    MAC_score$gold < 33.91 ~ 0,
    TRUE ~ ((MAC_score$gold-33.91)/(58.96-33.91))*100)

MAC_score$score_pct_arthropoda <-
  case_when(
    MAC_score$pct_arthropoda < 74.01 ~ 100,
    MAC_score$pct_arthropoda > 95.38 ~ 0,
    TRUE ~ ((95.38-MAC_score$pct_arthropoda)/(95.38-74.01))*100)

MAC_score$score_pct_chironomidae <-
  case_when(
    MAC_score$pct_chironomidae < 23.01 ~ 100,
    MAC_score$pct_chironomidae > 42.68 ~ 0,
    TRUE ~ ((42.68-MAC_score$pct_chironomidae)/(42.68-23.01))*100)
  
MAC_score$score_pct_collect <-
  case_when(
    MAC_score$pct_collect < 73.25 ~ 100,
    MAC_score$pct_collect > 96.19 ~ 0,
    TRUE ~ ((96.19 - MAC_score$pct_collect)/(96.19-73.25))*100)

MAC_score$score_pct_mod_tol_4_6 <-
  case_when(
    MAC_score$pct_mod_tol_4_6 < 52.49 ~ 100,
    MAC_score$pct_mod_tol_4_6 > 75.48 ~ 0,
    TRUE ~ ((75.48-MAC_score$pct_mod_tol_4_6)/(75.48-52.49))*100)

MAC_score$score_pct_odonata <-
  case_when(
    MAC_score$pct_odonata > 1.82 ~ 100,
    MAC_score$pct_odonata < 0.43 ~ 0,
    TRUE ~ ((MAC_score$pct_odonata - 0.43)/(1.82-0.43))*100)

MAC_score$score_pct_predator <-
  case_when(
    MAC_score$pct_predator > 5.16 ~ 100,
    MAC_score$pct_predator < 4.01 ~ 0,
    TRUE ~ ((MAC_score$pct_predator - 4.01)/(5.16-4.01))*100)

MAC_score$score_pct_scrape <-
  case_when(
    MAC_score$pct_scrape > 12.07 ~ 100,
    MAC_score$pct_scrape < 0 ~ 0,
    TRUE ~ ((MAC_score$pct_scrape-0)/(12.07-0))*100)

MAC_score$score_pct_urban_intol <-
  case_when(
    MAC_score$pct_urban_intol < 65.28 ~ 100,
    MAC_score$pct_urban_intol > 71.78 ~ 0,
    TRUE ~ ((71.78-MAC_score$pct_urban_intol)/(71.78-65.28))*100)

MAC_score$score_rich_climb <-
  case_when(
    MAC_score$rich_climb > 2.45 ~ 100,
    MAC_score$rich_climb < 0 ~ 0,
    TRUE ~ ((MAC_score$rich_climb-0)/(2.45-0))*100)

MAC_score$score_rich_tol <-
  case_when(
    MAC_score$rich_tol > 4.0 ~ 100,
    MAC_score$rich_tol < 2.16 ~ 0,
    TRUE ~ ((MAC_score$rich_tol-2.16)/(4.0-2.16))*100)

# SEP Metric Scoring
SEP_score$score_aspt_mod <-
  case_when(
    SEP_score$aspt_mod > 4.65 ~ 100,
    SEP_score$aspt_mod < 3.18 ~ 0,
    TRUE ~ ((SEP_score$aspt_mod-3.18)/(4.65-3.18))*100)

SEP_score$score_gold <-
  case_when(
    SEP_score$gold > 50.3 ~ 100,
    SEP_score$gold < 32.16 ~ 0,
    TRUE ~ ((SEP_score$gold-32.16)/(50.3-32.16))*100)

SEP_score$score_hbi <-
  case_when(
    SEP_score$hbi < 5.12 ~ 100,
    SEP_score$hbi > 5.86 ~ 0,
    TRUE ~ ((5.86-SEP_score$hbi)/(5.86-5.12))*100)

SEP_score$pct_cling <-
  case_when(
    SEP_score$pct_cling > 43.05 ~ 100,
    SEP_score$pct_cling < 0 ~ 0,
    TRUE ~ ((SEP_score$pct_cling-0)/(43.05-0))*100)

SEP_score$score_pct_cote <-
  case_when(
    SEP_score$pct_cote > 26.37 ~ 100,
    SEP_score$pct_cote < 1.62 ~ 0,
    TRUE ~ ((SEP_score$pct_cote-1.62)/(26.37-1.62))*100)

SEP_score$score_pct_ept_rich_no_tol <-
  case_when(
    SEP_score$pct_ept_rich_no_tol > 24.73 ~ 100,
    SEP_score$pct_ept_rich_no_tol < 16.49 ~ 0,
    TRUE ~ ((SEP_score$pct_ept_rich_no_tol-16.49)/(24.73-16.49))*100)

SEP_score$score_rich_cling <-
  case_when(
    SEP_score$rich_cling > 4.93 ~ 100,
    SEP_score$rich_cling < 1.03 ~ 0,
    TRUE ~ ((SEP_score$rich_cling-1.03)/(4.93-1.03))*100)

SEP_score$score_rich_collect <-
  case_when(
    SEP_score$rich_collect > 8.0 ~ 100,
    SEP_score$rich_collect < 4.62 ~ 0,
    TRUE ~ ((SEP_score$rich_collect-4.62)/(8.0-4.62))*100)

SEP_score$score_rich_ephemeroptera <-
  case_when(
    SEP_score$rich_ephemeroptera > 2.0 ~ 100,
    SEP_score$rich_ephemeroptera < 0 ~ 0,
    TRUE ~ ((SEP_score$rich_ephemeroptera-0)/(2.0-0))*100)

SEP_score$score_rich_filter <-
  case_when(
    SEP_score$rich_filter  > 2.31 ~ 100,
    SEP_score$rich_filter < 0.04 ~ 0,
    TRUE ~ ((SEP_score$rich_filter-0.04)/(2.31-0.04))*100)

SEP_score$score_rich_trichoptera <-
  case_when(
    SEP_score$rich_trichoptera > 2.0 ~ 100,
    SEP_score$rich_trichoptera < 0.41 ~ 0,
    TRUE ~ ((SEP_score$rich_trichoptera-0.41)/(2.0-0.41))*100)

# BLUE Metric Scoring
BLUE_score$score_pct_burrow <-
  case_when(
    BLUE_score$pct_burrow < 12.65 ~ 100,
    BLUE_score$pct_burrow > 40.71 ~ 0,
    TRUE ~ ((40.71-BLUE_score$pct_burrow)/(40.71-12.65))*100)

BLUE_score$score_pct_ept_rich_no_tol <-
  case_when(
    BLUE_score$pct_ept_rich_no_tol > 52.78 ~ 100,
    BLUE_score$pct_ept_rich_no_tol < 34.74 ~ 0,
    TRUE ~ ((BLUE_score$pct_ept_rich_no_tol-34.74)/(52.78-34.74))*100)

BLUE_score$score_pct_scrape <-
  case_when(
    BLUE_score$pct_scrape > 15.23 ~ 100,
    BLUE_score$pct_scrape < 0 ~ 0,
    TRUE ~ ((BLUE_score$pct_scrape-0)/(15.23-0))*100)

BLUE_score$score_pct_systellognatha <-
  case_when(
    BLUE_score$pct_systellognatha > 6.06 ~ 100,
    BLUE_score$pct_systellognatha < 0 ~ 0,
    TRUE ~ ((BLUE_score$pct_systellognatha-0)/(6.06-0))*100)

BLUE_score$score_rich_ephemeroptera <-
  case_when(
    BLUE_score$rich_ephemeroptera > 4.0 ~ 100,
    BLUE_score$rich_ephemeroptera < 0.33 ~ 0,
    TRUE ~ ((BLUE_score$rich_ephemeroptera-0.33)/(4.0-0.33))*100)

# NRV Metric Scoring
NRV_score$score_aspt_mod <-
  case_when(
    NRV_score$aspt_mod > 6.27 ~ 100,
    NRV_score$aspt_mod < 4.59 ~ 0,
    TRUE ~ ((NRV_score$aspt_mod-4.59)/(6.27-4.59))*100)

NRV_score$score_pct_ephemeroptera_no_baetid <-
  case_when(
    NRV_score$pct_ephemeroptera_no_baetid > 23.69 ~ 100,
    NRV_score$pct_ephemeroptera_no_baetid < 6.01 ~ 0,
    TRUE ~ ((NRV_score$pct_ephemeroptera_no_baetid-6.01)/(23.69-6.01))*100)

NRV_score$score_pct_ept_rich <-
  case_when(
    NRV_score$pct_ept_rich > 64.7 ~ 100,
    NRV_score$pct_ept_rich < 47.53 ~ 0,
    TRUE ~ ((NRV_score$pct_ept_rich-47.53)/(64.7-47.53))*100)

NRV_score$score_rich_collect <-
  case_when(
    NRV_score$rich_collect > 8.0 ~ 100,
    NRV_score$rich_collect < 7.14 ~ 0,
    TRUE ~ ((NRV_score$rich_collect-7.14)/(8.0-7.14))*100)

NRV_score$score_rich_ephemeroptera <-
  case_when(
    NRV_score$rich_ephemeroptera > 4.0 ~ 100,
    NRV_score$rich_ephemeroptera < 1.8 ~ 0,
    TRUE ~ ((NRV_score$rich_ephemeroptera-1.8)/(4.0-1.8))*100)

NRV_score$score_rich_predator <-
  case_when(
    NRV_score$rich_predator > 3.0 ~ 100,
    NRV_score$rich_predator < 0.32 ~ 0,
    TRUE ~ ((NRV_score$rich_predator-0.32)/(3.0-0.32))*100)

NRV_score$score_rich_shred <-
  case_when(
    NRV_score$rich_shred > 2.0 ~ 100,
    NRV_score$rich_shred < 1.27 ~ 0,
    TRUE ~ ((NRV_score$rich_shred-1.27)/(2.0-1.27))*100)

NRV_score$score_rich_trichoptera <-
  case_when(
    NRV_score$rich_trichoptera > 3.0 ~ 100,
    NRV_score$rich_trichoptera < 2.53 ~ 0,
    TRUE ~ ((NRV_score$rich_trichoptera-2.53)/(3.0-2.53))*100)

# UNP Metric Scoring
UNP_score$score_hbi <-
  case_when(
    UNP_score$hbi < 4.01 ~ 100,
    UNP_score$hbi > 5.74 ~ 0,
    TRUE ~ ((5.74-UNP_score$hbi)/(5.74-4.01))*100)

UNP_score$score_pct_collect <-
  case_when(
    UNP_score$pct_collect < 80.32 ~ 100,
    UNP_score$pct_collect > 90.97 ~ 0,
    TRUE ~ ((90.97-UNP_score$pct_collect)/(90.97-80.32))*100)

UNP_score$score_pct_ept_rich_no_tol <-
  case_when(
    UNP_score$pct_ept_rich_no_tol > 37.09 ~ 100,
    UNP_score$pct_ept_rich_no_tol < 22.1 ~ 0,
    TRUE ~ ((UNP_score$pct_ept_rich_no_tol-22.1)/(37.09-22.1))*100)

UNP_score$score_pct_filter <-
  case_when(
    UNP_score$pct_filter > 18.14 ~ 100,
    UNP_score$pct_filter < 1.91 ~ 0,
    TRUE ~ ((UNP_score$pct_filter - 1.91)/(18.14-1.91))*100)

UNP_score$score_pct_hexapoda <-
  case_when(
    UNP_score$pct_hexapoda > 99.5 ~ 100,
    UNP_score$pct_hexapoda < 96.59 ~ 0,
    TRUE ~ ((UNP_score$pct_hexapoda-96.59)/(99.5-96.59))*100)

UNP_score$score_pct_pisciforma <-
  case_when(
    UNP_score$pct_pisciforma > 9.43 ~ 100,
    UNP_score$pct_pisciforma < 0 ~ 0,
    TRUE ~ ((UNP_score$pct_pisciforma-0)/(9.43-0))*100)

UNP_score$score_pct_retreat_caddisfly <-
  case_when(
    UNP_score$pct_retreat_caddisfly > 10.97 ~ 100,
    UNP_score$pct_retreat_caddisfly < 2.89 ~ 0,
    TRUE ~ ((UNP_score$pct_retreat_caddisfly-2.89)/(10.97-2.89))*100)

UNP_score$score_pct_scrape <-
  case_when(
    UNP_score$pct_scrape > 5.72 ~ 100,
    UNP_score$pct_scrape < 0.22 ~ 0,
    TRUE ~ ((UNP_score$pct_scrape-0.22)/(5.72-0.22))*100)

UNP_score$score_pct_urban_intol <- 
  case_when(
    UNP_score$pct_urban_intol > 98.46 ~ 100,
    UNP_score$pct_urban_intol < 91.71 ~ 0,
    TRUE ~ ((UNP_score$pct_urban_intol-91.71)/(98.46-91.71))*100)

UNP_score$score_rich_ephemeroptera <-
  case_when(
    UNP_score$rich_ephemeroptera > 3.0 ~ 100,
    UNP_score$rich_ephemeroptera < 1.44 ~ 0,
    TRUE ~ ((UNP_score$rich_ephemeroptera-1.44)/(3.0-1.44))*100)

UNP_score$score_rich_plecoptera <-
  case_when(
    UNP_score$rich_plecoptera > 2.23 ~ 100,
    UNP_score$rich_plecoptera < 0 ~ 0,
    TRUE ~ ((UNP_score$rich_plecoptera-0)/(2.23-0))*100)

UNP_score$score_rich_trichoptera <-
  case_when(
    UNP_score$rich_trichoptera > 3.0 ~ 100,
    UNP_score$rich_trichoptera < 1.01 ~ 0,
    TRUE ~ ((UNP_score$rich_trichoptera-1.01)/(3.0-1.01))*100)

# SRV Metric Scoring
SRV_score$score_aspt_mod <-
  case_when(
    SRV_score$aspt_mod > 6.44 ~ 100,
    SRV_score$aspt_mod < 4.95 ~ 0,
    TRUE ~ ((SRV_score$aspt_mod-4.95)/(6.44-4.95))*100)

SRV_score$score_pct_burrow <-
  case_when(
    SRV_score$pct_burrow < 17.13 ~ 100,
    SRV_score$pct_burrow > 31.28 ~ 0,
    TRUE ~ ((31.28-SRV_score$pct_burrow)/(31.28-17.13))*100)

SRV_score$score_pct_cling <-
  case_when(
    SRV_score$pct_cling > 58.45 ~ 100,
    SRV_score$pct_cling < 40.16 ~ 100,
    TRUE ~ ((SRV_score$pct_cling-40.16)/(58.45-40.16))*100)

SRV_score$score_pct_cote <-
  case_when(
    SRV_score$pct_cote > 57.55 ~ 100,
    SRV_score$pct_cote < 43.78 ~ 0,
    TRUE ~ ((SRV_score$pct_cote-43.78)/(57.55-43.78))*100)

SRV_score$score_pct_ephemeroptera_no_baetid <-
  case_when(
    SRV_score$pct_ephemeroptera_no_baetid > 21.0 ~ 100,
    SRV_score$pct_ephemeroptera_no_baetid < 2.8 ~ 0,
    TRUE ~ ((SRV_score$pct_ephemeroptera_no_baetid-2.8)/(21.0-2.8))*100)

SRV_score$score_pct_ept_rich_no_tol <-
  case_when(
    SRV_score$pct_ept_rich_no_tol > 47.16 ~ 100,
    SRV_score$pct_ept_rich_no_tol < 34.94 ~ 0,
    TRUE ~ ((SRV_score$pct_ept_rich_no_tol-34.94)/(47.16-34.94))*100)

SRV_score$score_pct_predator <-
  case_when(
    SRV_score$pct_predator > 7.12 ~ 100,
    SRV_score$pct_predator < 2.24 ~ 0,
    TRUE ~ ((SRV_score$pct_predator-2.24)/(7.12-2.24))*100)

SRV_score$score_pct_scrape <-
  case_when(
    SRV_score$pct_scrape > 10.88 ~ 100,
    SRV_score$pct_scrape < 0.47 ~ 0,
    TRUE ~ ((SRV_score$pct_scrape-0.47)/(10.88-0.47))*100)

SRV_score$score_rich_cling <- 
  case_when(
    SRV_score$rich_cling > 8.43 ~ 100,
    SRV_score$rich_cling < 6.16 ~ 0,
    TRUE ~ ((SRV_score$rich_cling-6.16)/(8.43-6.16))*100)

SRV_score$score_rich_collect <-
  case_when(
    SRV_score$rich_collect > 7.97 ~ 100,
    SRV_score$rich_collect < 6.02 ~ 0,
    TRUE ~ ((SRV_score$rich_collect-6.02)/(7.97-6.02))*100)

SRV_score$score_rich_ephemeroptera <-
  case_when(
    SRV_score$rich_ephemeroptera > 4.0 ~ 100,
    SRV_score$rich_ephemeroptera < 3.66 ~ 0,
    TRUE ~ ((SRV_score$rich_ephemeroptera-3.66)/(4.0-3.66))*100)

SRV_score$score_rich_intol <-
  case_when(
    SRV_score$rich_intol > 8.0 ~ 100,
    SRV_score$rich_intol < 4.15 ~ 0,
    TRUE ~ ((SRV_score$rich_intol-4.15)/(8.0-4.15))*100)

SRV_score$score_rich_plecoptera <-
  case_when(
    SRV_score$rich_plecoptera > 3.0 ~ 100,
    SRV_score$rich_plecoptera < 0.65 ~ 0,
    TRUE ~ ((SRV_score$rich_plecoptera-0.65)/(3.0-0.65))*100)

SRV_score$score_rich_predator <-
  case_when(
    SRV_score$rich_predator > 3.0 ~ 100,
    SRV_score$rich_predator < 2.72 ~ 0,
    TRUE ~ ((SRV_score$rich_predator-2.72)/(3.0-2.72))*100)

SRV_score$score_rich_trichoptera <-
  case_when(
    SRV_score$rich_trichoptera > 3.0 ~ 100,
    SRV_score$rich_trichoptera < 2.79 ~ 0,
    TRUE ~ ((SRV_score$rich_trichoptera-2.79)/(3.0-2.79))*100)

# SGV Metric Scoring
SGV_score$score_aspt_mod <-
  case_when(
    SGV_score$aspt_mod > 6.44 ~ 100,
    SGV_score$aspt_mod < 3.94 ~ 0,
    TRUE ~ ((SGV_score$aspt_mod-3.94)/(6.44-3.94))*100)

SGV_score$score_gold <-
  case_when(
    SGV_score$gold > 80.89 ~ 100,
    SGV_score$gold < 64.95 ~ 0,
    TRUE ~ ((SGV_score$gold-64.95)/(80.89-64.95))*100)

SGV_score$score_pct_cling <-
  case_when(
    SGV_score$pct_cling > 62.23 ~ 100,
    SGV_score$pct_cling < 39.7 ~ 0,
    TRUE ~ ((SGV_score$pct_cling-39.7)/(62.23-39.7))*100)

SGV_score$score_pct_collect <-
  case_when(
    SGV_score$pct_collect < 71.17 ~ 100,
    SGV_score$pct_collect > 94.57 ~ 0,
    TRUE ~ ((94.57-SGV_score$pct_collect)/(94.57-71.17))*100)

SGV_score$score_pct_ephemeroptera_no_baetid <-
  case_when(
    SGV_score$pct_ephemeroptera_no_baetid > 23.94 ~ 100,
    SGV_score$pct_ephemeroptera_no_baetid < 5.23 ~ 0,
    TRUE ~ ((SGV_score$pct_ephemeroptera_no_baetid-5.23)/(23.94-5.23))*100)

SGV_score$score_pct_ept_no_hydro <-
  case_when(
    SGV_score$pct_ept_no_hydro > 83.21 ~ 100,
    SGV_score$pct_ept_no_hydro < 73.53 ~ 0,
    TRUE ~ ((SGV_score$pct_ept_no_hydro-73.53)/(83.21-73.53))*100)

SGV_score$score_pct_ept_rich_no_tol <-
  case_when(
    SGV_score$pct_ept_rich_no_tol > 32.97 ~ 100,
    SGV_score$pct_ept_rich_no_tol < 19.39 ~ 0,
    TRUE ~ ((SGV_score$pct_ept_rich_no_tol-19.39)/(32.97-19.39))*100)

SGV_score$score_pct_filter <-
  case_when(
    SGV_score$pct_filter > 22.79 ~ 100,
    SGV_score$pct_filter < 8.6 ~ 0,
    TRUE ~ ((SGV_score$pct_filter-8.6)/(22.79-8.6))*100)

SGV_score$score_pct_heptageniidae <-
  case_when(
    SGV_score$pct_heptageniidae > 7.89 ~ 100,
    SGV_score$pct_heptageniidae < 0 ~ 0,
    TRUE ~ ((SGV_score$pct_heptageniidae-0)/(7.89-0))*100)

SGV_score$score_pct_mod_tol_4_6 <-
  case_when(
    SGV_score$pct_mod_tol_4_6 < 52.49 ~ 100,
    SGV_score$pct_mod_tol_4_6 > 76.71 ~ 0,
    TRUE ~ ((76.71-SGV_score$pct_mod_tol_4_6)/(76.71-52.49))*100)

SGV_score$score_pct_non_hydrop_trichoptera <-
  case_when(
    SGV_score$pct_non_hydrop_trichoptera > 30.54 ~ 100,
    SGV_score$pct_non_hydrop_trichoptera < 10.57 ~ 0,
    TRUE ~ ((SGV_score$pct_non_hydrop_trichoptera-10.57)/(30.54-10.57))*100)

SGV_score$score_pct_pisciforma <-
  case_when(
    SGV_score$pct_pisciforma > 17.07 ~ 100,
    SGV_score$pct_pisciforma < 0.1 ~ 0,
    TRUE ~ ((SGV_score$pct_pisciforma-0.1)/(17.07-0.1))*100)

SGV_score$score_pct_predator <-
  case_when(
    SGV_score$pct_predator > 5.28 ~ 100,
    SGV_score$pct_predator < 0.02 ~ 0,
    TRUE ~ ((SGV_score$pct_predator-0.02)/(5.28-0.02))*100)

SGV_score$score_pct_pterygota <-
  case_when(
    SGV_score$pct_pterygota > 96.88 ~ 100,
    SGV_score$pct_pterygota < 90.2 ~ 0,
    TRUE ~ ((SGV_score$pct_pterygota-90.2)/(96.88-90.2))*100)

SGV_score$score_pct_scrape <-
  case_when(
    SGV_score$pct_scrape > 15.43 ~ 100,
    SGV_score$pct_scrape < 0.32 ~ 0,
    TRUE ~ ((SGV_score$pct_scrape-0.32)/(15.43-0.32))*100)

SGV_score$score_pct_sprawl <-
  case_when(
    SGV_score$pct_sprawl < 1.79 ~ 100,
    SGV_score$pct_sprawl > 4.54 ~ 0,
    TRUE ~ ((4.54-SGV_score$pct_sprawl)/(4.54-1.79))*100)

SGV_score$score_pct_tolerant_7_10 <-
  case_when(
    SGV_score$pct_tolerant_7_10 < 1.73 ~ 100,
    SGV_score$pct_tolerant_7_10 > 4.3 ~ 0,
    TRUE ~ ((4.3-SGV_score$pct_tolerant_7_10)/(4.3-1.73))*100)

SGV_score$score_rich_collect <-
  case_when(
    SGV_score$rich_collect > 8.0 ~ 100,
    SGV_score$rich_collect < 6.32 ~ 0,
    TRUE ~ ((SGV_score$rich_collect-6.32)/(8.0-6.32))*100)

SGV_score$score_rich_ephemeroptera <-
  case_when(
    SGV_score$rich_ephemeroptera > 3.94 ~ 100,
    SGV_score$rich_ephemeroptera < 0.61 ~ 0,
    TRUE ~ ((SGV_score$rich_ephemeroptera-0.61)/(3.94-0.61))*100)

SGV_score$score_rich_filter <-
  case_when(
    SGV_score$rich_filter > 3.0 ~ 100,
    SGV_score$rich_filter < 1.36 ~ 0,
    TRUE ~ ((SGV_score$rich_filter-1.36)/(3.0-1.36))*100)

SGV_score$score_rich_trichoptera <-
  case_when(
    SGV_score$rich_trichoptera > 2.54 ~ 100,
    SGV_score$rich_trichoptera < 0.08 ~ 0,
    TRUE ~ ((SGV_score$rich_trichoptera-0.08)/(2.54-0.08))*100)

# LNP Metric Scoring
LNP_score$score_aspt_mod <-
  case_when(
    LNP_score$aspt_mod > 6.62 ~ 100,
    LNP_score$aspt_mod < 2.7 ~ 0,
    TRUE ~ ((LNP_score$aspt_mod-2.7)/(6.62-2.7))*100)

LNP_score$score_pct_cling <-
  case_when(
    LNP_score$pct_cling > 62.77 ~ 100,
    LNP_score$pct_cling < 4.54 ~ 0,
    TRUE ~ ((LNP_score$pct_cling-4.54)/(62.77-4.54))*100)

LNP_score$score_pct_collect <-
  case_when(
    LNP_score$pct_collect < 72.82 ~ 100,
    LNP_score$pct_collect > 92.49 ~ 0,
    TRUE ~ ((92.49-LNP_score$pct_collect)/(92.49-72.82))*100)

LNP_score$score_pct_ept_rich <-
  case_when(
    LNP_score$pct_ept_rich > 52.32 ~ 100,
    LNP_score$pct_ept_rich < 38.07 ~ 0,
    TRUE ~ ((LNP_score$pct_ept_rich-38.07)/(52.32-38.07))*100)

LNP_score$score_pct_mod_tol_4_6 <-
  case_when(
    LNP_score$pct_mod_tol_4_6 < 46.54 ~ 100,
    LNP_score$pct_mod_tol_4_6 > 77.78 ~ 0,
    TRUE ~ ((77.78-LNP_score$pct_mod_tol_4_6)/(77.78-46.54))*100)

LNP_score$score_pct_predator <-
  case_when(
    LNP_score$pct_predator > 7.78 ~ 100,
    LNP_score$pct_predator < 1.43 ~ 0,
    TRUE ~ ((LNP_score$pct_predator-1.43)/(7.78-1.43))*100)

LNP_score$score_pct_pterygota <-
  case_when(
    LNP_score$pct_pterygota > 98.36 ~ 100,
    LNP_score$pct_pterygota < 92.09 ~ 0,
    TRUE ~ ((LNP_score$pct_pterygota-92.09)/(98.36-92.09))*100)

LNP_score$score_rich_cling <-
  case_when(
    LNP_score$rich_cling > 7.0 ~ 100,
    LNP_score$rich_cling < 3.69 ~ 0,
    TRUE ~ ((LNP_score$rich_cling-3.69)/(7.0-3.69))*100)

LNP_score$score_rich_ephemeroptera <-
  case_when(
    LNP_score$rich_ephemeroptera > 3.0 ~ 100,
    LNP_score$rich_ephemeroptera < 1.7 ~ 0,
    TRUE ~ ((LNP_score$rich_ephemeroptera-1.7)/(3.0-1.7))*100)

# PIED Metric Scoring
PIED_score$score_pct_cling <-
  case_when(
    PIED_score$pct_cling > 62.62 ~ 100,
    PIED_score$pct_cling < 32.33 ~ 0,
    TRUE ~ ((PIED_score$pct_cling-32.33)/(62.62-32.33))*100)

PIED_score$score_pct_ept_rich <-
  case_when(
    PIED_score$pct_ept_rich > 53.58 ~ 100,
    PIED_score$pct_ept_rich < 38.51 ~ 0,
    TRUE ~ ((PIED_score$pct_ept_rich-38.51)/(53.58-38.51))*100)

PIED_score$score_pct_intol_0_4 <-
  case_when(
    PIED_score$pct_intol_0_4 > 58.15 ~ 100,
    PIED_score$pct_intol_0_4 < 27.97 ~ 0,
    TRUE ~ ((PIED_score$pct_intol_0_4-27.97)/(58.15-27.97))*100)

PIED_score$score_pct_tolerant_7_10 <- 
  case_when(
    PIED_score$pct_tolerant_7_10 < 1.6 ~ 100,
    PIED_score$pct_tolerant_7_10 > 4.84 ~ 0,
    TRUE ~ ((4.84-PIED_score$pct_tolerant_7_10)/(4.84-1.6))*100)

PIED_score$score_pct_urban_intol <-
  case_when(
    PIED_score$pct_urban_intol > 95.41 ~ 100,
    PIED_score$pct_urban_intol < 86.93 ~ 0,
    TRUE ~ ((PIED_score$pct_urban_intol-86.93)/(95.41-86.93))*100)

PIED_score$score_rich_cling <-
  case_when(
    PIED_score$rich_cling > 8.0 ~ 100,
    PIED_score$rich_cling < 2.96 ~ 0,
    TRUE ~ ((PIED_score$rich_cling-2.96)/(8.0-2.96))*100)
  
PIED_score$score_rich_ephemeroptera <-  
  case_when(
    PIED_score$rich_ephemeroptera > 3.6 ~ 100,
    PIED_score$rich_ephemeroptera < 0.86 ~ 0,
    TRUE ~ ((PIED_score$rich_ephemeroptera-0.86)/(3.6-0.86))*100)

PIED_score$score_rich_trichoptera <-
  case_when(
    PIED_score$rich_trichoptera > 2.03 ~ 100,
    PIED_score$rich_trichoptera < 0.45 ~ 0,
    TRUE ~ ((PIED_score$rich_trichoptera-0.45)/(2.03-0.45))*100)
      
# Average scores to get the Chessie BIBI
NAPU_score$chessie_bibi <- rowMeans(NAPU_score[ ,25:43], na.rm = T)
NCA_score$chessie_bibi  <- rowMeans(NCA_score[ ,15:23], na.rm = T)
CA_score$chessie_bibi   <- rowMeans(CA_score[ ,19:31], na.rm = T)
MAC_score$chessie_bibi  <- rowMeans(MAC_score[ ,18:29], na.rm = T)
SEP_score$chessie_bibi  <- rowMeans(SEP_score[ ,17:26], na.rm = T)
BLUE_score$chessie_bibi <- rowMeans(BLUE_score[ ,11:15], na.rm = T) 
NRV_score$chessie_bibi  <- rowMeans(NRV_score[ ,14:21], na.rm = T)
SRV_score$chessie_bibi  <- rowMeans(SRV_score[ ,21:35], na.rm = T)
UNP_score$chessie_bibi  <- rowMeans(UNP_score[ ,18:29], na.rm = T)
SGV_score$chessie_bibi  <- rowMeans(SGV_score[ ,27:47], na.rm = T) 
LNP_score$chessie_bibi  <- rowMeans(LNP_score[ ,15:23], na.rm = T)
PIED_score$chessie_bibi <- rowMeans(PIED_score[ ,14:21], na.rm = T) 
      
# Assign narrative ratings
NAPU_score$rating <-
  case_when(
    NAPU_score$chessie_bibi < 17.93 ~ "very_poor",
    NAPU_score$chessie_bibi >= 17.93 & NAPU_score$chessie_bibi < 35.85 ~ "poor",
    NAPU_score$chessie_bibi >= 35.85 & NAPU_score$chessie_bibi < 47.42 ~ "fair",
    NAPU_score$chessie_bibi >= 47.42 & NAPU_score$chessie_bibi < 61.95 ~ "good",
    NAPU_score$chessie_bibi >= 61.95 ~ "excellent",
    TRUE ~ "ERROR")

NCA_score$rating <-
  case_when(
    NCA_score$chessie_bibi < 15.73 ~ "very_poor",
    NCA_score$chessie_bibi >= 15.73 & NCA_score$chessie_bibi < 31.46 ~ "poor",
    NCA_score$chessie_bibi >= 31.46 & NCA_score$chessie_bibi < 56.25 ~ "fair",
    NCA_score$chessie_bibi >= 56.25 & NCA_score$chessie_bibi < 78.66 ~ "good",
    NCA_score$chessie_bibi >= 78.66 ~ "excellent",
    TRUE ~ "ERROR")
  
CA_score$rating <-
  case_when(
    CA_score$chessie_bibi < 27.88 ~ "very_poor",
    CA_score$chessie_bibi >= 27.88 & CA_score$chessie_bibi < 55.77 ~ "poor",
    CA_score$chessie_bibi >= 55.77 & CA_score$chessie_bibi < 69.55 ~ "fair",
    CA_score$chessie_bibi >= 69.55 & CA_score$chessie_bibi < 78.56 ~ "good",
    CA_score$chessie_bibi >= 78.56 ~ "excellent",
    TRUE ~ "ERROR")

MAC_score$rating <-
  case_when(
    MAC_score$chessie_bibi < 22.72 ~ "very_poor",
    MAC_score$chessie_bibi >= 22.72 & MAC_score$chessie_bibi < 45.43 ~ "poor",
    MAC_score$chessie_bibi >= 45.43 & MAC_score$chessie_bibi < 63.18 ~ "fair",
    MAC_score$chessie_bibi >= 63.18 & MAC_score$chessie_bibi < 76.85 ~ "good",
    MAC_score$chessie_bibi >= 76.85 ~ "excellent",
    TRUE ~ "ERROR")
  
SEP_score$rating <-
  case_when(
    SEP_score$chessie_bibi < 16.28 ~ "very_poor",
    SEP_score$chessie_bibi >= 16.28 & SEP_score$chessie_bibi < 32.55 ~ "poor",
    SEP_score$chessie_bibi >= 32.55 & SEP_score$chessie_bibi < 56.66 ~ "fair",
    SEP_score$chessie_bibi >= 56.66 & SEP_score$chessie_bibi < 83.90 ~ "good",
    SEP_score$chessie_bibi >= 83.90 ~ "excellent",
    TRUE ~ "ERROR")
  
BLUE_score$rating <-
  case_when(
    BLUE_score$chessie_bibi < 30.87 ~ "very_poor",
    BLUE_score$chessie_bibi >= 30.87 & BLUE_score$chessie_bibi < 61.74 ~ "poor",
    BLUE_score$chessie_bibi >= 61.74 & BLUE_score$chessie_bibi < 82.21 ~ "fair",
    BLUE_score$chessie_bibi >= 82.21 & BLUE_score$chessie_bibi < 92.07 ~ "good",
    BLUE_score$chessie_bibi >= 92.07 ~ "excellent",
    TRUE ~ "ERROR")
  
NRV_score$rating <-
  case_when(
    NRV_score$chessie_bibi < 19.10 ~ "very_poor",
    NRV_score$chessie_bibi >= 19.10 & NRV_score$chessie_bibi < 38.21 ~ "poor",
    NRV_score$chessie_bibi >= 38.21 & NRV_score$chessie_bibi < 50.81 ~ "fair",
    NRV_score$chessie_bibi >= 50.81 & NRV_score$chessie_bibi < 75.00 ~ "good",
    NRV_score$chessie_bibi >= 75.00 ~ "excellent",
    TRUE ~ "ERROR")
  
SRV_score$rating <-
  case_when(
    SRV_score$chessie_bibi < 21.58 ~ "very_poor",
    SRV_score$chessie_bibi >= 21.58 & SRV_score$chessie_bibi < 43.16 ~ "poor",
    SRV_score$chessie_bibi >= 43.16 & SRV_score$chessie_bibi < 57.97 ~ "fair",
    SRV_score$chessie_bibi >= 57.97 & SRV_score$chessie_bibi < 71.93 ~ "good",
    SRV_score$chessie_bibi >= 71.93 ~ "excellent",
    TRUE ~ "ERROR")

UNP_score$rating <-
  case_when(
    UNP_score$chessie_bibi < 31.30 ~ "very_poor",
    UNP_score$chessie_bibi >= 31.30 & UNP_score$chessie_bibi < 62.59 ~ "poor",
    UNP_score$chessie_bibi >= 62.59 & UNP_score$chessie_bibi < 69.85 ~ "fair",
    UNP_score$chessie_bibi >= 69.85 & UNP_score$chessie_bibi < 80.52 ~ "good",
    UNP_score$chessie_bibi >= 80.52 ~ "excellent",
    TRUE ~ "ERROR")

SGV_score$rating <-
  case_when(
    SGV_score$chessie_bibi < 29.93 ~ "very_poor",
    SGV_score$chessie_bibi >= 29.93 & SGV_score$chessie_bibi < 59.86 ~ "poor",
    SGV_score$chessie_bibi >= 59.86 & SGV_score$chessie_bibi < 66.02 ~ "fair",
    SGV_score$chessie_bibi >= 66.02 & SGV_score$chessie_bibi < 76.48 ~ "good",
    SGV_score$chessie_bibi >= 76.48 ~ "excellent",
    TRUE ~ "ERROR")

LNP_score$rating <-
  case_when(
    LNP_score$chessie_bibi < 35.16 ~ "very_poor",
    LNP_score$chessie_bibi >= 35.16 & LNP_score$chessie_bibi < 70.32 ~ "poor",
    LNP_score$chessie_bibi >= 70.32 & LNP_score$chessie_bibi < 80.15 ~ "fair",
    LNP_score$chessie_bibi >= 80.15 & LNP_score$chessie_bibi < 91.30 ~ "good",
    LNP_score$chessie_bibi >= 91.30 ~ "excellent",
    TRUE ~ "ERROR")

PIED_score$rating <-
  case_when(
    PIED_score$chessie_bibi < 30.31 ~ "very_poor",
    PIED_score$chessie_bibi >= 30.31 & PIED_score$chessie_bibi < 60.62 ~ "poor",
    PIED_score$chessie_bibi >= 60.62 & PIED_score$chessie_bibi < 73.18 ~ "fair",
    PIED_score$chessie_bibi >= 73.18 & PIED_score$chessie_bibi < 81.83 ~ "good",
    PIED_score$chessie_bibi >= 81.83 ~ "excellent",
    TRUE ~ "ERROR")
      
# Merge dataframes
NAPU_sub <- NAPU_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
NCA_sub <- NCA_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
CA_sub <- CA_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
MAC_sub <- MAC_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
SEP_sub <- SEP_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
BLUE_sub <- BLUE_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
NRV_sub <- NRV_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
SRV_sub <- SRV_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
UNP_sub <- UNP_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
SGV_sub <- SGV_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
LNP_sub <- LNP_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)
PIED_sub <- PIED_score %>% select(unique_id, icprb_bioregion_id, station_id, sample_date, sample_time, chessie_bibi, rating)

bind1 <- rbind(NAPU_sub, NCA_sub)
bind2 <- rbind(bind1, CA_sub)
bind3 <- rbind(bind2, MAC_sub)
bind4 <- rbind(bind3, SEP_sub)
bind5 <- rbind(bind4, BLUE_sub)
bind6 <- rbind(bind5, NRV_sub)
bind7 <- rbind(bind6, SRV_sub)
bind8 <- rbind(bind7, UNP_sub)
bind9 <- rbind(bind8, SGV_sub)
bind10 <- rbind(bind9, LNP_sub)
bibi_ratings_combined <- rbind(bind10, PIED_sub)

bibi_ratings_combined <- bibi_ratings_combined %>%
  separate(unique_id, c('event_id', 'sample_num'), sep="_", remove = FALSE)

final_score <- bibi_ratings_combined %>%
  select(unique_id, chessie_bibi)
      
 return(final_score)
 }
```
**2.Program Repetition**
Now that the Chessie BIBI program has been incorporated into a function, it will be repeated a set number of times. The program is currently set to repeat 100 times. If you would prefer to repeat the program either greater than or fewer than 100 times, simple change `100` to your own number in `repeat_me <- 100`.This process will take several hours, so you may want to let the program run overnight.

The dataframe, **sc.df**, contains the reteated Chessie BIBI scores.
```{r}
non_rare <- taxa.df 

repeat.scores.df <- non_rare %>%              
 select(unique_id) %>%
 unique() %>%
 arrange(unique_id)

sc.df <- repeat.scores.df

repeat_me <- 100 

for(i in 1:repeat_me) {

  final_score <- ChessieBIBI_score()
   
  sc.df <-  full_join(sc.df, final_score, by = "unique_id")
  
  colnames(sc.df)[colnames(sc.df)=="chessie_bibi"] <- as.character(paste0("chessie_bibi", i))
}
```
**3. Prep the dataframes and assign narrative ratings**
Add the bioregion assignments to the repeated scores dataframe. Then, the Chessie BIBI scores in will be assigned a narrative rating.

If you have changed the number of repetitions to another number than 100 in the step above, you will need to change this line of code, near the end of the this code chunk to reflect your number: `names(ratings.df2)[3:102] <- paste0(c("rating"),1:100)`

The dataframe produced at the end of this step, **join.df**, contains the repeated scores and ratings.
```{r}
bioregions <- non_rare %>%
  select(icprb_bioregion_id, unique_id) %>%
  unique()

ratings.df <- sc.df

ratings.df2 <- full_join(ratings.df, bioregions, by = "unique_id")
ratings.df2 <- ratings.df2 %>%  
  select(icprb_bioregion_id, everything())

ratings.df2 <- ratings.df2 %>%
  mutate_at(vars(starts_with("chessie_")), list( 
  ~case_when(
    ratings.df2$icprb_bioregion_id == "napu" & . < 17.93 ~ "very_poor",
    ratings.df2$icprb_bioregion_id == "napu" & (. >= 17.93 & .< 35.85) ~ "poor",
    ratings.df2$icprb_bioregion_id == "napu" & (. >= 35.85 & .< 47.42) ~ "fair",
    ratings.df2$icprb_bioregion_id == "napu" & (. >= 47.42 & .< 61.95) ~ "good",
    ratings.df2$icprb_bioregion_id == "napu" & . >= 61.95 ~ "excellent",
    
    ratings.df2$icprb_bioregion_id == "nca" & . < 15.73 ~ "very_poor",
    ratings.df2$icprb_bioregion_id == "nca" & (. >= 15.73 & .< 31.46) ~ "poor",
    ratings.df2$icprb_bioregion_id == "nca" & (. >= 31.46 & .< 56.25) ~ "fair",
    ratings.df2$icprb_bioregion_id == "nca" & (. >= 56.25 & .< 78.66) ~ "good",
    ratings.df2$icprb_bioregion_id == "nca" & . >= 78.66 ~ "excellent",
    
    ratings.df2$icprb_bioregion_id == "ca" & . < 27.88 ~ "very_poor",
    ratings.df2$icprb_bioregion_id == "ca" & (. >= 27.88 & .< 55.77) ~ "poor",
    ratings.df2$icprb_bioregion_id == "ca" & (. >= 55.77 & .< 69.55) ~ "fair",
    ratings.df2$icprb_bioregion_id == "ca" & (. >= 69.55 & .< 78.56)~ "good",
    ratings.df2$icprb_bioregion_id == "ca" & . >= 78.56 ~ "excellent",
    
    ratings.df2$icprb_bioregion_id == "mac" & . < 22.72 ~ "very_poor",
    ratings.df2$icprb_bioregion_id == "mac" & (. >= 22.72 & .< 45.43) ~ "poor",
    ratings.df2$icprb_bioregion_id == "mac" & (. >= 45.43 & .< 63.18) ~ "fair",
    ratings.df2$icprb_bioregion_id == "mac" & (. >= 63.18 & .< 76.85) ~ "good",
    ratings.df2$icprb_bioregion_id == "mac" & . >= 76.85 ~ "excellent",
    
    ratings.df2$icprb_bioregion_id == "sep" & . < 16.28 ~ "very_poor",
    ratings.df2$icprb_bioregion_id == "sep" & (. >= 16.28 & .< 32.55) ~ "poor",
    ratings.df2$icprb_bioregion_id == "sep" & (. >= 32.55 & .< 56.66) ~ "fair",
    ratings.df2$icprb_bioregion_id == "sep" & (. >= 56.66 & .< 83.90) ~ "good",
    ratings.df2$icprb_bioregion_id == "sep" & . >= 83.90 ~ "excellent",
    
    ratings.df2$icprb_bioregion_id == "blue" & . < 30.87 ~ "very_poor",
    ratings.df2$icprb_bioregion_id == "blue" & (. >= 30.87 & .< 61.74) ~ "poor",
    ratings.df2$icprb_bioregion_id == "blue" & (. >= 61.74 & .< 82.21) ~ "fair",
    ratings.df2$icprb_bioregion_id == "blue" & (. >= 82.21 & .< 92.07) ~ "good",
    ratings.df2$icprb_bioregion_id == "blue" & . >= 92.07 ~ "excellent",
    
    ratings.df2$icprb_bioregion_id == "nrv" & . < 19.10 ~ "very_poor",
    ratings.df2$icprb_bioregion_id == "nrv" & (. >= 19.10 & .< 38.21) ~ "poor",
    ratings.df2$icprb_bioregion_id == "nrv" & (. >= 38.21 & .< 50.81)~ "fair",
    ratings.df2$icprb_bioregion_id == "nrv" & (. >= 50.81 & .< 75.00) ~ "good",
    ratings.df2$icprb_bioregion_id == "nrv" & . >= 75.00 ~ "excellent",
    
    ratings.df2$icprb_bioregion_id == "srv" & . < 21.58 ~ "very_poor",
    ratings.df2$icprb_bioregion_id == "srv" & (. >= 21.58 & .< 43.16) ~ "poor",
    ratings.df2$icprb_bioregion_id == "srv" & (. >= 43.16 & .< 57.97) ~ "fair",
    ratings.df2$icprb_bioregion_id == "srv" & (. >= 57.97 & .< 71.93) ~ "good",
    ratings.df2$icprb_bioregion_id == "srv" & . >= 71.93 ~ "excellent",
    
    ratings.df2$icprb_bioregion_id == "unp" & . < 31.30 ~ "very_poor",
    ratings.df2$icprb_bioregion_id == "unp" & (. >= 31.30 & .< 62.59) ~ "poor",
    ratings.df2$icprb_bioregion_id == "unp" & (. >= 62.59 & .< 69.85) ~ "fair",
    ratings.df2$icprb_bioregion_id == "unp" & (. >= 69.85 & .< 80.52) ~ "good",
    ratings.df2$icprb_bioregion_id == "unp" & . >= 80.52 ~ "excellent",
    
    ratings.df2$icprb_bioregion_id == "sgv" & . < 29.93 ~ "very_poor",
    ratings.df2$icprb_bioregion_id == "sgv" & (. >= 29.93 & .< 59.86) ~ "poor",
    ratings.df2$icprb_bioregion_id == "sgv" & (. >= 59.86 & .< 66.02) ~ "fair",
    ratings.df2$icprb_bioregion_id == "sgv" & (. >= 66.02 & .< 76.48) ~ "good",
    ratings.df2$icprb_bioregion_id == "sgv" & . >= 76.48 ~ "excellent",
    
    ratings.df2$icprb_bioregion_id == "lnp" & . < 35.16 ~ "very_poor",
    ratings.df2$icprb_bioregion_id == "lnp" & (. >= 35.16 & .< 70.32) ~ "poor",
    ratings.df2$icprb_bioregion_id == "lnp" & (. >= 70.32 & .< 80.15) ~ "fair",
    ratings.df2$icprb_bioregion_id == "lnp" & (. >= 80.15 & .< 91.30) ~ "good",
    ratings.df2$icprb_bioregion_id == "lnp" & . >= 91.30 ~ "excellent",
    
    ratings.df2$icprb_bioregion_id == "pied" & . < 30.31 ~ "very_poor",
    ratings.df2$icprb_bioregion_id == "pied" & (. >= 30.31 & .< 60.62) ~ "poor",
    ratings.df2$icprb_bioregion_id == "pied" & (. >= 60.62 & .< 73.18) ~ "fair",
    ratings.df2$icprb_bioregion_id == "pied" & (. >= 73.18 & .< 81.83) ~ "good",
    ratings.df2$icprb_bioregion_id == "pied" & . >= 81.83 ~ "excellent",
    TRUE ~ "ERROR")
    ))

names(ratings.df2)[3:102] <- paste0(c("rating"),1:100)

ChessieBIBI_repeated_ratings <- ratings.df2
ChessieBIBI_repeated_scores <- sc.df

join.df <- full_join(ChessieBIBI_repeated_scores, ChessieBIBI_repeated_ratings, by="unique_id")

join.df <- join.df %>%
  select(icprb_bioregion_id, everything())
```
**4. Analysis**
Determine the mean and median ratings of the program repetitions. This step will also calculate how often the mean or median rating occurs out of the 100 runs, as well as a comparison of whether or not the median and mean has the same rating.

As before, if you have changed the number of repetitions to another number than 100 in the step above, you will need to change the 2nd - 7th lines of code below to reflect your number of repetitions. The column numbers for lines 2-6 should be the column numbers that contain the narrative ratings. Line 7 should be the first to the last columns containing the scores.
```{r}
stats.df <- join.df
stats.df$very_poor <- rowSums(stats.df[ ,103:202] == "very_poor")
stats.df$poor = rowSums(stats.df[,103:202] == "poor")
stats.df$fair = rowSums(stats.df[,103:202] == "fair")
stats.df$good = rowSums(stats.df[,103:202] == "good")
stats.df$excellent = rowSums(stats.df[,103:202] == "excellent")
stats.df$median_score = apply(stats.df[ ,3:102],1, median)
stats.df$median_rating <- 
   case_when(
    stats.df$icprb_bioregion_id == "napu" & stats.df$median_score < 17.93 ~ "very_poor",
    stats.df$icprb_bioregion_id == "napu" & (stats.df$median_score >= 17.93 & stats.df$median_score < 35.85) ~ "poor",
    stats.df$icprb_bioregion_id == "napu" & (stats.df$median_score >= 35.85 & stats.df$median_score < 47.42) ~ "fair",
    stats.df$icprb_bioregion_id == "napu" & (stats.df$median_score >= 47.42 & stats.df$median_score < 61.95) ~ "good",
    stats.df$icprb_bioregion_id == "napu" & stats.df$median_score >= 61.95 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "nca" & stats.df$median_score < 15.73 ~ "very_poor",
    stats.df$icprb_bioregion_id == "nca" & (stats.df$median_score >= 15.73 & stats.df$median_score < 31.46) ~ "poor",
    stats.df$icprb_bioregion_id == "nca" & (stats.df$median_score >= 31.46 & stats.df$median_score < 56.25) ~ "fair",
    stats.df$icprb_bioregion_id == "nca" & (stats.df$median_score >= 56.25 & stats.df$median_score < 78.66) ~ "good",
    stats.df$icprb_bioregion_id == "nca" & stats.df$median_score >= 78.66 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "ca" & stats.df$median_score < 27.88 ~ "very_poor",
    stats.df$icprb_bioregion_id == "ca" & (stats.df$median_score >= 27.88 & stats.df$median_score < 55.77) ~ "poor",
    stats.df$icprb_bioregion_id == "ca" & (stats.df$median_score>= 55.77 & stats.df$median_score < 69.55) ~ "fair",
    stats.df$icprb_bioregion_id == "ca" & (stats.df$median_score>= 69.55 & stats.df$median_score < 78.56)~ "good",
    stats.df$icprb_bioregion_id == "ca" & stats.df$median_score>= 78.56 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "mac" & stats.df$median_score < 22.72 ~ "very_poor",
    stats.df$icprb_bioregion_id == "mac" & (stats.df$median_score>= 22.72 & stats.df$median_score < 45.43) ~ "poor",
    stats.df$icprb_bioregion_id == "mac" & (stats.df$median_score>= 45.43 & stats.df$median_score < 63.18) ~ "fair",
    stats.df$icprb_bioregion_id == "mac" & (stats.df$median_score>= 63.18 & stats.df$median_score < 76.85) ~ "good",
    stats.df$icprb_bioregion_id == "mac" & stats.df$median_score>= 76.85 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "sep" & stats.df$median_score < 16.28 ~ "very_poor",
    stats.df$icprb_bioregion_id == "sep" & (stats.df$median_score>= 16.28 & stats.df$median_score < 32.55) ~ "poor",
    stats.df$icprb_bioregion_id == "sep" & (stats.df$median_score>= 32.55 & stats.df$median_score < 56.66) ~ "fair",
    stats.df$icprb_bioregion_id == "sep" & (stats.df$median_score>= 56.66 & stats.df$median_score < 83.90) ~ "good",
    stats.df$icprb_bioregion_id == "sep" & stats.df$median_score>= 83.90 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "blue" & stats.df$median_score < 30.87 ~ "very_poor",
    stats.df$icprb_bioregion_id == "blue" & (stats.df$median_score >= 30.87 & stats.df$median_score < 61.74) ~ "poor",
    stats.df$icprb_bioregion_id == "blue" & (stats.df$median_score >= 61.74 & stats.df$median_score< 82.21) ~ "fair",
    stats.df$icprb_bioregion_id == "blue" & (stats.df$median_score >= 82.21 & stats.df$median_score < 92.07) ~ "good",
    stats.df$icprb_bioregion_id == "blue" & stats.df$median_score >= 92.07 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "nrv" & stats.df$median_score < 19.10 ~ "very_poor",
    stats.df$icprb_bioregion_id == "nrv" & (stats.df$median_score >= 19.10 & stats.df$median_score < 38.21) ~ "poor",
    stats.df$icprb_bioregion_id == "nrv" & (stats.df$median_score >= 38.21 & stats.df$median_score < 50.81)~ "fair",
    stats.df$icprb_bioregion_id == "nrv" & (stats.df$median_score >= 50.81 & stats.df$median_score < 75.00) ~ "good",
    stats.df$icprb_bioregion_id == "nrv" & stats.df$median_score >= 75.00 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "srv" & stats.df$median_score < 21.58 ~ "very_poor",
    stats.df$icprb_bioregion_id == "srv" & (stats.df$median_score>= 21.58 & stats.df$median_score < 43.16) ~ "poor",
    stats.df$icprb_bioregion_id == "srv" & (stats.df$median_score>= 43.16 & stats.df$median_score < 57.97) ~ "fair",
    stats.df$icprb_bioregion_id == "srv" & (stats.df$median_score>= 57.97 & stats.df$median_score < 71.93) ~ "good",
    stats.df$icprb_bioregion_id == "srv" & stats.df$median_score>= 71.93 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "unp" & stats.df$median_score < 31.30 ~ "very_poor",
    stats.df$icprb_bioregion_id == "unp" & (stats.df$median_score >= 31.30 & stats.df$median_score < 62.59) ~ "poor",
    stats.df$icprb_bioregion_id == "unp" & (stats.df$median_score >= 62.59 & stats.df$median_score < 69.85) ~ "fair",
    stats.df$icprb_bioregion_id == "unp" & (stats.df$median_score >= 69.85 & stats.df$median_score < 80.52) ~ "good",
    stats.df$icprb_bioregion_id == "unp" & stats.df$median_score >= 80.52 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "sgv" & stats.df$median_score < 29.93 ~ "very_poor",
    stats.df$icprb_bioregion_id == "sgv" & (stats.df$median_score >= 29.93 & stats.df$median_score < 59.86) ~ "poor",
    stats.df$icprb_bioregion_id == "sgv" & (stats.df$median_score >= 59.86 & stats.df$median_score < 66.02) ~ "fair",
    stats.df$icprb_bioregion_id == "sgv" & (stats.df$median_score >= 66.02 & stats.df$median_score < 76.48) ~ "good",
    stats.df$icprb_bioregion_id == "sgv" & stats.df$median_score >= 76.48 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "lnp" & stats.df$median_score < 35.16 ~ "very_poor",
    stats.df$icprb_bioregion_id == "lnp" & (stats.df$median_score >= 35.16 & stats.df$median_score < 70.32) ~ "poor",
    stats.df$icprb_bioregion_id == "lnp" & (stats.df$median_score >= 70.32 & stats.df$median_score < 80.15) ~ "fair",
    stats.df$icprb_bioregion_id == "lnp" & (stats.df$median_score >= 80.15 & stats.df$median_score < 91.30) ~ "good",
    stats.df$icprb_bioregion_id == "lnp" & stats.df$median_score >= 91.30 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "pied" & stats.df$median_score < 30.31 ~ "very_poor",
    stats.df$icprb_bioregion_id == "pied" & (stats.df$median_score >= 30.31 & stats.df$median_score < 60.62) ~ "poor",
    stats.df$icprb_bioregion_id == "pied" & (stats.df$median_score >= 60.62 & stats.df$median_score < 73.18) ~ "fair",
    stats.df$icprb_bioregion_id == "pied" & (stats.df$median_score >= 73.18 & stats.df$median_score < 81.83) ~ "good",
    stats.df$icprb_bioregion_id == "pied" & stats.df$median_score >= 81.83 ~ "excellent",
    TRUE ~ "ERROR")
  
##
# As stated above, if you chose to do more or less than 100 repetitions, the numbers in the next 3 lines of code also need to be changed.
##
stats.df$median_percent <- ((rowSums(stats.df[ , 103:202] == stats.df$median_rating))/100)*100
  stats.df$mean_score <- rowMeans((stats.df[ ,3:102]))
  stats.df$st_dev = apply(stats.df[ ,3:102],1, sd)
  stats.df$mean_rating <-
   case_when(
    stats.df$icprb_bioregion_id == "napu" & stats.df$mean_score < 17.93 ~ "very_poor",
    stats.df$icprb_bioregion_id == "napu" & (stats.df$mean_score >= 17.93 & stats.df$mean_score < 35.85) ~ "poor",
    stats.df$icprb_bioregion_id == "napu" & (stats.df$mean_score >= 35.85 & stats.df$mean_score < 47.42) ~ "fair",
    stats.df$icprb_bioregion_id == "napu" & (stats.df$mean_score >= 47.42 & stats.df$mean_score < 61.95) ~ "good",
    stats.df$icprb_bioregion_id == "napu" & stats.df$mean_score >= 61.95 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "nca" & stats.df$mean_score < 15.73 ~ "very_poor",
    stats.df$icprb_bioregion_id == "nca" & (stats.df$mean_score >= 15.73 & stats.df$mean_score < 31.46) ~ "poor",
    stats.df$icprb_bioregion_id == "nca" & (stats.df$mean_score >= 31.46 & stats.df$mean_score < 56.25) ~ "fair",
    stats.df$icprb_bioregion_id == "nca" & (stats.df$mean_score >= 56.25 & stats.df$mean_score < 78.66) ~ "good",
    stats.df$icprb_bioregion_id == "nca" & stats.df$mean_score >= 78.66 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "ca" & stats.df$mean_score < 27.88 ~ "very_poor",
    stats.df$icprb_bioregion_id == "ca" & (stats.df$mean_score >= 27.88 & stats.df$mean_score < 55.77) ~ "poor",
    stats.df$icprb_bioregion_id == "ca" & (stats.df$mean_score>= 55.77 & stats.df$mean_score < 69.55) ~ "fair",
    stats.df$icprb_bioregion_id == "ca" & (stats.df$mean_score>= 69.55 & stats.df$mean_score < 78.56)~ "good",
    stats.df$icprb_bioregion_id == "ca" & stats.df$mean_score>= 78.56 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "mac" & stats.df$mean_score < 22.72 ~ "very_poor",
    stats.df$icprb_bioregion_id == "mac" & (stats.df$mean_score>= 22.72 & stats.df$mean_score < 45.43) ~ "poor",
    stats.df$icprb_bioregion_id == "mac" & (stats.df$mean_score>= 45.43 & stats.df$mean_score < 63.18) ~ "fair",
    stats.df$icprb_bioregion_id == "mac" & (stats.df$mean_score>= 63.18 & stats.df$mean_score < 76.85) ~ "good",
    stats.df$icprb_bioregion_id == "mac" & stats.df$mean_score>= 76.85 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "sep" & stats.df$mean_score < 16.28 ~ "very_poor",
    stats.df$icprb_bioregion_id == "sep" & (stats.df$mean_score>= 16.28 & stats.df$mean_score < 32.55) ~ "poor",
    stats.df$icprb_bioregion_id == "sep" & (stats.df$mean_score>= 32.55 & stats.df$mean_score < 56.66) ~ "fair",
    stats.df$icprb_bioregion_id == "sep" & (stats.df$mean_score>= 56.66 & stats.df$mean_score < 83.90) ~ "good",
    stats.df$icprb_bioregion_id == "sep" & stats.df$mean_score>= 83.90 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "blue" & stats.df$mean_score < 30.87 ~ "very_poor",
    stats.df$icprb_bioregion_id == "blue" & (stats.df$mean_score >= 30.87 & stats.df$mean_score < 61.74) ~ "poor",
    stats.df$icprb_bioregion_id == "blue" & (stats.df$mean_score >= 61.74 & stats.df$mean_score< 82.21) ~ "fair",
    stats.df$icprb_bioregion_id == "blue" & (stats.df$mean_score >= 82.21 & stats.df$mean_score < 92.07) ~ "good",
    stats.df$icprb_bioregion_id == "blue" & stats.df$mean_score >= 92.07 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "nrv" & stats.df$mean_score < 19.10 ~ "very_poor",
    stats.df$icprb_bioregion_id == "nrv" & (stats.df$mean_score >= 19.10 & stats.df$mean_score < 38.21) ~ "poor",
    stats.df$icprb_bioregion_id == "nrv" & (stats.df$mean_score >= 38.21 & stats.df$mean_score < 50.81)~ "fair",
    stats.df$icprb_bioregion_id == "nrv" & (stats.df$mean_score >= 50.81 & stats.df$mean_score < 75.00) ~ "good",
    stats.df$icprb_bioregion_id == "nrv" & stats.df$mean_score >= 75.00 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "srv" & stats.df$mean_score < 21.58 ~ "very_poor",
    stats.df$icprb_bioregion_id == "srv" & (stats.df$mean_score>= 21.58 & stats.df$mean_score < 43.16) ~ "poor",
    stats.df$icprb_bioregion_id == "srv" & (stats.df$mean_score>= 43.16 & stats.df$mean_score < 57.97) ~ "fair",
    stats.df$icprb_bioregion_id == "srv" & (stats.df$mean_score>= 57.97 & stats.df$mean_score < 71.93) ~ "good",
    stats.df$icprb_bioregion_id == "srv" & stats.df$mean_score>= 71.93 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "unp" & stats.df$mean_score < 31.30 ~ "very_poor",
    stats.df$icprb_bioregion_id == "unp" & (stats.df$mean_score >= 31.30 & stats.df$mean_score < 62.59) ~ "poor",
    stats.df$icprb_bioregion_id == "unp" & (stats.df$mean_score >= 62.59 & stats.df$mean_score < 69.85) ~ "fair",
    stats.df$icprb_bioregion_id == "unp" & (stats.df$mean_score >= 69.85 & stats.df$mean_score < 80.52) ~ "good",
    stats.df$icprb_bioregion_id == "unp" & stats.df$mean_score >= 80.52 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "sgv" & stats.df$mean_score < 29.93 ~ "very_poor",
    stats.df$icprb_bioregion_id == "sgv" & (stats.df$mean_score >= 29.93 & stats.df$mean_score < 59.86) ~ "poor",
    stats.df$icprb_bioregion_id == "sgv" & (stats.df$mean_score >= 59.86 & stats.df$mean_score < 66.02) ~ "fair",
    stats.df$icprb_bioregion_id == "sgv" & (stats.df$mean_score >= 66.02 & stats.df$mean_score < 76.48) ~ "good",
    stats.df$icprb_bioregion_id == "sgv" & stats.df$mean_score >= 76.48 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "lnp" & stats.df$mean_score < 35.16 ~ "very_poor",
    stats.df$icprb_bioregion_id == "lnp" & (stats.df$mean_score >= 35.16 & stats.df$mean_score < 70.32) ~ "poor",
    stats.df$icprb_bioregion_id == "lnp" & (stats.df$mean_score >= 70.32 & stats.df$mean_score < 80.15) ~ "fair",
    stats.df$icprb_bioregion_id == "lnp" & (stats.df$mean_score >= 80.15 & stats.df$mean_score < 91.30) ~ "good",
    stats.df$icprb_bioregion_id == "lnp" & stats.df$mean_score >= 91.30 ~ "excellent",
    
    stats.df$icprb_bioregion_id == "pied" & stats.df$mean_score < 30.31 ~ "very_poor",
    stats.df$icprb_bioregion_id == "pied" & (stats.df$mean_score >= 30.31 & stats.df$mean_score < 60.62) ~ "poor",
    stats.df$icprb_bioregion_id == "pied" & (stats.df$mean_score >= 60.62 & stats.df$mean_score < 73.18) ~ "fair",
    stats.df$icprb_bioregion_id == "pied" & (stats.df$mean_score >= 73.18 & stats.df$mean_score < 81.83) ~ "good",
    stats.df$icprb_bioregion_id == "pied" & stats.df$mean_score >= 81.83 ~ "excellent",
    TRUE ~ "ERROR")
##  
# And the column numbers here need to be changed as well if you used more or less than 100 repetitions. 
##
  stats.df$mean_percent <- ((rowSums(stats.df[ , 103:202] == stats.df$mean_rating))/100)*100 
  stats.df$rating_compare <- stats.df$median_rating == stats.df$mean_rating

ChessieBIBI_final <- stats.df
```
**5. Final Dataframe**
Export the final data frame, **ChessieBIBI_final** as a CSV. Replace the text with your own.
```{r}
setwd("WORKING DIRECTORY")
write.csv(ChessieBIBI_final, file= "FILE NAME.csv")
```
