#*---------------------------------------------------------------------------------------------------
# | SCRIPT NAME: comparative_lake_and_reservoir
# | DATE CREATED: 2020-11-30
# | CREATED BY: Mathieu Auffray
# *--------------------------------------------------------------------------------------------------
# | OBJECTIVE: this script evaluates whether DAms impact is different from lake impact
#              in august 2020
# *--------------------------------------------------------------------------------------------------
# | DATA: comparative_lake_and_reservoir.csv
# *--------------------------------------------------------------------------------------------------

library(car)
library(ggplot2)
library(gridExtra)

# ---------------------------------------------------------------------------------------------------
# -----------------------------data import ----------------------------------------------------------
# ---------------------------------------------------------------------------------------------------

T <- read.table("comparative_lake_and_reservoir.csv", dec = ",", sep = ";", header = TRUE)

# data transformation

T$Fmax25 <- as.numeric(T$Fmax25)
T$Dmax25 <- as.numeric(T$Dmax25)
T$regulation <- as.character((T$regulation))

# ---------------------------------------------------------------------------------------------------
# -----------------------------check normality-------------------------------------------------------
# ---------------------------------------------------------------------------------------------------

xtabs(~regulation, data=T) # waterbodies regulated (1) = 5 / unregulated (0) = 4

shapiro.test(T$Mmo) # 0.49 
shapiro.test(T$Mmo[T$regulation=='1']) # 0.96
shapiro.test(T$Mmo[T$regulation=='0']) # 0.77

shapiro.test(T$Fmax25) # 0.55
shapiro.test(T$Fmax25[T$regulation=='1']) # 0.46
shapiro.test(T$Fmax25[T$regulation=='0']) # 0.57

shapiro.test(T$Dmax25) # 0.15
shapiro.test(T$Dmax25[T$regulation=='1']) # 0.23
shapiro.test(T$Dmax25[T$regulation=='0']) # 0.62

shapiro.test(T$Rrng) # 0.12
shapiro.test(T$Rrng[T$regulation=='1']) # 0.29
shapiro.test(T$Rrng[T$regulation=='0']) # 0.19

# ---------------------------------------------------------------------------------------------------
# ----------------------------- test de Student (t.test)---------------------------------------------
# ---------------------------------------------------------------------------------------------------

t.test(T$Mmo[T$regulation=='1'],T$Mmo[T$regulation=='0'], alternative = "two.sided")       # 1 = 3.72 / 2 = 3.40
t.test(T$Fmax25[T$regulation=='1'],T$Fmax25[T$regulation=='0'], alternative = "two.sided") # 1 = 7.40 / 2 = 0.75
t.test(T$Dmax25[T$regulation=='1'],T$Dmax25[T$regulation=='0'], alternative = "two.sided") # 1 = 4.60 / 2 = -0.25
t.test(T$Rrng[T$regulation=='1'],T$Rrng[T$regulation=='0'], alternative = "two.sided")     # 1 = -1.17 / 2 = -2.15

# ---------------------------------------------------------------------------------------------------
# ----------------------------- ANOVA ---------------------------------------------------------------
# ---------------------------------------------------------------------------------------------------

modelanov1 <- aov(T$Mmo~ T$regulation)
Anova(modelanov1) # Fvalue = 0.04 / Pr(>F) = 0.85  == > no difference

modelanov2 <- aov(T$Fmax25~ T$regulation)
Anova(modelanov2) # Fvalue = 1.98 / Pr(>F) = 0.20  == > no difference

modelanov3 <- aov(T$Dmax25~ T$regulation)
Anova(modelanov3) # Fvalue = 2.8 / Pr(>F) = 0.14  == > no difference

modelanov4 <- aov(T$Rrng~ T$regulation)
Anova(modelanov4) # Fvalue = 4.75 / Pr(>F) = 0.07  == > no difference

# ---------------------------------------------------------------------------------------------------
# -----------------------------illustrate the results------------------------------------------------
# ---------------------------------------------------------------------------------------------------

