#####################################################
# Marco FW Gauger
# 2021/02/23
#
# This is the script to obtain fig 2 of the manuscript submitted to PeerJ Gauger et al 2021 Diel influences on bottlenose dolphin acoustic detection in a coastal lagoon in the southwestern Gulf of California
#####################################################

ifelse("ggplot2" %in% rownames(installed.packages()) == FALSE, 
       install.packages(c("ggplot2","colorspace"), dependencies = T),
       ("ggplot2 is already installed"))
ifelse("ggpubr" %in% rownames(installed.packages()) == FALSE, 
       install.packages("ggpubr", dependencies = T),
       ("ggpubr is already installed"))
ifelse("PMCMRplus" %in% rownames(installed.packages()) == FALSE, 
       install.packages("PMCMRplus", dependencies = T),
       ("PMCMRplus is already installed"))
ifelse("PMCMR" %in% rownames(installed.packages()) == FALSE, 
       install.packages("PMCMR", dependencies = T),
       ("PMCMR is already installed"))
ifelse("WRS2" %in% rownames(installed.packages()) == FALSE, 
       install.packages("WRS2", dependencies = T),
       ("WRS2 is already installed"))
ifelse("nortest" %in% rownames(installed.packages()) == FALSE, 
       install.packages("nortest", dependencies = T),
       ("nortest is already installed"))
ifelse("car" %in% rownames(installed.packages()) == FALSE, 
       install.packages(c("car","utf8"), dependencies = T),
       ("car is already installed"))
ifelse("multcompView" %in% rownames(installed.packages()) == FALSE, 
       install.packages("multcompView", dependencies = T),
       ("multcompView is already installed"))
ifelse("cowplot" %in% rownames(installed.packages()) == FALSE, 
       install.packages("cowplot", dependencies = T),
       ("cowplot is already installed"))
ifelse("cowplot" %in% rownames(installed.packages()) == FALSE, 
       install.packages("cowplot", dependencies = T),
       ("cowplot is already installed"))
ifelse("tibble" %in% rownames(installed.packages()) == FALSE, 
       install.packages("tibble", dependencies = T),
       ("tibble is already installed"))

require(rstatix)
require(ggpubr)
require(PMCMRplus)
require(PMCMR)
require(WRS2)
require(nortest)
require(car)
require(multcompView)
require(ggplot2)
require(cowplot)

analysed_data_hour <- read.csv("Supplemental Data S1.csv", sep=";")
str(analysed_data_hour)

analysed_data_hour$fcampaign <- as.factor(analysed_data_hour$month_campaign)
str(analysed_data_hour$cluster_hcpc)
nortest::ad.test(analysed_data_hour$dp10m_ELAP1)
car::leveneTest(analysed_data_hour$dp10m_ELAP1~ as.factor(analysed_data_hour$month))
car::leveneTest(dp10m_ELAP1~ fcampaign, analysed_data_hour)
car::leveneTest(dp10m_ELAP1~ as.factor(cluster_hcpc), analysed_data_hour)


unique(analysed_data_hour$month_campaign)
table(analysed_data_hour$date, analysed_data_hour$month)
analysed_data_hour$four_season = rep("spring", nrow(analysed_data_hour))
analysed_data_hour$four_season[analysed_data_hour$month %in% c(1,2)] <- "winter"
analysed_data_hour$four_season[analysed_data_hour$month %in% c(7,8,9)] <- "summer"
analysed_data_hour$four_season[analysed_data_hour$month %in% c(10,11)] <- "autumn"
analysed_data_hour$four_season[analysed_data_hour$dates %in% c("2019-03-28", "2019-03-29")] <- "spring"
analysed_data_hour$four_season[analysed_data_hour$dates %in% c("2017-06-21","2017-06-22","2017-06-23","2018-06-24", "2018-06-29")] <- "summer"
analysed_data_hour$four_season[analysed_data_hour$dates %in% c("2017-09-25", "2017-09-26","2017-09-27")] <- "autumn"
analysed_data_hour$four_season <- ordered(analysed_data_hour$four_season, levels = c("spring", "summer", "autumn","winter"))

table(analysed_data_hour$month,analysed_data_hour$four_season)

