#####################################################
# Marco FW Gauger
# 2021/02/23
#
# This is the script to obtain table S 1 of the supplementary material 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
#####################################################

Sys.setlocale("LC_ALL","English")

ifelse("mgcv" %in% rownames(installed.packages()) == FALSE, 
       install.packages("mgcv", dependencies = T),
       ("mgcv is already installed"))
ifelse("bbmle" %in% rownames(installed.packages()) == FALSE, 
       install.packages("bbmle", dependencies = T),
       ("bbmle is already installed"))
ifelse("plyr" %in% rownames(installed.packages()) == FALSE, 
       install.packages("plyr", dependencies = T),
       ("plyr is already installed"))
ifelse("AICcmodavg" %in% rownames(installed.packages()) == FALSE, 
       install.packages("AICcmodavg", dependencies = T),
       ("AICcmodavg is already installed"))


analysed_data_hour <- read.csv("Supplemental Data S1.csv", sep=";")
analysed_data_hour$cluster_hcpc <- as.factor(analysed_data_hour$cluster_hcpc)



require(mgcv)
require(bbmle)
require(plyr)
require(AICcmodavg)

dp10m_AIC_00 <-
  gam(dp10m_ELAP1~ 1,
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_01 <-
  gam(dp10m_ELAP1~ s(SST),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_02 <-
  gam(dp10m_ELAP1~ s(SST_BLAP),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_03 <-
  gam(dp10m_ELAP1~ s(delta_SST),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_04 <-
  gam(dp10m_ELAP1~ s(CHL_BLAP),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_05 <-
  gam(dp10m_ELAP1~ s(lunar_phase, bs="cc", k=6), data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_06 <-
  gam(dp10m_ELAP1~ s(flow), data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_07 <-
  gam(dp10m_ELAP1~ s(derivate_tide),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_08 <-
  gam(dp10m_ELAP1~  s(tide),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_09 <-
  gam(dp10m_ELAP1~ effort ,
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_10 <-
  gam(dp10m_ELAP1~ s(h, bs="cc", k=10),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_11 <-
  gam(dp10m_ELAP1~ s(depth, k=4),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_12 <-
  gam(dp10m_ELAP1~ s(distance_BLAP, k=4),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_13 <-
  gam(dp10m_ELAP1~ s(distance_coast, k = 4),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_14 <-
  gam(dp10m_ELAP1~ s(distance_mangrove, k = 4),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_15 <-
  gam(dp10m_ELAP1~ s(SST) + s(lunar_phase, bs="cc", k=6) ,
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_16 <-
  gam(dp10m_ELAP1~ s(SST)+ s(tide) ,
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_17 <-
  gam(dp10m_ELAP1~ s(SST) + s(lunar_phase, bs="cc", k=6) + s(tide),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_18 <-
  gam(dp10m_ELAP1~   s(SST) + s(lunar_phase, bs="cc", k=6) + s(flow),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_19 <-
  gam(dp10m_ELAP1~ s(SST) + s(lunar_phase, bs="cc", k=6) + s(derivate_tide),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_20 <-
  gam(dp10m_ELAP1~ s(SST) + s(lunar_phase, bs="cc", k=6) + te(flow),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_21 <-
  gam(dp10m_ELAP1~ s(SST) + s(lunar_phase, bs="cc", k=6) + te(flow, tide, bs=c("tp","tp"), k=c(10,10)),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_22 <-
  gam(dp10m_ELAP1~ s(SST) + s(lunar_phase, bs="cc", k=6) + te(derivate_tide, tide, bs=c("tp","tp"), k=c(10,10)),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_23 <-
  gam(dp10m_ELAP1~ s(lunar_phase, bs="cc", k=6) + te(flow, tide, bs=c("tp","tp"), k=c(10,10)) + s(SST, k=10, bs="tp", m=2),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_24 <-
  gam(dp10m_ELAP1~ s(SST) + s(lunar_phase, bs="cc", k=6) + s(derivate_tide) + s(h, bs="cc",k=6),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_25 <-
  gam(dp10m_ELAP1~ s(SST) + s(lunar_phase, bs="cc", k=6) + s(derivate_tide) + s(tide) + s(h, bs="cc",k=6),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_26 <-
  gam(dp10m_ELAP1~ s(SST) + s(lunar_phase, bs="cc", k=6) + te(derivate_tide, tide, bs=c("tp","tp"), k=c(10,10)) + effort,
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))


dp10m_AIC_27 <-
  gam(dp10m_ELAP1~ s(SST, k=10, bs="tp", m=2) + s(lunar_phase, bs="cc", k=6) + te(derivate_tide, tide, bs=c("tp","tp"), k=c(10,10)) + s(h, bs="cc") + effort,
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_28 <-
  gam(dp10m_ELAP1~ s(lunar_phase, bs="cc", k=6) + te(derivate_tide, tide, bs=c("tp","tp"), k=c(10,10)) + s(SST, k=10, bs="tp", m=2) + s(h, bs="cc") + effort + s(depth, k=4),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))


dp10m_AIC_29 <-
  gam(dp10m_ELAP1~ s(lunar_phase, bs="cc", k=6) + te(derivate_tide, tide, bs=c("tp","tp"), k=c(10,10)) + s(SST, k=10, bs="tp", m=2) + s(h, bs="cc") + effort + s(depth, k=4) + s(delta_SST),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_30 <-
  gam(dp10m_ELAP1~ s(lunar_phase, bs="cc", k=6) + te(derivate_tide, tide, bs=c("tp","tp"), k=c(10,10)) + s(SST, k=10, bs="tp", m=2) + s(h, bs="cc") + effort + s(depth, k=4) + s(distance_BLAP, k=4),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))


dp10m_AIC_31 <-
  gam(dp10m_ELAP1~ s(lunar_phase, bs="cc", k=6) + te(derivate_tide, tide, bs=c("tp","tp"), k=c(10,10)) + s(SST, k=10, bs="tp", m=2) + s(h, bs="cc") + effort + s(depth, k=4) + s(distance_BLAP, k=4) + s(distance_coast, k=5),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_32 <-
  gam(dp10m_ELAP1~ s(lunar_phase, bs="cc", k=6) + te(derivate_tide, tide, bs=c("tp","tp"), k=c(10,10)) + s(SST, k=10, bs="tp", m=2) + s(h, bs="cc") + effort + s(depth, k=4) + s(distance_BLAP, k=4) + s(distance_coast, k=5) + s(distance_mangrove, k=5),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_AIC_33 <-
  gam(dp10m_ELAP1~ s(lunar_phase, bs="cc", k=6) + te(derivate_tide, tide, bs=c("tp","tp"), k=c(10,10)) + s(SST, k=10, bs="tp", m=2) + s(h, bs="cc") + effort + s(depth, k=4) + s(distance_BLAP, k=4) + s(distance_coast, k=5) + s(delta_SST),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

################

jj <- bbmle::AIC(dp10m_AIC_00,
                 dp10m_AIC_01, 
                 dp10m_AIC_02, 
                 dp10m_AIC_03,
                 dp10m_AIC_04,
                 dp10m_AIC_05,
                 dp10m_AIC_06,
                 dp10m_AIC_07,
                 dp10m_AIC_08,
                 dp10m_AIC_09,
                 dp10m_AIC_10,
                 dp10m_AIC_11,
                 dp10m_AIC_12,
                 dp10m_AIC_13,
                 dp10m_AIC_14,
                 dp10m_AIC_15,
                 dp10m_AIC_16,
                 dp10m_AIC_17,
                 dp10m_AIC_18,
                 dp10m_AIC_19,
                 dp10m_AIC_20,
                 dp10m_AIC_21,
                 dp10m_AIC_22,
                 dp10m_AIC_23,
                 dp10m_AIC_24,
                 dp10m_AIC_25,
                 dp10m_AIC_26,
                 dp10m_AIC_27,
                 dp10m_AIC_28,
                 dp10m_AIC_29,
                 dp10m_AIC_30,
                 dp10m_AIC_31,
                 dp10m_AIC_32,
                 dp10m_AIC_33)
View(jj)

models <- list(dp10m_AIC_00,dp10m_AIC_01,dp10m_AIC_02,dp10m_AIC_03,dp10m_AIC_04,dp10m_AIC_05,dp10m_AIC_06,dp10m_AIC_07,dp10m_AIC_08,dp10m_AIC_09,dp10m_AIC_10,dp10m_AIC_11,dp10m_AIC_12,dp10m_AIC_13,dp10m_AIC_14,dp10m_AIC_15,dp10m_AIC_16,dp10m_AIC_17,dp10m_AIC_18,dp10m_AIC_19,dp10m_AIC_20,dp10m_AIC_21,dp10m_AIC_22,dp10m_AIC_23,dp10m_AIC_24,dp10m_AIC_25,dp10m_AIC_26,dp10m_AIC_27,dp10m_AIC_28,dp10m_AIC_29,dp10m_AIC_30,dp10m_AIC_31,dp10m_AIC_32)

mof <- 
  data.frame(cbind(model = as.character(sapply(1:length(models), function(x) deparse(formula(models[[x]])))),ldply(models, function(x) cbind(AIC = AIC(x)))
))
mof <- cbind(mof, (bbmle::AICtab(models, weights=T, sort=F)))

mof$comment <- c(rep("",nrow(mof)-3),
                 "most parsimonious model",
                 "depth and distance_mangrove not significant",
                 "delta_SST not significant")

#View(mof)

write.csv(mof, "clipboard")