---
title: "Vieques Macrobes Paper"
author: "Jason Baer & Jessica Carilli"
date: "2023-12-19"
output: html_document
---

Survival analysis - prepare data

```{r}
library(ggplot2)
library(forcats) 
library(dplyr)
library(gridExtra) 
library(grid)
library(RColorBrewer)
library(data.table) 
library(survival)
library(ggfortify) 
library(survminer) 
library(gtsummary) 
library(ggsurvfit)

#add data
dataC<-read.table("Coral_paper_dataset_Sept_2023.txt",header=T)

#filter corals that died
dead<-dataC %>%
  filter(overall_status==1)
dead<-as.data.frame(dead)

#filter corals that went missing
missing<-dataC %>%
  filter(overall_status==2)
missing<-as.data.frame(missing)

#filter live corals
live<-dataC %>%
  filter(!is.na(overall_status),overall_status==0)
live<-as.data.frame(live)

#filter to keep first entry per dead coral 
deadfirst<-dead %>%
  group_by(Nubbin_ID,Treatment,Plate_No,Species,Project_stage,site_collect_simple,overall_status)%>%
  summarize(weeks=min(weeks,na.rm=T)) 
deadfirst<-as.data.frame(deadfirst)

#filter down to one entry per missing coral: penultimate time at which the missing coral was alive
missfirst<-missing %>%
  group_by(Nubbin_ID,Treatment,Plate_No,Species,Project_stage,site_collect_simple,overall_status)%>%
  summarize(weeks=min(weeks,na.rm=T))
missfirst<-as.data.frame(missfirst)

#filter down to one entry per living coral (latest date at which coral was observed live)  
livelast<-live %>%
  group_by(Nubbin_ID,Treatment,Plate_No,Species,Project_stage,site_collect_simple,overall_status)%>%
  summarize(weeks=max(weeks,na.rm=T))
livelast<-as.data.frame(livelast)

#now find Nubbin_IDs of corals that died and should be removed from the live category
matches<-unique(grep(paste(dead$Nubbin_ID,collapse="|"),livelast$Nubbin_ID,value=TRUE))
alsomatches<-livelast$Nubbin_ID %in% matches
livefiltered<-livelast[!alsomatches,]
check<-unique(grep(paste(dead$Nubbin_ID,collapse="|"),livefiltered$Nubbin_ID,value=TRUE)) 
check #if empty = successfully removed double entries for dead corals
#now remove any missing corals from live
match2<-unique(grep(paste(missing$Nubbin_ID,collapse="|"),livefiltered$Nubbin_ID,value=TRUE))
alsomatch2<-livefiltered$Nubbin_ID %in% match2
livefiltered2<-livefiltered[!alsomatch2,]
check<-unique(grep(paste(missing$Nubbin_ID,collapse="|"),livefiltered2$Nubbin_ID,value=TRUE)) 
check #if empty = successfully removed double entries for missing corals
#can also run these to double check that there is only one entry per coral: dead, missing, or live
#check<-unique(grep(paste(deadfirst$Nubbin_ID,collapse="|"),livefiltered2$Nubbin_ID,value=TRUE)) 
#check<-unique(grep(paste(missfirst$Nubbin_ID,collapse="|"),livefiltered2$Nubbin_ID,value=TRUE)) 
#check<-unique(grep(paste(missfirst$Nubbin_ID,collapse="|"),deadfirst$Nubbin_ID,value=TRUE)) 

#now complete table for survival analysis
dataSurvival<-rbind(livefiltered2,deadfirst,missfirst)
dataSurvival$weeks<-as.numeric(dataSurvival$weeks) 
#convert status and Treatments into factors
dataSurvival$status<-as.factor(dataSurvival$overall_status)
dataSurvival$Treatment<-as.factor(dataSurvival$Treatment)
#order dataset so that time is increasing
dataSurvivalsorted<-dataSurvival%>% arrange(weeks)

#save data as .csv 
write.csv(dataSurvivalsorted,file="dataSurvivalsorted.csv")
```

Kaplan Meier overall survivorship analysis and plot

```{r}
library(dplyr)
library(survival)
library(survminer)

dataCompete<- read.csv("dataSurvivalsorted.csv")
#make sure values are in correct format
dataCompete$weeks<-as.numeric(dataCompete$weeks) 
dataCompete$Treatment<-as.factor(dataCompete$Treatment)

#convert missing values from 2 to 1 because this plot is just comparing living vs. not
missing1<-which(dataCompete$status=="2")
dataCompete_survival<-dataCompete
dataCompete_survival$status[missing1]<-1

coral_surv <- survfit(Surv(weeks,status) ~ Treatment, data = dataCompete_survival) #
summary(coral_surv) #show survival table
result <- survdiff(Surv(weeks,status) ~ Treatment, data = dataCompete_survival)
result #show Chi-squared result testing for significant difference in survival by treatment

#survival plot
coral_surv_plot <- ggsurvplot(
  coral_surv,
  data=dataCompete_survival,
  size = 0.75,                 # change line size
  palette =
    c("darkcyan", "wheat3"),# custom color palettes
  conf.int = TRUE,          # Add confidence interval
  pval = TRUE,              # Add p-value
  risk.table = FALSE,
  legend = c("right"),
  legend.title = "Treatment",
  legend.labs =
    c("Ark", "Control"),    # Change legend labels
  xlim = c(0, 82), 
  ylim = c(0,1),
  ggtheme = theme_light(base_size = 14),      
  xlab = "Time (Weeks since deployment)",
  ylab = "Proportion Coral Survival"
) 

coral_surv_plot
ggsave("Coral_surv.jpg", plot=coral_surv_plot$plot, width = 8, height = 6, units = "in", dpi = 1200)

```

Competing risks (death, falling off) analysis and plot

```{r}
library(dplyr)
library(survival)
library(survminer)
library(ggsurvfit)
library(cmprsk)

dataCompete<- read.csv("dataSurvivalsorted.csv")
#make sure values are in correct format
dataCompete$weeks<-as.numeric(dataCompete$weeks) 
dataCompete$status<-as.factor(dataCompete$status)
dataCompete$Treatment<-as.factor(dataCompete$Treatment)

coral_fit <- survfit(Surv(weeks, status) ~ Treatment, data = dataCompete)
summary(coral_fit) #provides risk table for each risk type

#test for significant differences in outcome between treatments using Gray's test
library(tidycmprsk)
cuminc_data <- cuminc(Surv(weeks, status) ~ Treatment, dataCompete)
cuminc_data #no significant differences (p=0.1) for missingness, but yes (p<0.001) for death

#competing risks plot
options("ggsurvfit.switch-color-linetype" = TRUE)
coral_surv_plot2 <- cuminc(Surv(weeks, status) ~ Treatment, dataCompete) %>%
  ggcuminc(outcome = c("1", "2")) +
  scale_linetype_manual(name = "Outcome", values = 1:2, labels = c("Dead", "Missing")) +
  scale_color_manual(values = c("darkcyan", "wheat3")) +
  scale_y_continuous(
    limits = c(0, 0.5),
    labels = scales::percent, 
    expand = c(0.01, 0)) +
  labs(x = "Time (Weeks since deployment)", y = "% Coral Dead or Missing")

coral_surv_plot2
ggsave("Coral_surv_compete.jpg", coral_surv_plot2, width = 8, height = 6, units = "in", dpi = 1200)

```

Plot number of corals that fell off plates or died at each time point

```{r}
library(ggplot2)
library(dplyr)

#add data
dataC<-read.table("Coral_paper_dataset_Sept_2023.txt",header=T)
dataC$Date_approx<-as.Date(dataC$Date_approx)

#filter corals that died
dead<-dataC %>%
  filter(overall_status==1)
dead<-as.data.frame(dead)

#filter corals that went missing
missing<-dataC %>%
  filter(overall_status==2)
missing<-as.data.frame(missing)

#filter keep first entry per missing coral
missoccurences<-missing %>%
  group_by(Nubbin_ID,Treatment,Project_stage,overall_status)%>%
  summarize(missing=min(Date_approx,na.rm=T))
missfirst<-as.data.frame(missfirst)

#count number of missing entries per time
misscounts<-missoccurences %>%
  group_by(Treatment,Project_stage,missing)%>%
  summarize(count=n_distinct(Nubbin_ID))
misscounts<-as.data.frame(misscounts)
#add zero entry for May 2022 Arks
aM22miss<-c("Ark",1,"2022-05-01",0)
misscounts<-rbind(misscounts,aM22miss)
misscounts$count<-as.numeric(misscounts$count)

#plot
custom_breaks <- as.Date(c("2022-02-01", "2022-05-01","2022-08-01","2022-12-01","2023-06-01"))
missing_plot<- ggplot(misscounts, aes(x = missing, y = count, fill = Treatment, alpha = factor(Project_stage))) +
  geom_bar(stat = "identity", position = "dodge") +
  scale_fill_manual(values = c("Ark" = "darkcyan","Control"="wheat3")) +
  scale_alpha_manual(values = c("1"=1,"2" = 0.5)) +
  labs(
    x = "Monitoring Date",
    y = "Number of New Missing Corals",
    fill = "Treatment",
    alpha = "Stage"
  ) +
  scale_x_date(breaks = custom_breaks,date_labels = "%b %Y") +
  theme_light(base_size = 14)
missing_plot + ylim(0,15)
ggsave("missing_plot.jpg", missing_plot, width = 8, height = 6, units = "in", dpi = 1200)

#test significant changes in rate of missing corals with time
arktime1<-misscounts%>%
  filter(Treatment=="Ark"&Project_stage==1)
  arktime1<-as.data.frame(arktime1)
missarktime1<-lm(arktime1$count~arktime1$missing) 
summary(missarktime1) #p-value: 0.0539
arktime2<-misscounts%>%
  filter(Treatment=="Ark"&Project_stage==2)
  arktime2<-as.data.frame(arktime2)
missarktime2<-lm(arktime2$count~arktime2$missing) 
summary(missarktime2) #p-value: 0.2611
conttime1<-misscounts%>%
  filter(Treatment=="Control"&Project_stage==1)
  conttime1<-as.data.frame(conttime1)
missconttime1<-lm(conttime1$count~conttime1$missing) 
summary(missconttime1) #p-value: 0.406
conttime2<-misscounts%>%
  filter(Treatment=="Control"&Project_stage==2)
  conttime2<-as.data.frame(conttime2)
missconttime2<-lm(conttime2$count~conttime2$missing) 
summary(missconttime2) #p-value: 0.56

#do the same for death
#filter to keep first entry per dead coral 
deadoccurences<-dead %>%
  group_by(Nubbin_ID,Treatment,Project_stage,overall_status)%>%
  summarize(died=min(Date_approx,na.rm=T)) 
deadoccurences<-as.data.frame(deadoccurences)

#count number of dead entries per time
deadcounts<-deadoccurences %>%
  group_by(Treatment,Project_stage,died)%>%
  summarize(count=n_distinct(Nubbin_ID))
deadcounts<-as.data.frame(deadcounts)

#plot
dead_plot<- ggplot(deadcounts, aes(x = died, y = count, fill = Treatment, alpha = factor(Project_stage))) +
  geom_bar(stat = "identity", position = "dodge") +
  scale_fill_manual(values = c("Ark" = "darkcyan","Control"="wheat3")) +
  scale_alpha_manual(values = c("1"=1,"2" = 0.5)) +
  labs(
    x = "Monitoring Date",
    y = "Number of New Dead Corals",
    fill = "Treatment",
    alpha = "Stage"
  ) +
  scale_x_date(breaks = custom_breaks,date_labels = "%b %Y") +
  theme_light(base_size = 14)
dead_plot + ylim(0,15)
ggsave("dead_plot.jpg", dead_plot, width = 8, height = 6, units = "in", dpi = 1200)

#test significant changes in rate of dead corals with time
arktime1<-deadcounts%>%
  filter(Treatment=="Ark"&Project_stage==1)
  arktime1<-as.data.frame(arktime1)
deadarktime1<-lm(arktime1$count~arktime1$died) 
summary(deadarktime1) #p-value: 0.064
arktime2<-deadcounts%>%
  filter(Treatment=="Ark"&Project_stage==2)
  arktime2<-as.data.frame(arktime2)
deadarktime2<-lm(arktime2$count~arktime2$died) 
summary(deadarktime2) #p-value: 0.5032
conttime1<-deadcounts%>%
  filter(Treatment=="Control"&Project_stage==1)
  conttime1<-as.data.frame(conttime1)
deadconttime1<-lm(conttime1$count~conttime1$died) 
summary(deadconttime1) #p-value: 0.0995
conttime2<-deadcounts%>%
  filter(Treatment=="Control"&Project_stage==2)
  conttime2<-as.data.frame(conttime2)
deadconttime2<-lm(conttime2$count~conttime2$died) 
summary(deadconttime2) #p-value: 0.04881

```

