###### One-way ##########

Data <- read_csv("~/SBC_stat.csv")

library(lmPerm)
mod <- aovp(ratio ~ Site, data=Data) # one-way
summary(mod)

Data$Site = factor(Data$Site, ordered=FALSE, levels=c("Na16", "Na17", "Ha02", "Ha16", "Ei", "M.biserialis")) # Order factors, otherwise R will alphabetize them
boxplot(ratio ~ Site,data = Data,ylab ="Length to width ratio",xlab ="Site") # "Length to width ratio"

# Data = Data[-(151:188),] # for without M. biserialis option

#Permutation test
library(coin)
independence_test(ratio ~ Site,data = Data) # one-way

Data$Site = factor(Data$Site,levels = c("M.biserialis", "Ei", "Ha02", "Ha16", "Na16", "Na17", "Na18")) ### Order groups by median
# Data$Site = factor(Data$Site) ### Order groups by median

#Table output with pairwisePermutationTest
library(FSA)
headtail(Data)
library(rcompanion)
PT = pairwisePermutationTest(ratio ~ Site, data=Data, method="fdr")
PT
library(multcompView)
cldList(comparison = PT$Comparison, p.value = PT$p.adjust, threshold  = 0.05)
PM = pairwisePermutationMatrix(ratio ~ Site,data=Data, method="fdr")
PM


# The differeance between adjusted and unadjusted: the adjusted relates to the familywise comparisons


########### 2-way  ########

Data2 <- read_csv("~/SBC_stat_2way.csv")

mod2 <- aovp(val~ Site + measurement, data=Data2) # two-way
mod2 <- aovp(val~ Site + Error(measurement), data=Data2) # Nested
summary(mod2)


Data2$Site = factor(Data2$Site, ordered=FALSE, levels=unique(Data2$Site)) # Order factors, otherwise R will alphabetize them
Data2$measurement = factor(Data2$measurement, ordered=FALSE, levels=unique(Data2$measurement)) # Order factors, otherwise R will alphabetize them
boxplot(val ~ measurement+Site,data = Data2,ylab ="value",xlab ="Measurement + Site")
# boxplot(val ~ Site,data = Data,ylab ="value",xlab ="Site")
# boxplot(val ~ measurement,data = Data,ylab ="value",xlab ="Measurement")

library(rcompanion)
PT = pairwisePermutationTest(val ~ Site_and_measurement,data=Data2,method="fdr") # two-way
PT = pairwisePermutationTest(val ~ Site + Error(measurement),data=Data2,method="fdr") # Nested
PT
library(rcompanion)
library(multcompView)
cldList(comparison = PT$Comparison, p.value = PT$p.adjust, threshold  = 0.05)
PM = pairwisePermutationMatrix(val ~ Site_and_measurement,data=Data2,method="fdr")
PM

# Just one factor:
measure = "ratio"  # long, short, or ratio
fact = which(Data2$measurement==measure) # long, short, or ratio
mod1 <- aovp(val~Site, data=Data2[fact,]) # one-way
summary(mod1)

Data$Site = factor(Data2$Site, ordered=FALSE, levels=unique(Data2[fact,]$Site)) # Order factors, otherwise R will alphabetize them
boxplot(val ~ Site,data = Data2[fact,],ylab =measure,xlab ="Site")
PT = pairwisePermutationTest(val ~ Site, data=Data2[fact,], method="fdr")
PT

## 2-way plots ##

library(ggplot2)

Site = Data2$Site
Measurement = Data2$measurement
val = Data2$val

data = data.frame(Site, Measurement, val)

p = ggplot(data, aes(x=Site, y=val, fill=Measurement)) +
  geom_boxplot() +
  ggtitle('Measurement comparisons among sites') + # for the main title
  ylab('Normalized measurement') +
  xlab('Site')
p + scale_fill_manual(values = c("#999999", "#E69F00", "#56B4E9"))

##### warts #####

library(readr)

warts <- read_csv("~/warts.csv")

mod <- aovp(Density ~ Site, data=warts) # one-way
summary(mod)

PT = pairwisePermutationTest(Density ~ Site, data=warts, method="fdr")
PT

warts$Site = factor(warts$Site, ordered=FALSE, levels=c("Na16", "Na17", "Ha02", "Ha16", "Ei")) # levels=unique(warts$Site)) # Order factors, otherwise R will alphabetize them
boxplot(Density ~ Site,data = warts,ylab ="Warts density (# um-1)",xlab ="Site") # 


##### double box plot ######

source('~boxplotdou.R') # files were taken from the WWW: https://tomizonor.wordpress.com/2013/03/15/double-box-plot/
source('~/test.has.overlap.R') # files were taken from the WWW
# Site = factor(Site, ordered=FALSE, levels=unique(Site)) # Order factors, otherwise R will alphabetize them

long  = which(Data2$measurement=="long")
short = which(Data2$measurement=="short")

# all data
x = data.frame(Data2$Site[long], Data2$true_val[long])
y = data.frame(Data2$Site[short], Data2$true_val[short])

boxplotdou(x, y)

# library(wesanderson)
ratio = which(Data2$measurement=='ratio')
p = ggplot(Data2[ratio,], aes(x=year, y=true_val, fill=Site)) +
  geom_boxplot() +
  ggtitle('Length to width comparisons among sites') + # for the main title
  ylab('Normalized measurement') +
  xlab('Site')
p + scale_fill_brewer(palette="Dark2") # http://www.sthda.com/english/wiki/ggplot2-colors-how-to-change-colors-automatically-and-manually