#*------------------------------------------------------------------
# | SCRIPT NAME: GAM_analysis
# | DATE CREATED: 2020-11-30
# | CREATED BY: Mathieu Auffray
# *----------------------------------------------------------------
# | OBJECTIVE: this script buid a model GAM of water temperature for
#              the throughout the Gatineau Park
# *------------------------------------------------------------------
# | DATA: measured_site.csv
#         predicted_site.csv
# *------------------------------------------------------------------

library(mgcv)
library(gratia)
library(DHARMa)
library(Metrics)


library(units)



#library(ggeffects)
# library(openxlsx)


library(udunits2)
library(gridExtra)
#library(visreg)
# library(ggplot2)
# library(sjPlot)

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

measure <- read.table("measured_site.csv", dec = ",", sep = ";", header = TRUE)
predict <- read.table("predicted_site.csv", dec = ",", sep = ";", header = TRUE)


# data transformation

str(measure)
str(predict)

measure$id          <- as.factor(measure$id)
measure$WS          <- as.factor(measure$WS)
measure$site_name   <- as.factor(measure$site_name)

predict$id          <- as.factor(predict$id)
predict$WS          <- as.factor(predict$WS)
predict$site_name   <- as.factor(predict$site_name)


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

##  data "measure"

shapiro.test(measure$E) # 0.34

shapiro.test(measure$Slope) # 0.0002  - We try the LOG transformation
shapiro.test(log(measure$Slope)) # NA - slope values are equal to 0. the LOG transformation therefore does not work. We therefore try the SQRT transformation
shapiro.test(sqrt(measure$Slope)) # 0.13

shapiro.test(measure$Dist) # 0.00046 - We try the LOG transformation
shapiro.test(log(measure$Dist)) # 0.002 - We try the SQRT transformation
shapiro.test(sqrt(measure$Dist)) # 0.06

shapiro.test(measure$DA) # 4e-12 - We try the LOG transformation
shapiro.test(log(measure$DA)) # 0.02 - We try the SQRT transformation
shapiro.test(sqrt(measure$DA))# 2.497e-09 - doesn't improve normality. We choose the LOG transformation

shapiro.test(measure$Wb_surf) # 0.0006703 - We try the LOG transformation
shapiro.test(log(measure$Wb_surf)) # 0.09261

shapiro.test(measure$TA) # 0.001288 - We try the LOG transformation
shapiro.test(log(measure$TA)) # 0.001267  - We try the SQRT transformation
shapiro.test(sqrt(measure$TA)) # 0.001279 - not transformation

shapiro.test(measure$Mmo) # 0.7941

# data transformation to obtain normality

measure$E_sqrt <- sqrt(measure$E)
measure$Slope_sqrt <- sqrt(measure$Slope)
measure$Dist_sqrt <- sqrt(measure$Dist)
measure$DA_log <- log10(measure$DA)
measure$Wb_surf_log <- log(measure$Wb_surf)

predict$E_sqrt <- sqrt(predict$E)
predict$Slope_sqrt <- sqrt(predict$Slope)
predict$Dist_sqrt <- sqrt(predict$Dist)
predict$DA_log <- sqrt(predict$DA)
predict$Wb_surf_log <- log(predict$Wb_surf)
## assign the value "54.005" to the variable "dist_lakeSQRT" so that the rivers section
## without an upstream lake has no effect on the model. The value corresponds to the distance where the
## lake will no longer have an effect on the water temperature

predict_temp1 <- predict[-which(is.na(predict$Dist_sqrt)),]
predict_temp2 <- predict[which(is.na(predict$Dist_sqrt)),]
predict_temp2$Dist_sqrt <- 54.005
predict <- rbind(predict_temp1,predict_temp2)

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

# model building

library(MASS)

