setwd("~/Desktop/Mo'orea/")

library(ggplot2)
library(xlsx)
library(reshape2)
data = read.xlsx("Plastic.xlsx", sheetName="FS2")
data2 = read.xlsx("Plastic.xlsx", sheetName="LS2")

##### Define Functions ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### 
#

## Summarizes data.
## Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%).
##   data: a data frame.
##   measurevar: the name of a column that contains the variable to be summariezed
##   groupvars: a vector containing names of columns that contain grouping variables
##   na.rm: a boolean that indicates whether to ignore NA's
##   conf.interval: the percent range of the confidence interval (default is 95%)
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
  library(plyr)
  
  # New version of length which can handle NA's: if na.rm==T, don't count them
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
  }
  
  # This does the summary. For each group's data frame, return a vector with
  # N, mean, and sd
  datac <- ddply(data, groupvars, .drop=.drop,
                 .fun = function(xx, col) {
                   c(N    = length2(xx[[col]], na.rm=na.rm),
                     mean = mean   (xx[[col]], na.rm=na.rm),
                     sd   = sd     (xx[[col]], na.rm=na.rm)
                   )
                 },
                 measurevar
  )
  
  # Rename the "mean" column    
  datac <- rename(datac, c("mean" = measurevar))
  
  datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
  
  # Confidence interval multiplier for standard error
  # Calculate t-statistic for confidence interval: 
  # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  
  return(datac)
}


DistPlastic = ggplot(data, aes(x=Distance_From_Pop, y=Total_Plastic)) +
  geom_point() +
  geom_smooth(method="lm") +
  theme_bw()
DistPlastic

DistWB = ggplot(data, aes(x=Distance_From_Pop, y=WB)) +
  geom_point() +
  geom_smooth(method="lm") +
  theme_bw()
DistWB

DistPercentPlastic = ggplot(data, aes(x=Distance_From_Pop, y=Percent_Plastic)) +
  geom_point() +
  geom_smooth(method="lm") +
  theme_bw()
DistPercentPlastic

fill <- "#4271AE"
line <- "#1F3552"

p10=ggplot(data, aes(x=Category,y=Total_Plastic)) +
geom_boxplot(fill = fill, colour = line,
             alpha = 0.7)
p10 <- p10 + scale_x_discrete(name = "Beach Site Type") +
  scale_y_continuous(name = "Total Plastic\n (plastic pieces collected in 5 min)",
                     breaks = seq(0, 175, 25),
                     limits=c(0, 160))

p10 <- p10 + theme_bw()
p10

p11=ggplot(data2, aes(x=Island_Side,y=Total_Plastic)) +
  geom_boxplot(fill = fill, colour = line,
               alpha = 0.7)
p11 <- p11 + scale_x_discrete(name = "Category") +
  scale_y_continuous(name = "Total Plastic",
                     breaks = seq(0, 175, 25),
                     limits=c(0, 160))
p11 <- p11 + ggtitle("Amount of Plastic on Beaches of Moorea")
p11 <- p11 + theme_bw()
p11

anova = aov(Total_Plastic~Island_Side, data=data)
anova
summary(aov(Total_Plastic~Category, data=data))


kruskal.test(Total_Plastic~Category, data=data)
kruskal.test(Total_Plastic~Island_Side, data=data)

lm.out= lm(WB~Total_Plastic, data=data)
summary(lm.out)

CMsizeP = ggplot(data2, aes(x=CM_size, y=Num_Plastic.)) +
  geom_point() +
  geom_smooth(method="lm") +
  theme_bw()
CMsizeP

lm.me= lm(Num_Plastic~CM_size., data=data2)
summary(lm.me)

p12 <- ggplot(data2, aes(x = Treatment_Time, y = Num_Plastic., fill = )) + 
  geom_boxplot(alpha=0.7) +  ylab="Amount of Plastic Particles Consumed",
xlab="Treatment Time")


data(data2)
data2$Rock_ID <- factor(data2$Rock_ID, 
                           labels = c("1", "2", "3", "4", "5", "6", "7","8","9"))

p13 <- ggplot(data2, aes(x = Rock_ID, y = Num_Plastic.)) +  
 geom_boxplot() 
p13


kruskal.test(Num_Plastic.~Treatment_Time, data=data2)


ttest = t.test(data2$Num_B,data2$Num_G)
ttest

kruskal.test(Num_Plastic.~Treatment_Time, data=data2)

kruskal.test(Num_G~Treatment_Time, data=data2)

boxplot(Num_Plastic.~Treatment_Time, data=data2,
        ylab="Amount of Plastic Particles Consumed",
        xlab="Treatment Time")

p4<-ggplot(data2,aes(x=barTT,y=barNP),fill=factor(barColor), color=factor(vs)) +  
  stat_summary(fun.y=mean,position="stack",geom="bar")+
  title(main="Plastic Eaten By Color", sub="sub-title", 
        xlab="x-axis label", ylab="y-axis label") 

p4


