#*---------------------------------------------------------------------------------------------------
# | SCRIPT NAME: glmm_TMB_dams_impact
# | DATE CREATED: 2020-11-26
# | CREATED BY: Mathieu Auffray
# *--------------------------------------------------------------------------------------------------
# | OBJECTIVE: this script evaluates the fit of the linear mixed model. It evaluates the effect of 
#              the position (upstream-downstream) in relation to the reservoir (i.e. waterbodie regulated 
#              by a dam) on the water temperature in Gatineau Park in August 2019 and 2020 by means of four
#              thermal regime measurements:
#                   -  monthly average temperature in August in ° C (Mmo)
#                   -  frequency of hot events (> 25 ° C) in number of days (Fmax25)
#                   -  Duration of hot events (>> 25 ° C) in number of days (Dmax25)
#                   -  Daily rate of change in ° C (Rrng)
# *--------------------------------------------------------------------------------------------------
# | DATA: analysis\hypothesis_1_dams_impact\thermal_regime_dams.csv
# *--------------------------------------------------------------------------------------------------

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

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

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

# data transformation

T$year <- as.factor(T$year)
T$Fmax25 <- as.numeric((T$Fmax25))
T$Dmax25 <- as.numeric((T$Dmax25))

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

## check pairing downstream/upstream

xtabs(~reservoir_name + Pos, data=T) # Each reservoir has as many sites upstream and downstream

## shapiro test

shapiro.test(T$Mmo) # 0.11 
shapiro.test(T$Mmo[T$Pos =='upstream']) # 0.05
shapiro.test(T$Mmo[T$Pos =='downstream']) # 0.63

shapiro.test(T$Fmax25) # no and the transformation does not allow to obtain normality

shapiro.test(T$Dmax25) # no and the transformation does not allow to obtain normality

shapiro.test(T$Rrng) # 0.14

# ---------------------------------------------------------------------------------------------------
# -----------------------------building model--------------------------------------------------------
# ---------------------------------------------------------------------------------------------------

# ----------- building model monthly average temperature in August in ° C (Mmo) ---------------------

mod1 <- glmmTMB( Mmo ~ Pos  + (1|year/reservoir_name),family = "gaussian", data=T)
summary(mod1)

# --------------------------- model verification ----------------------------------------------------

shapiro.test(resid(mod1))
leveneTest(resid(mod1),T$Pos) # Fvalue = 0.08 / Pr(>F) = 0.78

# ----------- building model frequency of hot events (> 25 ° C) in number of days (Fmax25) ---------

mod2 <- glmmTMB( Fmax25~ Pos  + (1|year/reservoir_name),family = "poisson", data=T)
summary(mod2)

# --------------------------- model verification ----------------------------------------------------

shapiro.test(resid(mod2)) # 0.37
leveneTest(resid(mod2),T$Pos) # Fvalue = 3.9 / Pr(>F) = 0.07

# -------------building model Duration of hot events (>> 25 ° C) in number of days (Dmax25) ---------

mod3 <- glmmTMB( Dmax25~ Pos  + (1|year/reservoir_name),family = "poisson", data=T)
summary(mod3)

# --------------------------- model verification ----------------------------------------------------

shapiro.test(resid(mod3)) # 0.02
leveneTest(resid(mod3),T$Pos) # Fvalue = 2.3 / Pr(>F) = 0.16

# ----------- building model Daily rate of change in ° C (Rrng) -------------------------------------

mod4 <- glmmTMB( Rrng ~ Pos  + (1|year/reservoir_name),family = "gaussian", data=T)
summary(mod4)

# --------------------------- model verification ----------------------------------------------------

shapiro.test(resid(mod4)) # 0.32
leveneTest(resid(mod4),T$Pos) # Fvalue = 0.29 / Pr(>F) = 0.59

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

p1 <- ggplot(cbind(T, pred = predict(mod1)), aes(x = Pos, y = Mmo, fill = Pos))+
          stat_boxplot(geom ='errorbar', width = 0.4)+
          geom_boxplot(fatten = 2, outlier.size= 2)+
          theme_classic()+
          scale_fill_manual(values = c('upstream'='steelblue3', 'downstream'='tomato3'))+
          xlab("")+
          ylab("\n mean temperature (°C) \n")+
          expand_limits(y = c(16, 24))+
          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(cbind(T, pred = predict(mod2)), aes(x = Pos, y = Fmax25, fill = Pos))+
          stat_boxplot(geom ='errorbar', width = 0.4)+
          geom_boxplot(fatten = 2, outlier.size= 2)+
          theme_classic()+
          scale_fill_manual(values = c('upstream'='steelblue3', 'downstream'='tomato3'))+
          xlab("")+
          ylab(" \n Number of days where daily maximum water temperature was above 25 °C")+
          expand_limits(y = c(0, 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(cbind(T, pred = predict(mod3)), aes(x = Pos, y = Dmax25, fill = Pos))+
          stat_boxplot(geom ='errorbar', width = 0.4)+
          geom_boxplot(fatten = 2, outlier.size= 2)+
          theme_classic()+
          scale_fill_manual(values = c('upstream'='steelblue3', 'downstream'='tomato3'))+
          xlab("")+
          ylab("Maximum number of consecutive days where daily maximum water temperature  \n was above 25°C")+
          expand_limits(y = c(0, 15))+
          stat_summary(fun.y = mean, shape = 8, size = .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(cbind(T, pred = predict(mod4)), aes(x = Pos, y = Rrng, fill = Pos))+
          stat_boxplot(geom ='errorbar', width = 0.4)+
          geom_boxplot(fatten = 2, outlier.size= 2)+
          theme_classic()+
          scale_fill_manual(values = c('upstream'='steelblue3', 'downstream'='tomato3'))+
          xlab("")+
          ylab(" \n Mean daily range (°C) \n")+
          expand_limits(y = c(1, 6))+
          stat_summary(fun.y = mean, shape = 8, size = .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 = "GLMTMB_result.png", plot = p5, width = 15, height = 16, dpi = 1200, units = "cm")