Plot average living surface area and volume of coral per coral plate depending on treatment and time

```{r}
library(ggplot2)
library(dplyr)

#read coral data into R
dataC<-read.table("Coral_paper_dataset_Sept_2023.txt",header=T)

#assign morphologies to corals
branching<-c("Porites_porites","Acropora_cervicornis","Porites_furcata")
massive<-c("Orbicella_faveolata","Orbicella_franksii","Porites_astreoides","Siderastrea_radians","Siderastrea_siderea","Agaricia")
whichbranch<-which(dataC$Species %in% branching)
whichmass<-which(dataC$Species %in% massive)
dataC$Morphology<-rep(1,length(dataC$Species))
dataC$Morphology[whichbranch]<-"branching"
dataC$Morphology[whichmass]<-"massive"
dataC$surface_area<-rep(1,length(dataC$Species))
dataC$volume<-rep(1,length(dataC$Species))

#calculate surface area of massive/encrusting corals using Dome: π(h^2+r^2 ) 
dataC$surface_area[whichmass]<- ((pi)*((dataC$Height[whichmass]^2)+(((dataC$Max_horiz[whichmass]+dataC$X90_deg_horiz[whichmass])/2)^2)))

#calculate surface area of branching corals using Cylinder with top: 2πrh+πr^2
#Multiplied by adjustment factor of 0.44
dataC$surface_area[whichbranch]<-((4/3)*(pi)*((dataC$Height[whichbranch]/2)*(dataC$Max_horiz[whichbranch]/2)*(dataC$X90_deg_horiz[whichbranch]/2)))*(0.44)

#correct surface area for any partial mortality
dataC$surface_area_live<-dataC$surface_area*((100-dataC$partial_mortality.)/100)

#calculate volume of massive/encrusting corals using Dome: 1/6 πh(3r^2+h^2 )
dataC$volume[whichmass]<- (1/6)*(pi)*(dataC$Height[whichmass])*((3*(((dataC$Max_horiz[whichmass]+dataC$X90_deg_horiz[whichmass])/2)^2))+(dataC$Height[whichmass]^2))

#calculate volume of branching corals using Ellipse: 4/3 π(h/2*x/2*y/2)
dataC$volume[whichbranch]<-(4/3)*(pi)*((dataC$Height[whichbranch]/2)*(dataC$Max_horiz[whichbranch]/2)*(dataC$X90_deg_horiz[whichbranch]/2))

#correct for any partial mortality
dataC$volume_live<-dataC$volume*((100-dataC$partial_mortality.)/100)

#sum living coral tissue surface area per plate
dataPlateSA<-dataC %>%
  group_by(Treatment,Plate_No,approx_month) %>%
  summarize(livingSA=sum(surface_area_live,na.rm=TRUE))
dataPlateSA<-as.data.frame(dataPlateSA)
dataPlateSA<-na.omit(dataPlateSA)

#sum living coral volume per plate 
dataPlateV<-dataC %>%
  group_by(Treatment,Plate_No,approx_month) %>%
  summarize(livingV=sum(volume_live,na.rm=TRUE))
dataPlateV<-as.data.frame(dataPlateV)
dataPlateV<-na.omit(dataPlateV)

#stats - test for significant differences in living surface area/plate by treatment at each timepoint
  #3 months deployed
statsSA3mo<-dataPlateSA %>%
  filter(approx_month ==3) 
statsSA3mo <-as.data.frame(statsSA3mo)
  shapiro.test(statsSA3mo$livingSA[statsSA3mo$Treatment=="Ark"]) #normal
  shapiro.test(statsSA3mo$livingSA[statsSA3mo$Treatment=="Control"]) #normal
  t.test(statsSA3mo$livingSA~statsSA3mo$Treatment)
#t = 3.5581, df = 76.835, p-value = 0.000644
  #6 months deployed
statsSA6mo<-dataPlateSA %>%
  filter(approx_month ==6) 
statsSA6mo <-as.data.frame(statsSA6mo)
  shapiro.test(statsSA6mo$livingSA[statsSA6mo$Treatment=="Ark"]) #normal
  shapiro.test(statsSA6mo$livingSA[statsSA6mo$Treatment=="Control"]) #not normal
  #square-root-transform to see if normality results
  statsSA6mo$livingSA_sq<-sqrt(statsSA6mo$livingSA)
  shapiro.test(statsSA6mo$livingSA_sq[statsSA6mo$Treatment=="Control"]) #normal
  shapiro.test(statsSA6mo$livingSA_sq[statsSA6mo$Treatment=="Ark"]) #normal
  t.test(statsSA6mo$livingSA_sq~statsSA6mo$Treatment)
#t = 3.6261, df = 62.346, p-value = 0.0005806
  #9 months deployed
statsSA9mo<-dataPlateSA %>%
  filter(approx_month ==9)
  statsSA9mo <-as.data.frame(statsSA9mo)
  shapiro.test(statsSA9mo$livingSA[statsSA9mo$Treatment=="Ark"]) #not normal
  shapiro.test(statsSA9mo$livingSA[statsSA9mo$Treatment=="Control"]) #not normal
  #square-root-transform to see if normality results
  statsSA9mo$livingSA_sq<-sqrt(statsSA9mo$livingSA)
  shapiro.test(statsSA9mo$livingSA_sq[statsSA9mo$Treatment=="Control"]) #normal
  shapiro.test(statsSA9mo$livingSA_sq[statsSA9mo$Treatment=="Ark"]) #normal
  t.test(statsSA9mo$livingSA_sq~statsSA9mo$Treatment)
#t = 3.5326, df = 30.426, p-value = 0.001336
  #12 months deployed
statsSA12mo<-dataPlateSA %>%
  filter(approx_month ==12)
  statsSA12mo <-as.data.frame(statsSA12mo)
  shapiro.test(statsSA12mo$livingSA[statsSA12mo$Treatment=="Ark"]) #not normal
  shapiro.test(statsSA12mo$livingSA[statsSA12mo$Treatment=="Control"]) #not normal
  #square-root-transform to see if normality results
  statsSA12mo$livingSA_sq<-sqrt(statsSA12mo$livingSA)
  shapiro.test(statsSA12mo$livingSA_sq[statsSA12mo$Treatment=="Control"]) #not normal
  shapiro.test(statsSA12mo$livingSA_sq[statsSA12mo$Treatment=="Ark"]) #not normal
  wilcox.test(statsSA12mo$livingSA_sq~statsSA12mo$Treatment)
#W = 1190, p-value = 1.756e-05 
  #19 months deployed
statsSA19mo<-dataPlateSA %>%
  filter(approx_month ==19)
  statsSA19mo <-as.data.frame(statsSA19mo)
  shapiro.test(statsSA19mo$livingSA[statsSA19mo$Treatment=="Ark"]) #not normal
  shapiro.test(statsSA19mo$livingSA[statsSA19mo$Treatment=="Control"]) #not normal
  #square-root-transform to see if normality results
  statsSA19mo$livingSA_sq<-sqrt(statsSA19mo$livingSA)
  shapiro.test(statsSA19mo$livingSA_sq[statsSA19mo$Treatment=="Control"]) #not normal
  shapiro.test(statsSA19mo$livingSA_sq[statsSA19mo$Treatment=="Ark"]) #normal
  wilcox.test(statsSA19mo$livingSA_sq~statsSA19mo$Treatment)
#W = 291, p-value = 0.0141

#stats - test for significant differences in living volume/plate by treatment at each timepoint
  #3 months deployed
coral3v <- dataPlateV %>% 
  filter(approx_month==3)
coral3v <-as.data.frame(coral3v)
  shapiro.test(coral3v$livingV[coral3v$Treatment=="Ark"]) #not normal
  shapiro.test(coral3v$livingV[coral3v$Treatment=="Control"]) #not normal
  #square-root-transform to see if normality results
  coral3v$livingV_sq<-sqrt(coral3v$livingV)
  shapiro.test(coral3v$livingV_sq[coral3v$Treatment=="Ark"]) #normal
  shapiro.test(coral3v$livingV_sq[coral3v$Treatment=="Control"])#normal
  t.test(coral3v$livingV_sq ~ coral3v$Treatment)
#t = 4.1523, df = 69.489, p-value = 9.204e-05
  #6 months deployed
coral6v <- dataPlateV %>% 
  filter(approx_month==6)
coral6v <-as.data.frame(coral6v)
  shapiro.test(coral6v$livingV[coral6v$Treatment=="Ark"]) #not normal
  shapiro.test(coral6v$livingV[coral6v$Treatment=="Control"]) #not normal
  #square-root-transform to see if normality results
  coral6v$livingV_sq<-sqrt(coral6v$livingV)
  shapiro.test(coral6v$livingV_sq[coral6v$Treatment=="Ark"]) #normal
  shapiro.test(coral6v$livingV_sq[coral6v$Treatment=="Control"]) #normal
  t.test(coral6v$livingV_sq ~ coral6v$Treatment)
#t = 3.3104, df = 66.422, p-value = 0.001508
  #9 months deployed
coral9v <- dataPlateV %>% 
  filter(approx_month==9)
coral9v <-as.data.frame(coral9v)
  shapiro.test(coral9v$livingV[coral9v$Treatment=="Ark"]) #not normal
  shapiro.test(coral9v$livingV[coral9v$Treatment=="Control"]) #not normal
  #square-root-transform to see if normality results
  coral9v$livingV_sq<-sqrt(coral9v$livingV)
  shapiro.test(coral9v$livingV_sq[coral9v$Treatment=="Ark"]) #normal
  shapiro.test(coral9v$livingV_sq[coral9v$Treatment=="Control"]) #not normal
  wilcox.test(coral9v$livingV_sq ~ coral9v$Treatment)
#W = 318, p-value = 0.001041
  #12 months deployed
coral12v <- dataPlateV %>% 
  filter(approx_month==12)
coral12v <-as.data.frame(coral12v)
  shapiro.test(coral12v$livingV[coral12v$Treatment=="Ark"]) #not normal
  shapiro.test(coral12v$livingV[coral12v$Treatment=="Control"]) #not normal
  #square-root-transform to see if normality results
  coral12v$livingV_sq<-sqrt(coral12v$livingV)
  shapiro.test(coral12v$livingV_sq[coral12v$Treatment=="Ark"]) #not normal
  shapiro.test(coral12v$livingV_sq[coral12v$Treatment=="Control"]) #not normal
  wilcox.test(coral12v$livingV_sq ~ coral12v$Treatment)
#W = 1209, p-value = 7.333e-06
  #19 months deployed
coral19v <- dataPlateV %>% 
  filter(approx_month==19) 
coral19v <-as.data.frame(coral19v)
  shapiro.test(coral19v$livingV[coral19v$Treatment=="Ark"]) #not normal
  shapiro.test(coral19v$livingV[coral19v$Treatment=="Control"]) #not normal
  #square-root-transform to see if normality results
  coral19v$livingV_sq<-sqrt(coral19v$livingV)
  shapiro.test(coral19v$livingV_sq[coral19v$Treatment=="Ark"]) #not normal
  shapiro.test(coral19v$livingV_sq[coral19v$Treatment=="Control"]) #not normal
  wilcox.test(coral19v$livingV_sq ~ coral19v$Treatment)
#W = 291, p-value = 0.0141

#Plots - surface area living coral
Plate_SA <- dataPlateSA %>% ggplot(aes(x = as.factor(approx_month), y = livingSA, color = Treatment, group = Treatment)) +
  stat_summary(
    fun.data = mean_se,  # Calculates mean and standard error
    geom = "errorbar",   # Draws error bars
    width = 0.2         # Width of error bars
  ) +
  stat_summary(
    fun.y = "mean",  # Calculates the mean
    geom = "point",  # Draws a point for the mean
    size = 3,  # Adjust the size of the point
  ) +
  stat_summary(
    fun.y = "mean",  # Calculates the mean
    geom = "line",  # Draws a line between points
    size = 0.5,  # Adjust the size of the lines
  ) +
  theme_test(base_size = 14) +
  xlab("Months Deployed") +
  ylab(expression(Average~Living~Coral~Tissue~(cm^{"2"})~per~Plate)) +
  scale_color_manual(values = c("darkcyan", "wheat4")) +
  labs(fill = "Treatment") +
 # facet_wrap(~ Treatment, scales = "free") + 
  coord_cartesian(ylim = c(0, 800))
Plate_SA
ggsave("Plate_surf_area.jpg", Plate_SA, width = 8, height = 6, units = "in", dpi = 1200)

#volume living coral
Plate_Vol <- dataPlateV %>% ggplot(aes(x = as.factor(approx_month), y = livingV, color = Treatment, group = Treatment)) +
  stat_summary(
    fun.data = mean_se,  # Calculates mean and standard error
    geom = "errorbar",   # Draws error bars
    width = 0.2         # Width of error bars
  ) +
  stat_summary(
    fun.y = "mean",  # Calculates the mean
    geom = "point",  # Draws a point for the mean
    size = 3,  # Adjust the size of the point
  ) +
  stat_summary(
    fun.y = "mean",  # Calculates the mean
    geom = "line",  # Draws a line between points
    size = 0.5,  # Adjust the size of the lines
  ) +
  theme_test(base_size = 14) +
  xlab("Months Deployed") +
  ylab(expression(Average~Living~Coral~Volume~(cm^{"3"})~per~Plate)) +
  scale_color_manual(values = c("darkcyan", "wheat4")) +
  labs(fill = "Treatment") +
 # facet_wrap(~ Treatment, scales = "free") + 
  coord_cartesian(ylim = c(0, 2000))
Plate_Vol
ggsave("Plate_vol.jpg", Plate_Vol, width = 8, height = 6, units = "in", dpi = 1200)

```