DF <- read.table(text="Treatment BLUE GREEN YELLOW  
84   11    0    2
108  15    8    0
", header=TRUE)

library(reshape2)
DF1 <- melt(DF, id.var="Treatment")

DF1$Treatment = as.factor(DF1$Treatment)
p <- ggplot(DF1, aes(x = Treatment, y = value, fill = variable)) +
  geom_bar(stat = "identity") +
  scale_fill_manual("variable", values=c("deepskyblue2","green4","yellow2")) +
  scale_x_discrete(breaks=c(84,108)) +
  title(main="Total Plastic Eaten by Treatment and Color", ylab="Amount of Plastic Consumed",
  xlab="Treatment Time in Hours") 
p

binom.test(8,19)
binom.test(11,15)

#######  RESULTS ##############

####### Part One: Island-Wide Field Study ######

#beach Category#
#graph for beach category# 

p10=ggplot(data, aes(x=Category,y=Total_Plastic)) +
  geom_boxplot(fill = fill, colour = line,
               alpha = 0.7)
p10 <- p10 + scale_x_discrete(name = "Beach Site Type") +
  scale_y_continuous(name = "Total Plastic \n(plastic pieces collected in 5 min)",
                     breaks = seq(0, 175, 25),
                     limits=c(0, 160))

p10 <- p10 + theme_bw()
p10
#stats for beach category#
library(PMCMR)

kruskal.test(Total_Plastic~Category, data=data)
posthoc.kruskal.nemenyi.test(g = data$Category, x = data$Total_Plastic, method="Tukey")

#Turb graph#

k45=ggplot(data, aes(x=Turb,y=Total_Plastic)) +
  geom_boxplot(fill = fill, colour = line,
               alpha = 0.7)
k45 <- k45 + scale_x_discrete(name = "Environment Type") +
  scale_y_continuous(name = "Total Plastic \n (plastic pieces collected in 5 min)",
                     breaks = seq(0, 175, 25),
                     limits=c(0, 160))
k45 <- k45 + theme_bw()
k45

#stats for Turb: significant difference shown#
kruskal.test(Total_Plastic~Turb, data=data)
posthoc.kruskal.nemenyi.test(g = data$Turb, x = data$Total_Plastic, method="Tukey")
#distance from population centers graphs#

DistPercentPlastic = ggplot(data, aes(x=Distance_From_Pop, y=Percent_Plastic)) +
  geom_point() +
  geom_smooth(method="lm") +
  xlab("Distance from Population (km)")+
  ylab("Percentage of Plastic of Pollution on Beach")+
  theme_bw()
DistPercentPlastic

DistTotalPlastic = ggplot(data, aes(x=Distance_From_Pop, y=Total_Plastic)) +
  geom_point() +
  geom_smooth(method="lm") +
  xlab("Distance from Population (km)")+
  ylab("Number of Plastic Pieces \n(pieces found in 5 min)")+
  theme_bw()
DistTotalPlastic

#distance from population centers stats#

lm.out= lm(Distance_From_Pop~Total_Plastic, data=data)
summary(lm.out)

coefficients(lm(Distance_From_Pop~Total_Plastic, data=data))

lm.out2 = lm(Distance_From_Pop~Percent_Plastic, data=data)
summary(lm.out2)

coefficients(lm(Distance_From_Pop~Percent_Plastic, data=data))

#linear regressions to show that WB~%~Total Plastic#

lm.out= lm(WB~Total_Plastic, data=data)
summary(lm.out)

lm.o=lm(WB~Percent_Plastic, data=data)
summary(lm.o)

lm.q=lm(Total_Plastic~Percent_Plastic, data=data)
summary(lm.q)

###### Part Two: Lab Study of CM #######

## prove that rocks didn't statistically vary ###

lm.me= lm(Num_Plastic~CM_size., data=data2)
summary(lm.me)


### graph of average plastic eaten per treatment ### 

boxplot(Num_Plastic.~Treatment_Time, data=data2,
        ylab="Average Numbenof Plastic Particles Consumed",
        xlab="Treatment Time")

#stats for plastic eaten/treatment: NOT STAT SIGNIFICANT#

wilcox.test(Num_Plastic.~Treatment_Time, data=data2)

#### graph of colors of plastic present: ##
##NEEDS HEADER AND CONTROL#

DF <- read.table(text="Treatment BLUE GREEN YELLOW  
84   11    0    2
108  15    8    0
", header=TRUE)

library(reshape2)
DF1 <- melt(DF, id.var="Treatment")

DF1$Treatment = as.factor(DF1$Treatment)
p <- ggplot(DF1, aes(x = Treatment, y = value, fill = variable)) +
  geom_bar(stat = "identity") +
  scale_fill_manual("variable", values=c("deepskyblue2","green4","yellow2")) +
  scale_x_discrete(breaks=c(84,108)) +
  guides(fill=guide_legend(title="Plastic Color")) +
  ylab("Number of Plastic Particles \nConsumed by Corallimorphs") +
  xlab("Treatment Time (hours)")+
  theme_bw()
p

#stats for green bead amount over treatments #
#statistically significant!# 

wilcox.test(Num_G~Treatment_Time, data=data2)