set.seed(1010) # Set Seed so that same sample can be reproduced in future also ## 101-->good, 70% of dataframe
measure_b <- measure[which(is.na(measure$Dist_sqrt)),]
measure_c <- measure[-which(is.na(measure$Dist_sqrt)),]
measure_b$Dist_sqrt <- 54.005 ## assigning the value "54.005" to variable "dist_sqrt" so that the section of rivers with no lake upstream have no effect on the model

measure_d <- rbind(measure_b, measure_c)
measure_d <- measure_d[-which(is.na(measure_d$Wb_surf_log)),]

measure_d$WS_chiff <- as.numeric(measure_d$WS)
train <- measure_d[which(measure_d$WS_chiff<4), ]# training dataset with 70% included
test  <- measure_d[-which(measure_d$WS_chiff<4), ]# training dataset with 70% excluded

## modele with 70% data included

modelGAM <- gam(Mmo ~
                 s(E, k = 3, bs = "ts")+
                 s(Slope_sqrt, k = 3, bs = "ts")+
                 s(Dist_sqrt, pc = 54.005,k = 3, bs = "ts")+
                 s(DA_log, k =3, bs = "ts")+
                 s(Wb_surf, k = 3, bs = "ts")+
                 s(TA, k = 3, bs = "ts")+
                 s(WS, bs = "re")
               ,select = T,data = train, family = gaussian(link= "identity"), method = 'REML')


summary(modelGAM)
draw(modelGAM,residuals = T)
qq.gam(modelGAM,rep=100,level=1,pch=19,cex=.8)

## Model with all the data

modelGAM_2 <- gam(Mmo ~
                  s(E, k = 3, bs = "ts")+
                  s(Slope_sqrt, k = 3, bs = "ts")+
                  s(Dist_sqrt, pc = 54.005,k = 3, bs = "ts")+
                  s(DA_log, k =3, bs = "ts")+
                  s(Wb_surf, k = 3, bs = "ts")+
                  s(TA, k = 3, bs = "ts")+
                  s(WS, bs = "re")
                ,select = T,data = measure_d, family = gaussian(link= "identity"), method = 'REML')


summary(modelGAM_2)
draw(modelGAM_2,residuals = T)
qq.gam(modelGAM_2,rep=100,level=1,pch=19,cex=.8)

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

testDispersion(modelGAM)## 0.65
concurvity(modelGAM,full = F)##  concurvity of nonlinear parameters
appraise(modelGAM)
simulationOutput <- simulateResiduals(fittedModel = modelGAM, plot = T)

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

library(mgcViz)
b <- mgcViz::getViz(modelGAM)
o <- plot( mgcViz::sm(b, 3) )
plot1 <- o + mgcViz::l_fitLine(colour = "tomato3", size = 1)+
  mgcViz::l_ciLine(level = 0.95, colour = "grey40", linetype = 2, size = 1) + 
  mgcViz::l_points(shape = 19, size = 2.5, alpha = 1) + 
  ggplot2::theme_classic(base_size=18)+
  ggplot2::labs(x="Distance to an upstream lake \n (m expressed in square root)", y="Partial effect")+
  ggplot2::scale_x_continuous(breaks=seq(0,70,by=10))+
  ggplot2::scale_y_continuous(breaks=seq(-4,6,by=1))+
  ggplot2::ggtitle("a)")+
  ggplot2::theme(plot.title = element_text(color="black", hjust = 0.5))+
  ggplot2::theme(axis.line = element_line(size = 1))+
  ggplot2::theme(axis.ticks = element_line(size = 1))

p1<-plot1$ggObj