Turf algae on plates

```{r}
library(ggplot2)
library(dplyr)

#read coral plate data into R
dataC<-read.table("Coral_paper_dataset_Sept_2023.txt",header=T)

dataC$Treatment<-as.factor(dataC$Treatment)
dataC$Plate_No<-as.numeric(dataC$Plate_No)
dataC$approx_month<-as.numeric(dataC$approx_month)

#calculate mean and standard deviation of turf/macroalgae cover on coral plates by treatment and number of months plates had been deployed
dataPlateT<-dataC %>%
  group_by(Treatment,approx_month,Plate_No) %>% 
  summarize(Turf=mean(turf_percent,na.rm=TRUE)) %>% #make one entry/timepoint/plate
  group_by(Treatment,approx_month) %>%
  summarize(Turf=mean(Turf,na.rm=TRUE)) #calculate mean across plates
dataPlateT<-as.data.frame(dataPlateT)
dataPlateT<-na.omit(dataPlateT)
#add error
dataPlateTsd<-dataC %>%
  group_by(Treatment,approx_month,Plate_No) %>%
  summarize(Turfs=mean(turf_percent,na.rm=TRUE)) %>% #make one entry/timepoint/plate
  group_by(Treatment,approx_month) %>%
  summarize(TurfSD=sd(Turfs,na.rm=TRUE)) #calculate standard dev. across plates
dataPlateTsd<-as.data.frame(dataPlateTsd)
TurfSD<-dataPlateTsd$TurfSD
dataPlateT<-cbind(dataPlateT,TurfSD)
colnames(dataPlateT)[4]="TurfSD"

#Plot
Plate_turf <- dataPlateT %>% ggplot(aes(x = as.factor(approx_month), y = Turf, color = Treatment, group = Treatment)) +
  geom_point()+
  geom_line()+
  geom_errorbar(aes(ymin = (Turf-TurfSD), ymax = (Turf+TurfSD)), width = 0.2) +
  theme_test(base_size = 14) +
  xlab("Months Deployed") +
  ylab("% Cover Turf Algae per Plate") +
  scale_color_manual(values = c("cyan4", "wheat4")) +
  labs(color = "Treatment") 
  #ylim(0, 100) 
  #facet_wrap(~ Treatment, scales = "free")
Plate_turf
ggsave("Turf_plate.jpg", Plate_turf, width = 8, height = 6, units = "in", dpi = 1200)

#test for statistical differences in turf coverage for coral plates deployed for similar amounts of time between treatments

  #turf 3 mo deployment
turf3<-dataC %>% 
  filter(approx_month==3) %>%
  group_by(Treatment,Plate_No) %>%
  summarize(Turf=mean(turf_percent,na.rm=TRUE))
  shapiro.test(turf3$Turf[turf3$Treatment=="Ark"]) #normal
  shapiro.test(turf3$Turf[turf3$Treatment=="Control"]) #not normal
  #square-root-transform to see if normality results
  turf3$Turf_sq<-sqrt(turf3$Turf)
  shapiro.test(turf3$Turf_sq[turf3$Treatment=="Ark"]) #normal
  shapiro.test(turf3$Turf_sq[turf3$Treatment=="Control"]) #not normal
wilcox.test(turf3$Turf_sq~turf3$Treatment) 
#W = 33.5, p-value = 1.256e-13
  #turf 6 mo deployment
turf6<-dataC %>% 
  filter(approx_month==6) %>%
  group_by(Treatment, Plate_No) %>%
  summarize(Turf=mean(turf_percent,na.rm=TRUE))
  shapiro.test(turf6$Turf[turf6$Treatment=="Ark"]) #not normal
  shapiro.test(turf6$Turf[turf6$Treatment=="Control"]) #not normal
  #square-root-transform to see if normality results
  turf6$Turf_sq<-sqrt(turf6$Turf)
  shapiro.test(turf6$Turf_sq[turf6$Treatment=="Ark"]) #not normal
  shapiro.test(turf6$Turf_sq[turf6$Treatment=="Control"]) #not normal
wilcox.test(turf6$Turf_sq~turf6$Treatment) 
#W = 62, p-value = 3.648e-12
  #turf 9 mo deployment
turf9<-dataC %>% 
  filter(approx_month==9) %>%
  group_by(Treatment,Plate_No) %>%
  summarize(Turf=mean(turf_percent,na.rm=TRUE))
  shapiro.test(turf9$Turf[turf9$Treatment=="Ark"]) #normal
  shapiro.test(turf9$Turf[turf9$Treatment=="Control"]) #not normal
  #square-root-transform to see if normality results
  turf9$Turf_sq<-sqrt(turf9$Turf)
  shapiro.test(turf9$Turf_sq[turf9$Treatment=="Ark"]) #normal
  shapiro.test(turf9$Turf_sq[turf9$Treatment=="Control"]) #not normal
wilcox.test(turf9$Turf_sq~turf9$Treatment) 
#W = 0, p-value = 5.136e-08
  #turf 12 mo deployment
turf12<-dataC %>% 
  filter(approx_month==12) %>%
  group_by(Treatment, Plate_No) %>%
  summarize(Turf=mean(turf_percent,na.rm=TRUE))
  shapiro.test(turf12$Turf[turf12$Treatment=="Ark"]) #not normal
  shapiro.test(turf12$Turf[turf12$Treatment=="Control"]) #not normal
  #square-root-transform to see if normality results
  turf12$Turf_sq<-sqrt(turf12$Turf)
  shapiro.test(turf12$Turf_sq[turf12$Treatment=="Ark"]) #normal
  shapiro.test(turf12$Turf_sq[turf12$Treatment=="Control"]) #not normal
wilcox.test(turf12$Turf_sq~turf12$Treatment) 
#W = 39, p-value = 7.883e-13
  #turf 19 mo deployment
turf19<-dataC %>% 
  filter(approx_month==19) %>%
  group_by(Treatment, Plate_No) %>%
  summarize(Turf=mean(turf_percent,na.rm=TRUE))
  shapiro.test(turf19$Turf[turf19$Treatment=="Ark"]) #normal
  shapiro.test(turf19$Turf[turf19$Treatment=="Control"]) #not normal
  #square-root-transform to see if normality results
  turf19$Turf_sq<-sqrt(turf19$Turf)
  shapiro.test(turf19$Turf_sq[turf19$Treatment=="Ark"]) #normal
  shapiro.test(turf19$Turf_sq[turf19$Treatment=="Control"]) #not normal
wilcox.test(turf19$Turf_sq~turf19$Treatment) #W = 0, p-value = 4.723e-08

#test for changes in turf with time for each treatment (After time 0)
#Arks
turftimeA<-dataC %>%
  filter(Treatment == "Ark" & approx_month !=0) %>%
  group_by(approx_month,Plate_No) %>%
  summarize(Turf=mean(turf_percent,na.rm=TRUE))
turftimeA<-as.data.frame(turftimeA)
turfarktrend<-lm(turftimeA$Turf~turftimeA$approx_month)
summary(turfarktrend) 
#sig reducing trend with time F-statistic: 6.392 on 1 and 158 DF,  p-value: 0.01245
#Controls
turftimeC<-dataC %>%
  filter(Treatment == "Control" & approx_month!=0)%>%
  group_by(approx_month,Plate_No) %>%
  summarize(Turf=mean(turf_percent,na.rm=TRUE))
turftimeC<-as.data.frame(turftimeC)
turfconttrend<-lm(turftimeC$Turf~turftimeC$approx_month)
summary(turfconttrend)
#F-statistic: 2.557 on 1 and 152 DF,  p-value: 0.1119

```