boxplot(analysed_data_hour$dp10m_ELAP1~analysed_data_hour$four_season)
kruskal.test(analysed_data_hour$dp10m_ELAP1~analysed_data_hour$four_season)
PMCMRplus::kwAllPairsNemenyiTest(analysed_data_hour$dp10m_ELAP1~analysed_data_hour$four_season, dist = "Chisq")
analysed_data_hour$two_season = rep("June-Nov", nrow(analysed_data_hour))
analysed_data_hour$two_season[analysed_data_hour$month %in% c("01","02","03","04","05")] <- "Jan-May"
table(analysed_data_hour$cluster_hcpc)
table(analysed_data_hour$month,analysed_data_hour$cluster_hcpc)

analysed_data_hour$month <- formatC(analysed_data_hour$month, flag="0", format="fg", digits=1)

out <- kwAllPairsNemenyiTest(dp10m_ELAP1~as.factor(month),
                             analysed_data_hour, dist=c("Chisquare"))
letters_out <- PMCMR::get.pvalues(out)
letters_out<- multcompView::multcompLetters(letters_out, threshold=0.05)

myletters_df <- data.frame(month_campaign=names(letters_out$Letters),letter = letters_out$Letters )
head(myletters_df)
l = data.frame(new = letters[1:length(unique(myletters_df$letter))], old = unique(myletters_df$letter))
myletters_df$unique_letter <- rep(0,nrow(myletters_df))
for(k in 1:nrow(l)){
myletters_df$unique_letter[myletters_df$letter %in% l$old[k]] <- k
}
myletters_df$new_letters <- LETTERS[myletters_df$unique_letter]
myletters_df$ty = 1
myletters_df$ty[myletters_df$unique_letter > 1] <- 2
myletters_df$month <- as.factor(substr(myletters_df$month_campaign,1,2))
myletters_df_month <- myletters_df

(boxpplot_month <- ggplot(analysed_data_hour, aes(y=dp10m_ELAP1, x=month, color= month))+
geom_boxplot(outlier.alpha = 0, alpha=0.25)+geom_jitter(width=0.25)+
#stat_summary(fun=mean, colour="black", geom="point", shape=18, size=3) +
theme_classic()+
#ylim(0,9)+
theme(legend.position="none")+
theme(axis.text.x = element_text(angle=40, hjust=1, vjust=1))+
ylab(expression("dp10m h"^-1))+ xlab("month"))

sigs =
data.frame(first=c("01","03","03","05","05","06","06","06","07","07","07"),
second=c("07","09","10","06","07","09","10","11","09","10","11"),
add=(c(1:11)/11)*2+6.5,
color=c(1,1,2,1,2,2,2,1,2,2,2))
(boxpplot_month <- boxpplot_month+
annotate("segment", x = sigs$first, xend = sigs$second,
y = sigs$add, yend = sigs$add,
colour = as.factor(sigs$color))+
annotate("segment", x = sigs$first, xend = sigs$first,
y = sigs$add, yend = sigs$add-.07,
colour = as.factor(sigs$color))+
annotate("segment", x = sigs$second, xend = sigs$second,
y = sigs$add, yend = sigs$add-.07,
colour = as.factor(sigs$color)))
levels(analysed_data_hour$month_campaign)
analysed_data_hour$month_campaign <- ordered(analysed_data_hour$month_campaign, levels = sort(unique(analysed_data_hour$month_campaign)))





out <- PMCMRplus::kwAllPairsNemenyiTest(dp10m_ELAP1~as.factor(month_campaign),
analysed_data_hour, dist=c("Chisquare"))
letters_out <- PMCMR::get.pvalues(out)
sigs_letters_out <- cbind(letters_out[letters_out < 0.05])
(sigs_letters_out <- 
  data.frame(first      = substr(rownames(sigs_letters_out),1,5),
             second     = substr(rownames(sigs_letters_out),7,12),
             p.value    = sigs_letters_out))
sigs_letters_out$col <- rep(1, nrow(sigs_letters_out))
sigs_letters_out$col[sigs_letters_out$p.value < 0.01] <- 2
sigs_letters_out$add <- 
  ((1:nrow(sigs_letters_out))/nrow(sigs_letters_out))*2+6.5


