library(readr)
library(tidyverse)
library(vegan)
library(stringr)
library(ggplot2)
library(FSA)
#so that I can create a pdf using R.markdown on my PC, I need to install the following
devtools::install_github('yihui/tinytex')
tinytex::install_tinytex()

#import all bird data for 2016 and 2017
June_16 <- read_csv("NEON.D07.GRSM.DP1.10003.001.brd_countdata.2016-06.expanded.20180418T200555Z.csv")
may_17 <- read_csv("NEON.D07.GRSM.DP1.10003.001.brd_countdata.2017-05.expanded.20180418T200734Z.csv")
june_17 <- read_csv("NEON.D07.GRSM.DP1.10003.001.brd_countdata.2017-06.expanded.20180418T200639Z.csv")
may_18 <- read_csv("NEON.D07.GRSM.DP1.10003.001.brd_countdata.2018-05.expanded.20181206T143246Z.csv")

########################################
#DATA MANIPULATION
########################################
#combine all date from both years into one bird file and add a column for year from startDate
#str_sub is a function that extracts and replaces substrings from a character vector
brd_dat_all <- bind_rows(June_16, june_17, may_17, may_18)
brd_dat_all$Year <- str_sub(brd_dat_all$startDate, 1, 4)

#filter out the plotID values that are not replicated before and after the fire
brd_dat <- filter(brd_dat_all, plotID %in% c("GRSM_001","GRSM_002","GRSM_003","GRSM_004",
                  "GRSM_006","GRSM_007","GRSM_010","GRSM_013","GRSM_015", "GRSM_020"))

#convert plotID and taxonID to factors
brd_dat$plotID <- as.factor(brd_dat$plotID)
brd_dat$taxonID <- as.factor(brd_dat$taxonID)
summary(brd_dat)

#get rid of rows where targetTaxaPresent = N
brd_dat <- filter (brd_dat, targetTaxaPresent=="Y")

####################################
#CREATE DATA INPUT FOR VEGAN
####################################
#first filter by year to get the pre- and post-burn datasets separate
bird_16 <- filter(brd_dat, Year==2016)
bird_17 <- filter(brd_dat, Year==2017)
bird_18 <- filter(brd_dat, Year==2018)

#then create a matrix of plotID by species using the piping function (%>%) 
#to carry out multiple steps at a time using the date in the step before

#number of each species detected per plot
brd_matrix_16 <- bird_16 %>% group_by(plotID, taxonID) %>%
  summarise(count=n()) %>%
  spread(taxonID,count) 
brd_matrix_16[is.na(brd_matrix_16)] <- 0 #convert NAs to 0
#vegan package does not like first column of plotID, so get rid of it!
brd_matrix_16 <- brd_matrix_16[,-1] 
brd_matrix_16 #look at it - does it look ok?

#repeat for 2017 and 2018
brd_matrix_17 <- bird_17 %>% group_by(plotID, taxonID) %>%
  summarise(count=n()) %>%
  spread(taxonID,count)
brd_matrix_17[is.na(brd_matrix_17)] <- 0 
brd_matrix_17 <- brd_matrix_17[,-1]  
brd_matrix_17 #look at it - does it look ok?

brd_matrix_18 <- bird_18 %>% group_by(plotID, taxonID) %>%
  summarise(count=n()) %>%
  spread(taxonID,count)
brd_matrix_18[is.na(brd_matrix_18)] <- 0 
brd_matrix_18 <- brd_matrix_18[,-1]  
brd_matrix_18 #look at it - does it look ok?

#now we are ready to calculate some diversity indices for each plot
#http://cc.oulu.fi/~jarioksa/softhelp/vegan/html/diversity.html

#calculate shannon diversity and richness for each plot in each year

shd_16 <- diversity(brd_matrix_16, index = "shannon", MARGIN = 1, base = exp(1))
shd_17 <- diversity(brd_matrix_17, index = "shannon", MARGIN = 1, base = exp(1))
shd_18 <- diversity(brd_matrix_18, index = "shannon", MARGIN = 1, base = exp(1))

rich_16 <- specnumber(brd_matrix_16, MARGIN = 1)
rich_17 <- specnumber(brd_matrix_17, MARGIN = 1)
rich_18 <- specnumber(brd_matrix_18, MARGIN = 1)

####################################
#Create dataset for running ANOVA to test hypotheses about changes in richness or diversity pre- and post-fire
###################################

#create a new dataframe for each year (x_16 and x_17) with the above vectors from vegan and a few more created below
plotID <- c("GRSM_001","GRSM_002","GRSM_003","GRSM_004","GRSM_006","GRSM_007","GRSM_010","GRSM_013","GRSM_015", "GRSM_020")
burn <- c("y","n","y","n","n","y","n","n","n","n")
year <- rep(2016, 10)
year2 <- rep(2017, 10)
year3 <- rep(2018, 10)
x_16 <- data.frame(plotID, year, burn, shd_16, rich_16)
#make the column names the same in both so we can combine
colnames(x_16)<- c("plotID","year","burn","shd", "rich") 
x_17 <- data.frame(plotID, year2, burn, shd_17, rich_17)
colnames(x_17)<- c("plotID","year", "burn","shd", "rich")
x_18 <- data.frame(plotID, year3, burn, shd_18, rich_18)
colnames(x_18)<- c("plotID","year", "burn","shd", "rich")

x <- bind_rows(x_16, x_17, x_18) #combine x_16 and x_17 and x_18 into one dataframe
x$year <- as.factor(x$year) # make year a factor

#run ANOVA to see if treatment (burn), year, or an interaction of burn and year affect richness and shannon diversity
fit_rich <- aov(rich~burn+year+burn*year, x)
summary(fit_rich)

fit_shd <- aov(shd~burn+year+burn*year, x)
summary(fit_shd)

#to calculate SD
rich_means <- Summarize(rich ~ burn + year, data=x)
diversity_means <- Summarize(shd ~ burn + year, data=x)

#plot the data to visualize differences
ggplot(rich_means, aes(x=burn, y=mean, fill=year))+
  geom_bar(stat="identity", position=position_dodge()) +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=0.2, position=position_dodge(0.9))+
  xlab("Burn") +
  ylab("Bird Spp Richness")+
  theme_classic()

ggplot(diversity_means, aes(x=burn, y=mean, fill=year))+
  geom_bar(stat="identity", position=position_dodge()) +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=0.2, position=position_dodge(0.9))+
  xlab("Burn") +
  ylab("Bird Spp Diversity")+
  theme_classic()