Turf algae on plates - compare plates deployed with and without ARMS units

```{r}
library(ggplot2)
library(dplyr)

#read coral data into R
dataC<-read.table("Coral_paper_dataset_Sept_2023.txt",header=T)

dataC$Treatment<-as.factor(dataC$Treatment)
dataC$Project_stage<-as.factor(dataC$Project_stage)
dataC$Plate_No<-as.numeric(dataC$Plate_No)
dataC$approx_month<-as.numeric(dataC$approx_month)

#calculate mean and standard deviation of turf/macroalgae cover on coral plates by treatment, project stage (1 without, 2 with ARMS units) and number of months plates had been deployed
dataPlateT2<-dataC %>%
  group_by(Treatment,approx_month,Plate_No,Project_stage) %>%
  summarize(Turf=mean(turf_percent,na.rm=TRUE)) %>% #make one entry/timepoint/plate
  group_by(Treatment,approx_month,Project_stage) %>%
  summarize(Turf=mean(Turf,na.rm=TRUE)) #calculate mean across plates
dataPlateT2<-as.data.frame(dataPlateT2)
dataPlateT2<-na.omit(dataPlateT2)
#add error for turf
dataPlateTsd2<-dataC %>%
  group_by(Treatment,approx_month,Plate_No, Project_stage) %>% 
  summarize(Turf=mean(turf_percent,na.rm=TRUE)) %>% #make one entry/timepoint/plate
  group_by(Treatment,approx_month,Project_stage) %>%
  summarize(TurfSD=sd(Turf,na.rm=TRUE)) #calculate standard dev. across plates
dataPlateTsd2<-as.data.frame(dataPlateTsd2) 
TurfSD2<-dataPlateTsd2$TurfSD
dataPlateT2<-cbind(dataPlateT2,TurfSD2)
colnames(dataPlateT2)[5]="TurfSD"

#plot
Plate_turf5 <- dataPlateT2 %>% filter(Treatment == "Ark" & approx_month != 9 & approx_month != 12 & approx_month != 19) %>% ggplot(aes(x = as.factor(approx_month), y = Turf, color = factor(Project_stage), group = Project_stage)) +
  geom_point()+
  geom_line()+
  geom_errorbar(aes(ymin = (Turf-TurfSD), ymax = (Turf+TurfSD)), width = 0.2) +
  theme_test(base_size = 14) +
  xlab("Months Deployed") +
  ylab("% Cover Turf Algae per Plate") +
  scale_color_manual(values = c("black", "cyan4"),labels=c("- ARMS","+ ARMS")) +
  labs(color = "ARMS") 
  #ylim(0, 100) #+
  #facet_wrap(~ Phase_1, scales = "free")
Plate_turf5
ggsave("Turf_ARMS.jpg", Plate_turf5, width = 8, height = 6, units = "in", dpi = 1200)

#test for differences in turf for plates deployed over similar timeframes in Phase 1 (without ARMS) & 2 (with ARMS) at the Arks
  #3 mo deployment
turfphaseA3<-dataC %>%
  filter(Treatment == "Ark" & approx_month==3) %>%
  group_by(Project_stage, Plate_No) %>%
  summarize(Turf=mean(turf_percent,na.rm=TRUE))
turfphaseA3<-as.data.frame(turfphaseA3)
  shapiro.test(turfphaseA3$Turf[turfphaseA3$Project_stage==1]) #normal
  shapiro.test(turfphaseA3$Turf[turfphaseA3$Project_stage==2]) #normal
t.test(turfphaseA3$Turf~turfphaseA3$Project_stage) 
#t = -0.1812, df = 37.985, p-value = 0.8572
  #6 mo deployment
turfphaseA<-dataC %>%
  filter(Treatment == "Ark" & approx_month==6) %>%
  group_by(Project_stage, Plate_No) %>%
  summarize(Turf=mean(turf_percent,na.rm=TRUE))
turfphaseA<-as.data.frame(turfphaseA)
  shapiro.test(turfphaseA$Turf[turfphaseA$Project_stage==1]) #not normal
  shapiro.test(turfphaseA$Turf[turfphaseA$Project_stage==2]) #normal
  #square-root-transform to see if normality results
  turfphaseA$Turf_sq<-sqrt(turfphaseA$Turf)
  shapiro.test(turfphaseA$Turf_sq[turfphaseA$Project_stage==1]) #not normal
  shapiro.test(turfphaseA$Turf_sq[turfphaseA$Project_stage==2]) #normal
wilcox.test(turfphaseA$Turf~turfphaseA$Project_stage) 
#W = 318.5, p-value = 0.001348

#test for differences in turf for plates deployed for 6 months in Phase 1 (without ARMS) & 2 (with ARMS) at the Control sites
  #6 mo deployment
turfphaseC<-dataC %>%
  filter(Treatment == "Control" & approx_month==6) %>%
  group_by(Project_stage, Plate_No) %>%
  summarize(Turf=mean(turf_percent,na.rm=TRUE))
turfphaseC<-as.data.frame(turfphaseC)
  shapiro.test(turfphaseC$Turf[turfphaseC$Project_stage==1]) #not normal
  shapiro.test(turfphaseC$Turf[turfphaseC$Project_stage==2]) #not normal
  #square-root-transform to see if normality results
  turfphaseC$Turf_sq<-sqrt(turfphaseC$Turf)
  shapiro.test(turfphaseC$Turf_sq[turfphaseC$Project_stage==1]) #not normal
  shapiro.test(turfphaseC$Turf_sq[turfphaseC$Project_stage==2]) #not normal
wilcox.test(turfphaseC$Turf~turfphaseC$Project_stage) 
#W = 93.5, p-value = 0.01517

```

Fish abundance and biomass

```{r}
library(dplyr)
library(ggplot2)
library(lubridate)
library(wesanderson)
library(rstatix)

fish<-read.csv("fish_paper.csv",header=T)
#set formats
fish$Date<-as.Date(fish$Date,format="%m/%d/%y")
fish$Site<-as.numeric(fish$Site)
fish$Survey<-as.numeric(fish$Survey)
fish$Tropic_role<-as.factor(fish$Trophic_role)

#test for significant change with time in fish biomass (all trophic roles)
biomass_fish<-fish %>%
  group_by(Treatment,Site,Survey,Date)%>%
  summarize(sum(Biomass))
biomass_fish<-as.data.frame(biomass_fish)
colnames(biomass_fish)<-c("Treatment","Site","Survey","Date","biomass")
ark_biomass<-filter(biomass_fish,Treatment=="Ark")
control_biomass<-filter(biomass_fish,Treatment=="Control")
biomass_time_ark<-aov(biomass~Date,data=ark_biomass)
summary(biomass_time_ark)
#Df Sum Sq Mean Sq F value   Pr(>F)    
#Date         1 1450.7  1450.7   21.79 0.000363 ***
#Residuals   14  932.1    66.6                     
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
biomass_time_control<-aov(biomass~Date,data=control_biomass)
summary(biomass_time_control)
#Df Sum Sq Mean Sq F value Pr(>F)
#Date         1   21.3   21.27    0.17  0.689
#Residuals    9 1122.6  124.74 

#test for significant differences in fish biomass within timepoints with at least 3 surveys/treatment (Aug 2022 and June 2023)
#first filter data
aug_biomassT<-biomass_fish%>% 
  filter(Date=="2022-08-28")%>%
  group_by(Treatment,Site,Survey)%>%
  summarize(Biomass=sum(biomass))
jun23_biomassT<-biomass_fish%>% 
  filter(Date=="2023-06-04")%>%
  group_by(Treatment,Site,Survey)%>%
  summarize(Biomass=sum(biomass))

#test for normality
#Aug 2022
shapiro.test1<-aug_biomassT %>%
  filter(Treatment=="Ark")
shapiro.test(shapiro.test1$Biomass) #normal
shapiro.test2<-aug_biomassT %>%
  filter(Treatment=="Control")
shapiro.test(shapiro.test2$Biomass) #normal
#June 2023
shapiro.test3<-jun23_biomassT %>%
  filter(Treatment=="Ark")
shapiro.test(shapiro.test3$Biomass) #normal
shapiro.test4<-jun23_biomassT %>%
  filter(Treatment=="Control")
shapiro.test(shapiro.test4$Biomass) #normal

#test for significant differences by treatment
t.test(Biomass~Treatment,data=aug_biomassT)
#not signif: t = 1.9655, df = 2.5238, p-value = 0.1612
t.test(Biomass~Treatment,data=jun23_biomassT)
#not signif t = 0.7256, df = 4.9646, p-value = 0.5008

#test for significant change with time in fish abundance (all trophic roles)
abundance_fish<-fish %>%
  group_by(Treatment,Site,Survey,Date)%>%
  summarize(sum(Abundance))
abundance_fish<-as.data.frame(abundance_fish)
colnames(abundance_fish)<-c("Treatment","Site","Survey","Date","Abundance")
ark_abundance<-filter(abundance_fish,Treatment=="Ark")
control_abundance<-filter(abundance_fish,Treatment=="Control")
abundance_time_ark<-aov(Abundance~Date,data=ark_abundance)
summary(abundance_time_ark)
#Df Sum Sq Mean Sq F value   Pr(>F)    
#Date         1 312147  312147   15.01 0.00169 **
#Residuals   14  291204   20800                    
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
abundance_time_control<-aov(Abundance~Date,data=control_abundance)
summary(abundance_time_control)
#Df Sum Sq Mean Sq F value Pr(>F)
#Date         1   12081   12081   2.142  0.177
#Residuals    9 50753    5639 

#test for significant differences in fish abundance within timepoints with at least 3 surveys/treatment (Aug 2022 and June 2023)
#first filter data
aug_abundanceT<-abundance_fish%>% 
  filter(Date=="2022-08-28")%>%
  group_by(Treatment,Site,Survey)%>%
  summarize(Abundance=sum(Abundance))
jun23_abundanceT<-abundance_fish%>% 
  filter(Date=="2023-06-04")%>%
  group_by(Treatment,Site,Survey)%>%
  summarize(Abundance=sum(Abundance))

#test for normality
#Aug 2022
shapiro.test(aug_abundanceT$Abundance[aug_abundanceT$Treatment=="Ark"]) #normal
shapiro.test(aug_abundanceT$Abundance[aug_abundanceT$Treatment=="Control"]) #normal
#June 2023
shapiro.test(jun23_abundanceT$Abundance[jun23_abundanceT$Treatment=="Ark"]) #normal
shapiro.test(jun23_abundanceT$Abundance[jun23_abundanceT$Treatment=="Control"]) #normal

#test for significant differences by treatment
t.test(Abundance~Treatment,data=aug_abundanceT)
#t = 16.671, df = 4.7537, p-value = 2.119e-05
#mean in group Ark mean in group Control 
#             420.2500              150.6667
t.test(Abundance~Treatment,data=jun23_abundanceT)
#not signif t = 0.68966, df = 4.708, p-value = 0.5229
#mean in group Ark mean in group Control 
#                  356                   317 

#plot fish biomass and abundance with time by treatment
custom_labels <- c("Primary", "Secondary", "Planktivore", "Piscivore")

#biomass
biomass_fish$Date<-factor(biomass_fish$Date) 

Fish_bio <- biomass_fish %>%
  ggplot(aes(x = Date, y = biomass, color = Treatment, group = Treatment)) +
  stat_summary(
    fun.data = mean_se,  # Calculates mean and standard error
    geom = "errorbar",   # Draws error bars
    width = 0.2         # Width of error bars
  ) +
  stat_summary(
    fun.y = "mean",  # Calculates the mean
    geom = "point",  # Draws a point for the mean
    size = 3  # Adjust the size of the point
  ) +
  stat_summary(
    fun.y = "mean",  # Calculates the mean
    geom = "line",  # Draws lines between points
    size = 0.5  # Adjust the size of the line
  ) +
  theme_test(base_size = 14) +
  xlab("Monitoring Event") +
  ylab("Average Biomass of Fish per Survey (kg)") +
  scale_color_manual(values = c("darkcyan", "wheat4")) +
  labs(fill = "Treatment") +
  coord_cartesian(ylim = c(0, 50)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_x_discrete(labels = custom_labels)
Fish_bio
ggsave("Fish_bio.jpg", Fish_bio, width = 8, height = 6, units = "in", dpi = 1200)

#abundance
abundance_fish$Date<-factor(abundance_fish$Date) 

Fish_numbers <- abundance_fish %>%
  ggplot(aes(x = Date, y = Abundance, color = Treatment, group = Treatment)) +
  stat_summary(
    fun.data = mean_se,  # Calculates mean and standard error
    geom = "errorbar",   # Draws error bars
    width = 0.2         # Width of error bars
  ) +
  stat_summary(
    fun.y = "mean",  # Calculates the mean
    geom = "point",  # Draws a point for the mean
    size = 3  # Adjust the size of the point
  ) +
  stat_summary(
    fun.y = "mean",  # Calculates the mean
    geom = "line",  # Draws lines between points
    size = 0.5  # Adjust the size of the line
  ) +
  theme_test(base_size = 14) +
  xlab("Monitoring Event") +
  ylab("Average Number of Fish per Survey") +
  scale_color_manual(values = c("darkcyan", "wheat4")) +
  labs(fill = "Treatment") +
  coord_cartesian(ylim = c(0, 600)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_x_discrete(labels = custom_labels)
Fish_numbers
ggsave("Fish_numbers.jpg", Fish_numbers, width = 8, height = 6, units = "in", dpi = 1200)

```