letters_out<- multcompView::multcompLetters(letters_out, threshold=0.05)
out
myletters_df <- data.frame(month_campaign=names(letters_out$Letters),letter = letters_out$Letters )
head(myletters_df)
l = data.frame(new = letters[1:length(unique(myletters_df$letter))], old = unique(myletters_df$letter))
myletters_df$unique_letter <- rep(0,nrow(myletters_df))
for(k in 1:nrow(l)){
  myletters_df$unique_letter[myletters_df$letter %in% l$old[k]] <- k
}
myletters_df$new_letters <- LETTERS[myletters_df$unique_letter]
myletters_df$ty = 1
myletters_df$ty[myletters_df$unique_letter > 1] <- 2
myletters_df$month <- as.factor(substr(myletters_df$month_campaign,1,2))
myletters_df

analysed_data_hour$MONTH <- as.factor(analysed_data_hour$month)
analysed_data_hour$MONTH[analysed_data_hour$month_campaign == "10_13"] <- "10"
(boxpplot_campaign0 <- ggplot(analysed_data_hour, aes(y=dp10m_ELAP1, x=month_campaign, color= MONTH))+
geom_boxplot(outlier.alpha = 0, alpha=0.25)+geom_jitter(width=0.25)+
#stat_summary(fun=mean, colour="black", geom="point",shape=18, size=3) +
theme_classic()+
#ylim(0,9.0)+
theme(legend.position="none")+
theme(axis.text.x = element_text(angle=40, hjust=1, vjust=1))+
ylab(expression("dp10m h"^-1))+ xlab(""))
boxpplot_campaign <- 
  boxpplot_campaign0 +
  annotate("segment", x = sigs_letters_out$first,
           xend = sigs_letters_out$second,
           y = sigs_letters_out$add,
           yend = sigs_letters_out$add,
           colour = as.factor(sigs_letters_out$col))+
  annotate("segment", x = sigs_letters_out$first, 
           xend = sigs_letters_out$first,
           y = sigs_letters_out$add, 
           yend = sigs_letters_out$add-.05,
           colour = as.factor(sigs_letters_out$col))+
  annotate("segment", x = sigs_letters_out$second, 
           xend = sigs_letters_out$second,
           y = sigs_letters_out$add, 
           yend = sigs_letters_out$add-.05,
           colour = as.factor(sigs_letters_out$col))



#x11()
(boxpplot_campaign_month <- cowplot::plot_grid(boxpplot_campaign+ xlab("deployment"),boxpplot_month+ xlab("month"), nrow=2, ncol=1, axis="l", align="h", labels = c("A","B")))
table(myletters_df$month)




(j <- aggregate(analysed_data_hour$SST, list(analysed_data_hour$month_campaign),mean))
names(j) <- c("month_campaign","SST")
(j$new <- paste0("St.",
                 aggregate(analysed_data_hour$Station, list(analysed_data_hour$month_campaign),dplyr::first)$x,
                 "_",
                 aggregate(analysed_data_hour$date, list(analysed_data_hour$month_campaign),dplyr::first)$x))
j$new[j$month_campaign == "07_04"] <- "St.4_13.07.2018"

stat_sum_df <- function(fun, geom="crossbar", ...) {
  stat_summary(fun.data = fun, colour = "red", geom = geom, width = 0.2, ...)
}
(p_campaign <- boxpplot_campaign +
  #  geom_smooth(aes(SST/4~month_campaign), method="loess", col="blue")
    stat_summary(aes(y = SST-20), fun = "mean", colour = "black", bg="white", size = 3, geom = "point", shape=21)+stat_summary(aes(y = SST-20), fun = "mean", colour = "black", size = 1.5, geom = "point", shape='')+  scale_y_continuous(expression("dp10m h"^-1), sec.axis = sec_axis(~ .+20, name = "SST [°C]"))+ xlab("deployment"))
p_campaign + scale_x_discrete(breaks=j$month_campaign,
                   labels=j$new)

ggsave("Figure 3b.png",p_campaign + scale_x_discrete(breaks=j$month_campaign,
                                                      labels=j$new), 
       dpi=300, units="cm", width=15, height=10)