CalculatingES50.Rmd

---
title: "CalculatingES50"
author: "Natalie Posdaljian"
date: "8/12/2020"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Hurlbert's index, or ES50, is the expected number of distinct species in a random sample of 50 observations; it is a bias-indepenent indicator on marine biodiversity richness. It has the ability to somewhat correct for sampling bias.

Assumptions of ES50
- That indivdiuals are randomly distributed
- Sample size is sufficiently large
- Samples are taxonomically similar
- All samples have been taken in the same manner

```{r}
#loading robis
installed <- rownames(installed.packages())
if ( !("robis" %in% installed) ){
   if ( !("remotes" %in% installed) )install.packages("remotes")
   remotes::install_github("iobis/robis")
}

#load libraries
library(robis)
library('ggplot2')
library("rnaturalearth")
library("rnaturalearthdata")
library(roperators)
library(dplyr)
library(dggridR)
library(magrittr) 
library(vegan)
```

Load all occurences within a specificed polygon. (WARNING: using Casco Bay for now, but this will eventually be the MPA polygons)

```{r}
SpeciesOccurence <- occurrence(geometry = "POLYGON ((-70.2 43.5, -69.8 43.5, -69.8 43.9, -70.2 43.9, -70.2 43.5)) ")
```
Visualize and mutate the data
```{r}
##visualize the data
head(SpeciesOccurence)

#convert individual counts from character to numeric
SpeciesOccurence$individualCount <- as.numeric(SpeciesOccurence$individualCount)
SpeciesOccurence$individualCount[is.na(SpeciesOccurence$individualCount)] <- 1 #convert NANs to 1; I'm assuming that if it's listed, there was at least one count, even if it wasn't listed
SpeciesOccurence$Count <- 1 * SpeciesOccurence$individualCount

#function for extracting obis dates
extract_obis_date <- function(x = c("2016-01-02", "2016-01-03T05:06:07", "June 29, 1999")){
  as.Date(substring(x, 1, nchar("1234-56-78")), format = "%Y-%m-%d")
}

SpeciesOccurences$eventDate = extract_obis_date(SpeciesOccurence$eventDate) #change eventDate from character to date

#calculate the number of unique species
SpeciesCount <- aggregate(SpeciesOccurence$Count, by=list(Category=SpeciesOccurence$scientificName),FUN=sum)
head(SpeciesCount)
hist(SpeciesCount$x)
```

Calculate ES50 for all records
```{r}
ES50 = rarefy(SpeciesCount$x,50) #calculate ES50
print(ES50)
```

Calculate ES50 before and after a specified date
```{r}
MPA_established = as.Date(c("2007-06-22")) #what date was the MPA established

SpeciesOccurence_preMPA <- subset(SpeciesOccurence, SpeciesOccurence$eventDate < MPA_established) #Species Occurences before the MPA was established
SpeciesOccurence_postMPA <- subset(SpeciesOccurence, SpeciesOccurence$eventDate > MPA_established) #Species Occurences after the MPA was established

#Calculate unique species pre-MPA
SpeciesCount_preMPA <- aggregate(SpeciesOccurence_preMPA$Count, by=list(Category=SpeciesOccurence_preMPA$scientificName),FUN=sum)

#Calculate ES50 pre-MPA
ES50_preMPA = rarefy(SpeciesCount_preMPA$x,50) #calculate ES50
print(ES50_preMPA)

#Calculate unique species post-MPA
SpeciesCount_postMPA <- aggregate(SpeciesOccurence_postMPA$Count, by=list(Category=SpeciesOccurence_postMPA$scientificName),FUN=sum)

#Calculate ES50 post-MPA
ES50_postMPA = rarefy(SpeciesCount_postMPA$x,50) #calculate ES50
print(ES50_postMPA)
```