setwd ("~/Desktop/big_data")
library(scales); 
library(zoo)
library(FSA)
reefcols <- c("#8dd3c7", "#bebada", "#d9d9d9")



# Import sediment data
sed.all <- read.csv("Kbay_sediment_2015.csv")
head (sed.all)
str (sed.all)
# Set date to date class
sed.all$date <- as.Date(sed.all$date, format="%Y-%m-%d")
str (sed.all)

summary (sed.all$sediment_grams_per_day)
sd (sed.all$sediment_grams_per_day)

#reduce file to only the three reefs
rf44.sed <- subset(sed.all, site=="reef44")
rf25.sed <- subset(sed.all, site=="reef25")
HIMB.sed <- subset(sed.all, site=="HIMB")

summary (rf44.sed$sediment_grams_per_day)
summary (rf25.sed$sediment_grams_per_day)
summary (HIMB.sed$sediment_grams_per_day)

# plot sediment values for each reef
#pdf (file="~/Desktop/Sed_fig1.pdf", width= 7, height= 3) 
#par(mar=c(2,3,2,1), mgp=c(2,0.5,0))
plot(sediment_grams_per_day ~ date, rf44.sed, type="n", ylab="sediment (grams/day)", xaxt="n", xlab="", 
     ylim=c(0,.2))
lines (rf44.sed$date, rf44.sed$sediment_grams_per_day, type="o", col=reefcols[1])
lines (rf25.sed$date, rf25.sed$sediment_grams_per_day, type="o", col=reefcols[2])
lines (HIMB.sed$date, HIMB.sed$sediment_grams_per_day, type="o", col=reefcols[3])
title ("Sediment Rates in Kaneohe Bay", adj=0)
axis.Date(1, at=seq(min(sed.all$date), max(sed.all$date), by="1 mon"), format="%b '%y")
legend("topright", lty=1, col=reefcols, legend=c("44","25","HIMB"), lwd=2, bty="n", inset=c(0.07,0))

#dev.off()

#Statistical Analysis

#analysis for difference in sed among all the three reefs
reefs <- rbind (rf44.sed, rf25.sed, HIMB.sed)
fit <- aov(sediment_grams_per_day ~ site, data= reefs)
summary(fit)

TukeyHSD (fit)

#check for normality
shapiro.test (reefs$sediment_grams_per_day)
# not normal P<0.001

kruskal.test (sediment_grams_per_day ~ site, data= reefs)