Fish trophic roles

```{r}
library(dplyr)
library(ggplot2)
library(lubridate)
library(wesanderson)
library(rstatix)

fish<-read.csv("fish_paper.csv",header=T)
#set formats
fish$Date<-as.Date(fish$Date,format="%m/%d/%y")
fish$Site<-as.numeric(fish$Site)
fish$Survey<-as.numeric(fish$Survey)
fish$Tropic_role<-as.factor(fish$Trophic_role)

#abundance
abundance_trophic<-fish %>%
  group_by(Treatment,Date,Site,Survey,Trophic_role)%>%
  summarize(sum(Abundance))
abundance_trophic<-as.data.frame(abundance_trophic)
colnames(abundance_trophic)<-c("Treatment","Date","Site","Survey","Trophic_role","Abundance")

#Treatments don't have same # of surveys at each timepoint
#Can plot representative results from one survey from Ark 2 and Control 2 at each timepoint
abundance_plot<-filter(abundance_trophic,Site==2,Survey==1)
#convert trophic roles to factors and reorder so it goes bottom up
abundance_plot$Trophic_role <- factor(abundance_plot$Trophic_role, levels=c('Piscivore','Planktivore','Secondary','Primary'))

fish_trophic_time<-ggplot(abundance_plot,aes(x=Treatment,y=Abundance,fill=Trophic_role))+
  geom_bar(position="stack",stat="identity")+
  facet_grid(.~Date,labeller = function(x) format(x,'%b %Y'))+
  theme_test()+
  scale_fill_brewer(palette="Blues",direction = -1)+
  scale_y_continuous(expand = expansion(mult=c(0,0)),limits = c(0,500))+
  theme(axis.text.x = element_text(angle=45, hjust = 1))+
  guides(fill = guide_legend(title = "Trophic role"))+
  labs(x="Treatment and Monitoring Event",y="Number of Fish")
fish_trophic_time

ggsave("numbers_trophic_time.jpg", fish_trophic_time, width = 7, height = 6, units = "in", dpi = 1200)

#biomass
biomass_trophic_fish<-fish %>%
  group_by(Treatment,Date,Site,Survey,Trophic_role)%>%
  summarize(sum(Biomass))
biomass_trophic_fish<-as.data.frame(biomass_trophic_fish)
colnames(biomass_trophic_fish)<-c("Treatment","Date","Site","Survey","Trophic_role","Biomass")

#biomass to plot
biomass_trophic_plot<-filter(biomass_trophic_fish,Site==2,Survey==1)
#convert trophic roles to factors and reorder so it goes bottom up
biomass_trophic_plot$Trophic_role <- factor(biomass_trophic_plot$Trophic_role, levels=c('Piscivore','Planktivore','Secondary','Primary'))

bio_trophic_time<-ggplot(biomass_trophic_plot,aes(x=Treatment,y=Biomass,fill=Trophic_role))+
  geom_bar(position="stack",stat="identity")+
  facet_grid(.~Date,labeller = function(x) format(x,'%b %Y'))+
  theme_test()+
  scale_fill_brewer(palette="Blues",direction = -1)+
  scale_y_continuous(expand = expansion(mult=c(0,0)),limits=c(0,40))+
  theme(axis.text.x = element_text(angle=45, hjust = 1))+
  guides(fill = guide_legend(title = "Trophic role"))+
  labs(x="Treatment and Monitoring Event",y="Biomass of Fish (kg)")
bio_trophic_time

ggsave("biomass_trophic_time.jpg", bio_trophic_time, width = 7, height = 6, units = "in", dpi = 1200)

#test for significant differences in biomass of each trophic role for timepoints with at least 3 surveys
#first prep data for august
aug_biomass=filter(biomass_trophic_fish,Date=="2022-08-28")
aug_biomass_Pisc<-filter(aug_biomass,Trophic_role=="Piscivore")
#add zero observation for Control site 1 survey 1
cont2<-data.frame("Control","2022-08-28",1,1,"Piscivore",0)
colnames(cont2)<-c("Treatment","Date","Site","Survey","Trophic_role","Biomass")
aug_biomass_Pisc<-rbind(aug_biomass_Pisc,cont2)
aug_biomass_Plank<-filter(aug_biomass,Trophic_role=="Planktivore")
aug_biomass_Prim<-filter(aug_biomass,Trophic_role=="Primary")
aug_biomass_Sec<-filter(aug_biomass,Trophic_role=="Secondary")

#test normality
shapiro.test(aug_biomass_Pisc$Biomass[aug_biomass_Pisc$Treatment=="Ark"]) #normal
shapiro.test(aug_biomass_Pisc$Biomass[aug_biomass_Pisc$Treatment=="Control"]) #normal
shapiro.test(aug_biomass_Plank$Biomass[aug_biomass_Plank$Treatment=="Ark"]) #normal
shapiro.test(aug_biomass_Plank$Biomass[aug_biomass_Plank$Treatment=="Control"]) #normal
shapiro.test(aug_biomass_Prim$Biomass[aug_biomass_Prim$Treatment=="Ark"]) #not normal 
#try square-root transformation to see if becomes normal
aug_biomass_Prim$Biomass_sq<-sqrt(aug_biomass_Prim$Biomass)
shapiro.test(aug_biomass_Prim$Biomass_sq[aug_biomass_Prim$Treatment=="Ark"]) #normal
shapiro.test(aug_biomass_Prim$Biomass_sq[aug_biomass_Prim$Treatment=="Control"]) #normal
shapiro.test(aug_biomass_Sec$Biomass[aug_biomass_Sec$Treatment=="Ark"]) #normal
shapiro.test(aug_biomass_Sec$Biomass[aug_biomass_Sec$Treatment=="Control"]) #normal

#test for significant differences by treatment
t.test(aug_biomass_Pisc$Biomass~aug_biomass_Pisc$Treatment)
# data:  aug_biomass_Pisc$biomass by aug_biomass_Pisc$Treatment
# t = 11.667, df = 4.2934, p-value = 0.0002059
# sample estimates:
#   mean in group Ark mean in group Control 
#  15.5881653             0.6049642  
t.test(aug_biomass_Plank$Biomass~aug_biomass_Plank$Treatment)
# data:  aug_biomass_Plank$biomass by aug_biomass_Plank$Treatment
# t = -8.9991, df = 2.1135, p-value = 0.01017
# sample estimates:
#   mean in group Ark mean in group Control 
# 0.04009133            1.04264124 
t.test(aug_biomass_Prim$Biomass_sq~aug_biomass_Prim$Treatment)
#ns
t.test(aug_biomass_Sec$Biomass~aug_biomass_Sec$Treatment)
#ns

#prep June 2023 fish data
jun23_biomass=filter(biomass_trophic_fish,Date=="2023-06-04")
jun23_biomass_Pisc<-filter(jun23_biomass,Trophic_role=="Piscivore")
#need to add 1 observation of 0 for control
c12j<-data.frame("Control","2023-06-04",1,2,"Piscivore",0)
colnames(c12j)<-c("Treatment","Date","Site","Survey","Trophic_role","Biomass")
jun23_biomass_Pisc<-rbind(jun23_biomass_Pisc,c12j)
jun23_biomass_Plank<-filter(jun23_biomass,Trophic_role=="Planktivore")
#need to add 3 more observations of 0 for Arks...
a11j<-data.frame("Ark","2023-06-04",1,1,"Planktivore",0)
colnames(a11j)<-c("Treatment","Date","Site","Survey","Trophic_role","Biomass")
a21j<-data.frame("Ark","2023-06-04",2,1,"Planktivore",0)
colnames(a21j)<-c("Treatment","Date","Site","Survey","Trophic_role","Biomass")
a22j<-data.frame("Ark","2023-06-04",2,2,"Planktivore",0)
colnames(a22j)<-c("Treatment","Date","Site","Survey","Trophic_role","Biomass")
jun23_biomass_Plank<-rbind(jun23_biomass_Plank,a11j,a21j,a22j)
jun23_biomass_Prim<-filter(jun23_biomass,Trophic_role=="Primary")
jun23_biomass_Sec<-filter(jun23_biomass,Trophic_role=="Secondary")

#test for normality
shapiro.test(jun23_biomass_Pisc$Biomass[jun23_biomass_Pisc$Treatment=="Ark"]) #normal
shapiro.test(jun23_biomass_Pisc$Biomass[jun23_biomass_Pisc$Treatment=="Control"]) #normal
shapiro.test(jun23_biomass_Plank$Biomass[jun23_biomass_Plank$Treatment=="Ark"]) #not normal 
#try square-root transformation to see if becomes normal
jun23_biomass_Plank$Biomass_sq<-sqrt(jun23_biomass_Plank$Biomass)
shapiro.test(jun23_biomass_Plank$Biomass_sq[jun23_biomass_Plank$Treatment=="Ark"]) #not normal
shapiro.test(jun23_biomass_Plank$Biomass_sq[jun23_biomass_Plank$Treatment=="Control"]) #normal
shapiro.test(jun23_biomass_Prim$Biomass[jun23_biomass_Prim$Treatment=="Ark"]) #normal
shapiro.test(jun23_biomass_Prim$Biomass[jun23_biomass_Prim$Treatment=="Control"]) #not normal 
#try square-root transformation to see if becomes normal
jun23_biomass_Prim$Biomass_sq<-sqrt(jun23_biomass_Prim$Biomass)
shapiro.test(jun23_biomass_Prim$Biomass_sq[jun23_biomass_Prim$Treatment=="Control"]) #not normal 
shapiro.test(jun23_biomass_Sec$Biomass[jun23_biomass_Sec$Treatment=="Ark"]) #normal
shapiro.test(jun23_biomass_Sec$Biomass[jun23_biomass_Sec$Treatment=="Control"]) #normal

t.test(jun23_biomass_Pisc$Biomass~jun23_biomass_Pisc$Treatment)
#data:  jun23_biomass_Pisc$Biomass by jun23_biomass_Pisc$Treatment
#t = 3.0762, df = 4.9935, p-value = 0.02764
#sample estimates:
#    mean in group Ark mean in group Control 
#            24.632530              4.452348 
wilcox.test(jun23_biomass_Plank$Biomass_sq~jun23_biomass_Plank$Treatment) 
#W = 0, p-value = 0.0436
wilcox.test(jun23_biomass_Prim$Biomass_sq~jun23_biomass_Prim$Treatment)
#W = 0, p-value = 0.05714
t.test(jun23_biomass_Sec$Biomass~jun23_biomass_Sec$Treatment)
#t = -6.7305, df = 2.023, p-value = 0.02075
# mean in group Ark mean in group Control 
#             0.221647              4.732111 

#test for significant differences in abundance of each trophic role for timepoints with at least 3 surveys
#first prep data for august
aug_abundance=filter(abundance_trophic,Date=="2022-08-28")
aug_abundance_Pisc<-filter(aug_abundance,Trophic_role=="Piscivore")
#add zero observation for Control site 1 survey 1
cont2<-data.frame("Control","2022-08-28",1,1,"Piscivore",0)
colnames(cont2)<-c("Treatment","Date","Site","Survey","Trophic_role","Abundance")
aug_abundance_Pisc<-rbind(aug_abundance_Pisc,cont2)
aug_abundance_Plank<-filter(aug_abundance,Trophic_role=="Planktivore")
aug_abundance_Prim<-filter(aug_abundance,Trophic_role=="Primary")
aug_abundance_Sec<-filter(aug_abundance,Trophic_role=="Secondary")

#test normality
shapiro.test(aug_abundance_Pisc$Abundance[aug_abundance_Pisc$Treatment=="Ark"]) #not normal 
#try square-root transformation to see if becomes normal
aug_abundance_Pisc$Abundance_sq<-sqrt(aug_abundance_Pisc$Abundance)
shapiro.test(aug_abundance_Pisc$Abundance_sq[aug_abundance_Pisc$Treatment=="Ark"]) #not normal
shapiro.test(aug_abundance_Pisc$Abundance[aug_abundance_Pisc$Treatment=="Control"]) #normal
shapiro.test(aug_abundance_Plank$Abundance[aug_abundance_Plank$Treatment=="Ark"]) #normal
shapiro.test(aug_abundance_Plank$Abundance[aug_abundance_Plank$Treatment=="Control"]) #normal
shapiro.test(aug_abundance_Prim$Abundance[aug_abundance_Prim$Treatment=="Ark"]) #normal
shapiro.test(aug_abundance_Prim$Abundance[aug_abundance_Prim$Treatment=="Control"]) #normal
shapiro.test(aug_abundance_Sec$Abundance[aug_abundance_Sec$Treatment=="Ark"]) #normal
shapiro.test(aug_abundance_Sec$Abundance[aug_abundance_Sec$Treatment=="Control"]) #normal

#test for significant differences by treatment
wilcox.test(aug_abundance_Pisc$Abundance_sq~aug_abundance_Pisc$Treatment)
#W = 12, p-value = 0.0436
t.test(aug_abundance_Plank$Abundance~aug_abundance_Plank$Treatment)
#t = -10.787, df = 3.2401, p-value = 0.001198
# sample estimates:
#   mean in group Ark mean in group Control 
# 6.75000              60.33333 
t.test(aug_abundance_Prim$Abundance~aug_abundance_Prim$Treatment)
#ns
t.test(aug_abundance_Sec$Abundance~aug_abundance_Sec$Treatment)
#ns

#prep June 2023 fish data
jun23_abundance=filter(abundance_trophic,Date=="2023-06-04")
jun23_abundance_Pisc<-filter(jun23_abundance,Trophic_role=="Piscivore")
#need to add 1 observation of 0 for control
c12j<-data.frame("Control","2023-06-04",1,2,"Piscivore",0)
colnames(c12j)<-c("Treatment","Date","Site","Survey","Trophic_role","Abundance")
jun23_abundance_Pisc<-rbind(jun23_abundance_Pisc,c12j)
jun23_abundance_Plank<-filter(jun23_abundance,Trophic_role=="Planktivore")
#need to add 3 more observations of 0 for Arks...
a11j<-data.frame("Ark","2023-06-04",1,1,"Planktivore",0)
colnames(a11j)<-c("Treatment","Date","Site","Survey","Trophic_role","Abundance")
a21j<-data.frame("Ark","2023-06-04",2,1,"Planktivore",0)
colnames(a21j)<-c("Treatment","Date","Site","Survey","Trophic_role","Abundance")
a22j<-data.frame("Ark","2023-06-04",2,2,"Planktivore",0)
colnames(a22j)<-c("Treatment","Date","Site","Survey","Trophic_role","Abundance")
jun23_abundance_Plank<-rbind(jun23_abundance_Plank,a11j,a21j,a22j)
jun23_abundance_Prim<-filter(jun23_abundance,Trophic_role=="Primary")
jun23_abundance_Sec<-filter(jun23_abundance,Trophic_role=="Secondary")

#test for normality
shapiro.test(jun23_abundance_Pisc$Abundance[jun23_abundance_Pisc$Treatment=="Ark"]) #normal
shapiro.test(jun23_abundance_Pisc$Abundance[jun23_abundance_Pisc$Treatment=="Control"]) #not normal #try square-root transformation to see if becomes normal
jun23_abundance_Pisc$Abundance_sq<-sqrt(jun23_abundance_Pisc$Abundance)
shapiro.test(jun23_abundance_Pisc$Abundance_sq[jun23_abundance_Pisc$Treatment=="Control"]) #not normal
shapiro.test(jun23_abundance_Plank$Abundance[jun23_abundance_Plank$Treatment=="Ark"]) #not normal #try square-root transformation to see if becomes normal
jun23_abundance_Plank$Abundance_sq<-sqrt(jun23_abundance_Plank$Abundance)
shapiro.test(jun23_abundance_Plank$Abundance_sq[jun23_abundance_Plank$Treatment=="Ark"]) #not normal
shapiro.test(jun23_abundance_Plank$Abundance[jun23_abundance_Plank$Treatment=="Control"]) #normal
shapiro.test(jun23_abundance_Prim$Abundance[jun23_abundance_Prim$Treatment=="Ark"]) #normal
shapiro.test(jun23_abundance_Prim$Abundance[jun23_abundance_Prim$Treatment=="Control"]) #normal
shapiro.test(jun23_abundance_Sec$Abundance[jun23_abundance_Sec$Treatment=="Ark"]) #normal
shapiro.test(jun23_abundance_Sec$Abundance[jun23_abundance_Sec$Treatment=="Control"]) #normal

wilcox.test(jun23_abundance_Pisc$Abundance_sq~jun23_abundance_Pisc$Treatment)
#W = 12, p-value = 0.04975
#note Wilcox test not ideal with this type of dataset
wilcox.test(jun23_abundance_Plank$Abundance_sq~jun23_abundance_Plank$Treatment) 
#W = 0, p-value = 0.0436
t.test(jun23_abundance_Sec$Abundance~jun23_abundance_Sec$Treatment)
#t = -15.524, df = 4.9684, p-value = 2.116e-05
#mean in group Ark mean in group Control 
#              18.0000              117.3333 
t.test(jun23_abundance_Prim$Abundance~jun23_abundance_Prim$Treatment)
#t = -3.6125, df = 3.0074, p-value = 0.03629
#mean in group Ark mean in group Control 
# 31.25000              98.33333 

#representative plot - June 2023 data - biomass
custom_labels <- c("Primary", "Secondary", "Planktivore", "Piscivore")
jun23_biomass$Trophic_role <- factor(jun23_biomass$Trophic_role, levels = custom_labels)

jun<-ggplot(jun23_biomass, aes(x = Trophic_role, y = Biomass, fill = Treatment)) +
  geom_boxplot() +
  #facet_wrap(~ Trophic_role, scales = "free") +
  scale_fill_manual(values = c("Ark" = "darkcyan", "Control" = "wheat3")) +
  xlab("Trophic Category") +
  ylab("Fish Biomass (kg)") +
  theme_light(base_size = 14)
jun
ggsave("jun_trophic_biomass.jpg", jun, width = 8, height = 6, units = "in", dpi = 1200)

#representative plot - June 2023 data - abundance
custom_labels <- c("Primary", "Secondary", "Planktivore", "Piscivore")
jun23_abundance$Trophic_role <- factor(jun23_abundance$Trophic_role, levels = custom_labels)

junA<-ggplot(jun23_abundance, aes(x = Trophic_role, y = Abundance, fill = Treatment)) +
  geom_boxplot() +
  #facet_wrap(~ Trophic_role, scales = "free") +
  scale_fill_manual(values = c("Ark" = "darkcyan", "Control" = "wheat3")) +
  xlab("Trophic Category") +
  ylab("Number of Fish") +
  theme_light(base_size = 14)
junA

ggsave("jun_trophic_numbers.jpg", junA, width = 8, height = 6, units = "in", dpi = 1200)

```