b <- mgcViz::getViz(modelG1)
o <- plot( mgcViz::sm(b, 4) )
plot2 <- o + mgcViz::l_fitLine(colour = "tomato3", size = 1)+
  mgcViz::l_ciLine(level = 0.95, colour = "grey40", linetype = 2, size = 1) + 
  mgcViz::l_points(shape = 19, size = 2.5, alpha = 1) + 
  ggplot2::theme_classic(base_size=18)+
  ggplot2::labs(x="Drainage area \n (ha expressed in logarithm)", y="")+
  ggplot2::scale_x_continuous(breaks=seq(0,5,by=1))+
  ggplot2::scale_y_continuous(breaks=seq(-4,4,by=1))+
  ggplot2::ggtitle("b)")+
  ggplot2::theme(plot.title = element_text(color="black", hjust = 0.5))+
  ggplot2::theme(axis.line = element_line(size = 1))+
  ggplot2::theme(axis.ticks = element_line(size = 1))

p2<-plot2$ggObj

b <- mgcViz::getViz(modelG1)
o <- plot( mgcViz::sm(b, 5) )
plot3 <- o + mgcViz::l_fitLine(colour = "tomato3", size = 1)+
  mgcViz::l_ciLine(level = 0.95, colour = "grey40", linetype = 2, size = 1) + 
  mgcViz::l_points(shape = 19, size = 2.5, alpha = 1) + 
  ggplot2::theme_classic(base_size=18)+
  ggplot2::labs(x="Watershed area occupied \n by waterbodies (ha)", y="")+
  ggplot2::scale_x_continuous(breaks=seq(0,1,by=0.05))+
  ggplot2::scale_y_continuous(breaks=seq(-4,4,by=1))+
  ggplot2::ggtitle("c)")+
  ggplot2::theme(plot.title = element_text(color="black", hjust = 0.5))+
  ggplot2::theme(axis.line = element_line(size = 1))+
  ggplot2::theme(axis.ticks = element_line(size = 1))

p3<-plot3$ggObj

p4 <-grid.arrange(p1,p2,p3, ncol = 3)
ggsave(filename = "Gam_result.png", plot = p4, width = 30, height = 20, dpi = 300, units = "cm")

# ---------------------------------------------------------------------------------------------------
# -----------------------------predict model --------------------------------------------------------
# ---------------------------------------------------------------------------------------------------

reponse <- predict(modelGAM,test,type="response",exclude="s(WS)",newdata.guaranteed=TRUE,na.action = na.omit)
bb<-cbind(test,reponse)
fitcorr <- format(cor(bb$Mmo, bb$reponse)^2, digits=4)
op <- par(c(lty="solid", col="blue"))
plot(bb[['Mmo']], bb[['reponse']], asp=1, xlab="Mmo", ylab="Predicted",xlim=c(13,25),ylim=c(13,25))
prline <- lm(bb$reponse ~ bb$Mmo)
abline(prline)
par(c(lty="dashed", col="black"))
abline(0, 1)
legend("bottomright",  sprintf(" Pseudo R-square=%s ", fitcorr),  bty="n")
title(main=paste( round(prline$coefficients[1],2),"+", round(prline$coefficients[2],2),"*","Mmo", sep =""))
grid()

bb$ecart <- bb$Mmo - bb$reponse
hist(bb$ecart)

# compare measured and predicted datasets

summary(measure_d$DA_log)
summary(predict$DA_log)
summary(train$DA_log)

summary(measure_d$Slope_sqrt)
summary(predict$Slope_sqrt)
summary(train$Slope_sqrt)

summary(measure_d$Wb_surf)
summary(predict$Wb_surf)
summary(train$Wb_surf)

hist(predict$Wb_surf)
hist(measure_d$Wb_surf)


# determinate predictions using the modelGAM


bb<-cbind(predict,reponse)
bb<-cbind(bb,reponse)

summary(bb$reponse)
summary(measure_d$Mmo)


# ---------------------------------------------------------------------------------------------------
# -----------------------------illustrate the results------------------------------------------------
# ---------------------------------------------------------------------------------------------------
install.packages(Rcpp)
library(Rcpp)

install.packages(sf)
library(sf)

projcrs <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
df <- st_as_sf(x = bb, coords = c("lat", "long"),crs = projcrs)
st_write(df, "prediction.gpkg")