p1 <- ggplot(T, aes(x = regulation, y = Mmo, fill = regulation))+
        stat_boxplot(geom ='errorbar', width = 0.4)+
        geom_boxplot(fatten = 2, outlier.size= 2)+
        theme_classic()+
        scale_fill_manual(values = c('0'='turquoise4', '1'='tomato3'))+
        scale_x_discrete(breaks=c("0","1"),
                   labels=c("lake", "reservoir"))+
        xlab("")+
        ylab("\n ?? mean temperature (°C) \n")+
        expand_limits(y = c(0, 8))+
        stat_summary(fun.y = mean, shape = 8, size = 0.6)+
        theme(axis.text = element_text(color = 'black', size = 11))+
        theme(axis.title = element_text(size = 12))+
        theme(legend.position = "none")+
        theme(axis.line = element_line(size = .8))+
        theme(axis.ticks = element_line(size = .8))+
        ggtitle("a) ?? Mmo")+
        theme(plot.title = element_text(color="black", size=14, hjust = 0.5))

p2 <- ggplot(T, aes(x = regulation, y = Fmax25, fill = regulation))+
        stat_boxplot(geom ='errorbar', width = 0.4)+
        geom_boxplot(fatten = 2, outlier.size= 2)+
        theme_classic()+
        scale_fill_manual(values = c('0'='turquoise4', '1'='tomato3'))+
        scale_x_discrete(breaks=c("0","1"),
                   labels=c("lake", "reservoir"))+
        xlab("")+
        ylab(" \n  ?? Number of days where daily maximum \n water temperature was above 25 °C")+
        expand_limits(y = c(-5, 20))+
        stat_summary(fun.y = mean, shape = 8, size = 0.6)+
        theme(axis.text = element_text(color = 'black', size = 11))+
        theme(axis.title = element_text(size = 12))+
        theme(legend.position = "none")+
        theme(axis.line = element_line(size = .8))+
        theme(axis.ticks = element_line(size = .8))+
        ggtitle("b) ?? Fmax25")+
        theme(plot.title = element_text(color="black", size=14, hjust = 0.5))

p3 <- ggplot(T, aes(x = regulation, y = Dmax25, fill = regulation))+
        stat_boxplot(geom ='errorbar', width = 0.4)+
        geom_boxplot(fatten = 2, outlier.size= 2)+
        theme_classic()+
        scale_fill_manual(values = c('0'='turquoise4', '1'='tomato3'))+
        scale_x_discrete(breaks=c("0","1"),
                   labels=c("lake", "reservoir"))+
        xlab("")+
        ylab(" ?? Maximum number of consecutive days \n where daily maximum water temperature  \n was above 25°C")+
        expand_limits(y = c(-5, 15))+
        stat_summary(fun.y = mean, shape = 8, size = 0.6)+
        theme(axis.text = element_text(color = 'black', size = 11))+
        theme(axis.title = element_text(size = 12))+
        theme(legend.position = "none")+
        theme(axis.line = element_line(size = .8))+
        theme(axis.ticks = element_line(size = .8))+
        ggtitle("c) ?? Dmax25")+
        theme(plot.title = element_text(color="black", size=14, hjust = 0.5))

p4 <- ggplot(T, aes(x = regulation, y = Rrng, fill = regulation))+
        stat_boxplot(geom ='errorbar', width = 0.4)+
        geom_boxplot(fatten = 2, outlier.size= 2)+
        theme_classic()+
        scale_fill_manual(values = c('0'='turquoise4', '1'='tomato3'))+
        scale_x_discrete(breaks=c("0","1"),
                   labels=c("lake", "reservoir"))+
        xlab("")+
        ylab("  \n ?? Mean daily range (°C)  \n")+
        expand_limits(y = c(-3, 3))+
        stat_summary(fun.y = mean, shape = 8, size = 0.6)+
        theme(axis.text = element_text(color = 'black', size = 11))+
        theme(axis.title = element_text(size = 12))+
        theme(legend.position = "none")+
        theme(axis.line = element_line(size = .8))+
        theme(axis.ticks = element_line(size = .8))+
        ggtitle("d) ?? Rrng")+
        theme(plot.title = element_text(color="black", size=14, hjust = 0.5))

p5 <- grid.arrange(p1,p2,p3,p4, ncol = 2)
ggsave(filename = "result_figure_2.png", plot = p5, width = 15, height = 16, dpi = 1200, units = "cm")