Fish species richness

```{r}
library(dplyr)
library(ggplot2)
library(lubridate)
library(wesanderson)
library(rstatix)

fish<-read.csv("fish_paper.csv",header=T)
#set dates
fish$Date<-as.Date(fish$Date,format="%m/%d/%y")

#Fish species richness
richness<-fish %>%
  group_by(Treatment,Site,Survey,Date)%>%
  summarize(count=n_distinct(Scientific_name))
richness<-as.data.frame(richness)

#check for change with time
ark_rich<-filter(richness,Treatment=="Ark")
control_rich<-filter(richness,Treatment=="Control")
rich_time_ark<-aov(count~Date,data=ark_rich)
summary(rich_time_ark)
# Df Sum Sq Mean Sq F value Pr(>F)  
# Date         1  72.89   72.89   4.982 0.0425 *
#  Residuals   14 204.86   14.63                
#---
#  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
rich_time_control<-aov(count~Date,data=control_rich)
summary(rich_time_control)
# Df Sum Sq Mean Sq F value Pr(>F)
# Date         1   18.7   18.69   0.452  0.518
# Residuals    9  372.0   41.34


#plot richness 

custom_labels <- c("Nov 2021", "Feb 2022", "May 2022", "Aug 2022", "Dec 2022", "June 2023")
richness$Date<-as.factor(richness$Date)

Fish_rich <- richness %>%
  ggplot(aes(x = Date, y = count, color = Treatment, group = Treatment)) +
  stat_summary(
    fun.data = mean_se,  # Calculates mean and standard error
    geom = "errorbar",   # Draws error bars
    width = 0.2         # Width of error bars
  ) +
  stat_summary(
    fun.y = "mean",  # Calculates the mean
    geom = "point",  # Draws a point for the mean
    size = 3  # Adjust the size of the point
  ) +
  stat_summary(
    fun.y = "mean",  # Calculates the mean
    geom = "line",  # Draws lines between points
    size = 0.5  # Adjust the size of the line
  ) +
  theme_test(base_size = 14) +
  xlab("Monitoring Event") +
  ylab("Species Richness") +
  scale_color_manual(values = c("darkcyan", "wheat4")) +
  labs(fill = "Treatment") +
  coord_cartesian(ylim = c(0, 40)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_x_discrete(labels = custom_labels)
Fish_rich
ggsave("Fish_rich.jpg", Fish_rich, width = 8, height = 6, units = "in", dpi = 1200)

```

Fish species evenness

```{r}
library(dplyr)
library(ggplot2)
library(lubridate)
library(wesanderson)
library(rstatix)
library(vegan)
library(tidyr)

fish<-read.csv("fish_paper.csv",header=T)
#set dates
fish$Date<-as.Date(fish$Date,format="%m/%d/%y")

#Fish species richness
richness<-fish %>%
  group_by(Treatment,Site,Survey,Date)%>%
  summarize(count=n_distinct(Scientific_name))
richness<-as.data.frame(richness)

#Fish species evenness
#first calculate shannon diversity index
shan1<-fish %>%
  group_by(Treatment,Site,Survey,Date,Scientific_name)%>%
  summarize(numfish=sum(Abundance))
shan1<-as.data.frame(shan1)
shan1wide<-pivot_wider(shan1,names_from=Scientific_name,values_from=numfish,values_fill=0)
shan1wide <- as.data.frame(shan1wide)
shan2tot<-fish %>%
  group_by(Treatment,Site,Survey,Date) %>%
  summarize(totalfish=sum(Abundance))
pShannon<-shan1wide[,5:72]/shan2tot$totalfish
plnpShannon<-pShannon*(log(pShannon))
  shan3<-cbind(shan1wide[,1:4],plnpShannon)
cols_to_sum<-c(names(shan3[5:72]))
shan3$Shannon <- (rowSums(shan3[cols_to_sum],na.rm=TRUE))*(-1)
Shan_summary<-cbind(shan3[,1:4],shan3$Shannon)
Shan_even<-cbind(Shan_summary,richness) #visually check Shan_even and see that richness & Shannon diversity values are associated with same Treatment, Site, Survey, Date entries
#remove duplicate columns
Shan_even <- Shan_even[, !duplicated(t(Shan_even))]
Shan_even$Evenness <- shan3$Shannon/log(Shan_even$count)

#check for change in evenness with time
ark_even<-filter(Shan_even,Treatment=="Ark")
control_even<-filter(Shan_even,Treatment=="Control")
even_time_ark<-aov(Evenness~Date,data=ark_even)
summary(even_time_ark)
# Df Sum Sq Mean Sq F value Pr(>F)  
#Date         1 0.3844  0.3844   14.42 0.00254 **
#Residuals   12 0.3198  0.0266              
even_time_control<-aov(Evenness~Date,data=control_even)
summary(even_time_control)
# Df Sum Sq Mean Sq F value Pr(>F)
# Date         1 0.000148 0.0001479   0.254  0.626
#Residuals    9 0.005239 0.0005821 


#plot evenness

Shan_even$Date<-as.factor(Shan_even$Date)

Fish_even <- Shan_even %>%
  ggplot(aes(x = Date, y = Evenness, color = Treatment, group = Treatment)) +
  stat_summary(
    fun.data = mean_se,  # Calculates mean and standard error
    geom = "errorbar",   # Draws error bars
    width = 0.2         # Width of error bars
  ) +
  stat_summary(
    fun.y = "mean",  # Calculates the mean
    geom = "point",  # Draws a point for the mean
    size = 3  # Adjust the size of the point
  ) +
  stat_summary(
    fun.y = "mean",  # Calculates the mean
    geom = "line",  # Draws lines between points
    size = 0.5  # Adjust the size of the line
  ) +
  theme_test(base_size = 14) +
  xlab("Monitoring Event") +
  ylab("Species Evenness") +
  scale_color_manual(values = c("darkcyan", "wheat4")) +
  labs(fill = "Treatment") +
  coord_cartesian(ylim = c(0, 1)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_x_discrete(labels = custom_labels)
Fish_even
ggsave("Fish_even.jpg", Fish_even, width = 8, height = 6, units = "in", dpi = 1200)

```

Fish biomass and abundance before/after ARMS added

```{r}

library(dplyr)
library(ggplot2)
library(lubridate)
library(wesanderson)
library(rstatix)

fish<-read.csv("fish_paper.csv",header=T)
#set formats
fish$Date<-as.Date(fish$Date,format="%m/%d/%y")
fish$Site<-as.numeric(fish$Site)
fish$Survey<-as.numeric(fish$Survey)
fish$Tropic_role<-as.factor(fish$Trophic_role)

#abundance
abundance_trophic<-fish %>%
  group_by(Treatment,Date,Site,Survey,Trophic_role)%>%
  summarize(sum(Abundance))
abundance_trophic<-as.data.frame(abundance_trophic)
colnames(abundance_trophic)<-c("Treatment","Date","Site","Survey","Trophic_role","Abundance")

#biomass
biomass_trophic_fish<-fish %>%
  group_by(Treatment,Date,Site,Survey,Trophic_role)%>%
  summarize(sum(Biomass))
biomass_trophic_fish<-as.data.frame(biomass_trophic_fish)
colnames(biomass_trophic_fish)<-c("Treatment","Date","Site","Survey","Trophic_role","Biomass")

#first prep data - biomass

#before ARMS - filter by time
before_biomass=filter(biomass_trophic_fish,Date %in% c("2021-11-21","2022-02-27","2022-05-08"))
#add zero Piscivore observation for Control site 2 survey 1 May 2022
cont1<-data.frame("Control","2022-05-08",2,1,"Piscivore",0)
colnames(cont1)<-c("Treatment","Date","Site","Survey","Trophic_role","Biomass")
before_biomass<-rbind(before_biomass,cont1)
#filter by trophic role
before_biomass_Pisc<-filter(before_biomass,Trophic_role=="Piscivore")
before_biomass_Plank<-filter(before_biomass,Trophic_role=="Planktivore")
before_biomass_Prim<-filter(before_biomass,Trophic_role=="Primary")
before_biomass_Sec<-filter(before_biomass,Trophic_role=="Secondary")

#after ARMS - filter by time
after_biomass=filter(biomass_trophic_fish,Date %in% c("2022-08-28","2022-12-11","2023-06-04"))
#add zero Piscivore observations for Control (site 1 survey 1 Aug 2022, site 1 survey 2 June 2023)
cont2<-data.frame("Control","2022-08-28",1,1,"Piscivore",0)
colnames(cont2)<-c("Treatment","Date","Site","Survey","Trophic_role","Biomass")
cont3<-data.frame("Control","2023-06-04",1,2,"Piscivore",0)
colnames(cont3)<-c("Treatment","Date","Site","Survey","Trophic_role","Biomass")
#add zero Planktivore observations for Arks June 2023
a11j<-data.frame("Ark","2023-06-04",1,1,"Planktivore",0)
colnames(a11j)<-c("Treatment","Date","Site","Survey","Trophic_role","Biomass")
a21j<-data.frame("Ark","2023-06-04",2,1,"Planktivore",0)
colnames(a21j)<-c("Treatment","Date","Site","Survey","Trophic_role","Biomass")
a22j<-data.frame("Ark","2023-06-04",2,2,"Planktivore",0)
colnames(a22j)<-c("Treatment","Date","Site","Survey","Trophic_role","Biomass")
after_biomass<-rbind(after_biomass,cont2,cont3,a11j,a21j,a22j)
#filter by trophic roles
after_biomass_Pisc<-filter(after_biomass,Trophic_role=="Piscivore")
after_biomass_Plank<-filter(after_biomass,Trophic_role=="Planktivore")
after_biomass_Prim<-filter(after_biomass,Trophic_role=="Primary")
after_biomass_Sec<-filter(after_biomass,Trophic_role=="Secondary")

#test normality - biomass
  #before
shapiro.test(before_biomass_Pisc$Biomass[before_biomass_Pisc$Treatment=="Ark"]) #normal
shapiro.test(before_biomass_Pisc$Biomass[before_biomass_Pisc$Treatment=="Control"]) #normal
shapiro.test(before_biomass_Plank$Biomass[before_biomass_Plank$Treatment=="Ark"]) #normal
shapiro.test(before_biomass_Plank$Biomass[before_biomass_Plank$Treatment=="Control"]) #not normal
#try square-root transformation to see if becomes normal
before_biomass_Plank$Biomass_sq<-sqrt(before_biomass_Plank$Biomass)
shapiro.test(before_biomass_Plank$Biomass_sq[before_biomass_Plank$Treatment=="Control"]) #still not normal
shapiro.test(before_biomass_Prim$Biomass[before_biomass_Prim$Treatment=="Ark"]) #normal
shapiro.test(before_biomass_Prim$Biomass[before_biomass_Prim$Treatment=="Control"]) #normal
shapiro.test(before_biomass_Sec$Biomass[before_biomass_Sec$Treatment=="Ark"]) #normal
shapiro.test(before_biomass_Sec$Biomass[before_biomass_Sec$Treatment=="Control"]) #normal
  #after
shapiro.test(after_biomass_Pisc$Biomass[after_biomass_Pisc$Treatment=="Ark"]) #normal
shapiro.test(after_biomass_Pisc$Biomass[after_biomass_Pisc$Treatment=="Control"]) #not normal
#try square-root transformation to see if becomes normal
after_biomass_Pisc$Biomass_sq<-sqrt(after_biomass_Pisc$Biomass) 
shapiro.test(after_biomass_Pisc$Biomass_sq[after_biomass_Pisc$Treatment=="Control"]) #still not normal
shapiro.test(after_biomass_Plank$Biomass[after_biomass_Plank$Treatment=="Ark"]) #not normal
#try square-root transformation to see if becomes normal
after_biomass_Plank$Biomass_sq<-sqrt(after_biomass_Plank$Biomass)
shapiro.test(after_biomass_Plank$Biomass_sq[after_biomass_Plank$Treatment=="Ark"]) #normal
shapiro.test(after_biomass_Plank$Biomass[after_biomass_Plank$Treatment=="Control"]) #normal
shapiro.test(after_biomass_Prim$Biomass[after_biomass_Prim$Treatment=="Ark"]) #normal
shapiro.test(after_biomass_Prim$Biomass[after_biomass_Prim$Treatment=="Control"]) #normal
shapiro.test(after_biomass_Sec$Biomass[after_biomass_Sec$Treatment=="Ark"]) #not normal
#try square-root transformation to see if becomes normal
after_biomass_Sec$Biomass_sq<-sqrt(after_biomass_Sec$Biomass)
shapiro.test(after_biomass_Sec$Biomass_sq[after_biomass_Sec$Treatment=="Ark"]) #normal
shapiro.test(after_biomass_Sec$Biomass[after_biomass_Sec$Treatment=="Control"]) #normal

#test for significant differences by time within treatment - biomass
#Arks 
 #Piscivores
t.test(before_biomass_Pisc$Biomass[before_biomass_Pisc$Treatment=="Ark"],after_biomass_Pisc$Biomass[after_biomass_Pisc$Treatment=="Ark"])
# t = -4.9293, df = 13.572, p-value = 0.0002428
 #Planktivores
wilcox.test(before_biomass_Plank$Biomass_sq[before_biomass_Plank$Treatment=="Ark"],after_biomass_Plank$Biomass_sq[after_biomass_Plank$Treatment=="Ark"])
# W = 22, p-value = 0.8313
 #Primary consumers
t.test(before_biomass_Prim$Biomass[before_biomass_Prim$Treatment=="Ark"],after_biomass_Prim$Biomass[after_biomass_Prim$Treatment=="Ark"])
#t = -4.4781, df = 9.9599, p-value = 0.001194
 #Secondary consumers
#need to sqrt transform before before_biomass_Sec to compare to sqrt-transformed after data
before_biomass_Sec$Biomass_sq<-sqrt(before_biomass_Sec$Biomass)
t.test(before_biomass_Sec$Biomass_sq[before_biomass_Sec$Treatment=="Ark"],after_biomass_Sec$Biomass_sq[after_biomass_Sec$Treatment=="Ark"])
#t = -3.1128, df = 9.9938, p-value = 0.01102

#Control
 #Piscivores
wilcox.test(before_biomass_Pisc$Biomass[before_biomass_Pisc$Treatment=="Control"],after_biomass_Pisc$Biomass[before_biomass_Pisc$Treatment=="Control"])
#W = 5, p-value = 0.358
 #Planktivores
wilcox.test(before_biomass_Plank$Biomass[before_biomass_Plank$Treatment=="Control"],after_biomass_Plank$Biomass[before_biomass_Plank$Treatment=="Control"])
#W = 9, p-value = 1
 #Primary consumers
t.test(before_biomass_Prim$Biomass[before_biomass_Prim$Treatment=="Control"],after_biomass_Prim$Biomass[after_biomass_Prim$Treatment=="Control"])
#t = 0.65995, df = 2.462, p-value = 0.5657
 #Secondary consumers
t.test(before_biomass_Sec$Biomass[before_biomass_Sec$Treatment=="Control"],after_biomass_Sec$Biomass[after_biomass_Sec$Treatment=="Control"])
#t = 1.4253, df = 2.7278, p-value = 0.2578

#first prep data - abundance

#before ARMS - filter by time
before_Abundance=filter(abundance_trophic,Date %in% c("2021-11-21","2022-02-27","2022-05-08"))
#add zero Piscivore observation for Control site 2 survey 1 May 2022
cont1<-data.frame("Control","2022-05-08",2,1,"Piscivore",0)
colnames(cont1)<-c("Treatment","Date","Site","Survey","Trophic_role","Abundance")
before_Abundance<-rbind(before_Abundance,cont1)
#filter by trophic role
before_Abundance_Pisc<-filter(before_Abundance,Trophic_role=="Piscivore")
before_Abundance_Plank<-filter(before_Abundance,Trophic_role=="Planktivore")
before_Abundance_Prim<-filter(before_Abundance,Trophic_role=="Primary")
before_Abundance_Sec<-filter(before_Abundance,Trophic_role=="Secondary")

#after ARMS - filter by time
after_Abundance=filter(abundance_trophic,Date %in% c("2022-08-28","2022-12-11","2023-06-04"))
#add zero Piscivore observations for Control (site 1 survey 1 Aug 2022, site 1 survey 2 June 2023)
cont2<-data.frame("Control","2022-08-28",1,1,"Piscivore",0)
colnames(cont2)<-c("Treatment","Date","Site","Survey","Trophic_role","Abundance")
cont3<-data.frame("Control","2023-06-04",1,2,"Piscivore",0)
colnames(cont3)<-c("Treatment","Date","Site","Survey","Trophic_role","Abundance")
#add zero Planktivore observations for Arks June 2023
a11j<-data.frame("Ark","2023-06-04",1,1,"Planktivore",0)
colnames(a11j)<-c("Treatment","Date","Site","Survey","Trophic_role","Abundance")
a21j<-data.frame("Ark","2023-06-04",2,1,"Planktivore",0)
colnames(a21j)<-c("Treatment","Date","Site","Survey","Trophic_role","Abundance")
a22j<-data.frame("Ark","2023-06-04",2,2,"Planktivore",0)
colnames(a22j)<-c("Treatment","Date","Site","Survey","Trophic_role","Abundance")
after_Abundance<-rbind(after_Abundance,cont2,cont3,a11j,a21j,a22j)
#filter by trophic roles
after_Abundance_Pisc<-filter(after_Abundance,Trophic_role=="Piscivore")
after_Abundance_Plank<-filter(after_Abundance,Trophic_role=="Planktivore")
after_Abundance_Prim<-filter(after_Abundance,Trophic_role=="Primary")
after_Abundance_Sec<-filter(after_Abundance,Trophic_role=="Secondary")

#test normality - Abundance
#before
shapiro.test(before_Abundance_Pisc$Abundance[before_Abundance_Pisc$Treatment=="Ark"]) #not normal
#try square-root transformation to see if becomes normal
before_Abundance_Pisc$Abundance_sq<-sqrt(before_Abundance_Pisc$Abundance)
shapiro.test(before_Abundance_Pisc$Abundance_sq[before_Abundance_Pisc$Treatment=="Ark"]) #normal
shapiro.test(before_Abundance_Pisc$Abundance[before_Abundance_Pisc$Treatment=="Control"]) #normal
shapiro.test(before_Abundance_Plank$Abundance[before_Abundance_Plank$Treatment=="Ark"]) #normal
shapiro.test(before_Abundance_Plank$Abundance[before_Abundance_Plank$Treatment=="Control"]) #normal
shapiro.test(before_Abundance_Prim$Abundance[before_Abundance_Prim$Treatment=="Ark"]) #normal
shapiro.test(before_Abundance_Prim$Abundance[before_Abundance_Prim$Treatment=="Control"]) #normal
shapiro.test(before_Abundance_Sec$Abundance[before_Abundance_Sec$Treatment=="Ark"]) #normal
shapiro.test(before_Abundance_Sec$Abundance[before_Abundance_Sec$Treatment=="Control"]) #normal
#after
shapiro.test(after_Abundance_Pisc$Abundance[after_Abundance_Pisc$Treatment=="Ark"]) #normal
shapiro.test(after_Abundance_Pisc$Abundance[after_Abundance_Pisc$Treatment=="Control"]) #normal
shapiro.test(after_Abundance_Plank$Abundance[after_Abundance_Plank$Treatment=="Ark"]) #normal
shapiro.test(after_Abundance_Plank$Abundance[after_Abundance_Plank$Treatment=="Control"]) #normal
shapiro.test(after_Abundance_Prim$Abundance[after_Abundance_Prim$Treatment=="Ark"]) #normal
shapiro.test(after_Abundance_Prim$Abundance[after_Abundance_Prim$Treatment=="Control"]) #normal
shapiro.test(after_Abundance_Sec$Abundance[after_Abundance_Sec$Treatment=="Ark"]) #normal
shapiro.test(after_Abundance_Sec$Abundance[after_Abundance_Sec$Treatment=="Control"]) #normal

#test for significant differences by time within treatment - abundance
#Arks 
  #Piscivores
#sqrt transform after to compare to before Pisc
after_Abundance_Pisc$Abundance_sq<-sqrt(after_Abundance_Pisc$Abundance)
t.test(before_Abundance_Pisc$Abundance_sq[before_Abundance_Pisc$Treatment=="Ark"],after_Abundance_Pisc$Abundance_sq[before_Abundance_Pisc$Treatment=="Ark"])
# t = -3.3887, df = 13.552, p-value = 0.004596
  #Planktivores
t.test(before_Abundance_Plank$Abundance[before_Abundance_Plank$Treatment=="Ark"],after_Abundance_Plank$Abundance[before_Abundance_Pisc$Treatment=="Ark"])
# t = -3.1311, df = 11.261, p-value = 0.009308
  #Primary consumers
t.test(before_Abundance_Prim$Abundance[before_Abundance_Prim$Treatment=="Ark"],after_Abundance_Prim$Abundance[after_Abundance_Prim$Treatment=="Ark"])
#t = -4.1243, df = 10.965, p-value = 0.0017
  #Secondary consumers
t.test(before_Abundance_Sec$Abundance[before_Abundance_Sec$Treatment=="Ark"],after_Abundance_Sec$Abundance[after_Abundance_Sec$Treatment=="Ark"])
#t = -2.0161, df = 11.76, p-value = 0.06722

#Control
  #Piscivores
t.test(before_Abundance_Pisc$Abundance[before_Abundance_Pisc$Treatment=="Control"],after_Abundance_Pisc$Abundance[before_Abundance_Pisc$Treatment=="Control"])
#t = -1.9602, df = 5.2255, p-value = 0.1048
  #Planktivores
t.test(before_Abundance_Plank$Abundance[before_Abundance_Plank$Treatment=="Control"],after_Abundance_Plank$Abundance[before_Abundance_Plank$Treatment=="Control"])
#t = 0.40869, df = 5.0864, p-value = 0.6994
  #Primary consumers
t.test(before_Abundance_Prim$Abundance[before_Abundance_Prim$Treatment=="Control"],after_Abundance_Prim$Abundance[after_Abundance_Prim$Treatment=="Control"])
#t = -0.54575, df = 5.277, p-value = 0.6075
  #Secondary consumers
t.test(before_Abundance_Sec$Abundance[before_Abundance_Sec$Treatment=="Control"],after_Abundance_Sec$Abundance[after_Abundance_Sec$Treatment=="Control"])
#t = 0.9522, df = 4.8549, p-value = 0.386


```