#####################################################
# Marco FW Gauger
# 2021/02/23
#
# This is the script to obtain tab 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("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"))

require(mgcv)
require(bbmle)
########
analysed_data_hour <- read.csv("Supplemental Data S1.csv", sep=";")
analysed_data_hour$cluster_hcpc <- as.factor(analysed_data_hour$cluster_hcpc)
sum(analysed_data_hour$dp10m_ELAP1)
length(analysed_data_hour$dp10m_ELAP1[analysed_data_hour$dp10m_ELAP1>0])

final <-
  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") +
        s(h, k=10, 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)))
AIC(final)

Sys.setlocale("LC_ALL","English")
Sys.setenv(TZ='UTC')
dp10m_mod_cluster_S_h <-
  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") +
        s(h, k=10, cluster_hcpc, bs="fs",xt=list(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_mod_cluster_I_h <-
  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") +
        s(h, k=10, by = cluster_hcpc, bs="cc", m=1) + 
        effort +
        s(depth, k=4) +
        s(distance_BLAP, k=4) + 
        s(cluster_hcpc, bs="re"),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_mod_cluster_GS_h <-
  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") +
        s(h, k=6, bs="cc") + 
        s(h, k=10, cluster_hcpc, bs="fs",xt=list(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)))

# Note that s(cluster_hcpc, bs="re") has to be explicitly included here, as the h by cluster smoother does not include an intercept
dp10m_mod_cluster_GI_h <-
  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") +
        s(h, k=6, bs="cc", m=2) + 
        s(h, k=10, by = cluster_hcpc, bs="cc", m=1) + 
        effort +
        s(depth, k=4) +
        s(distance_BLAP, k=4) + 
        s(cluster_hcpc, bs="re"),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))


dp10m_mod_cluster_GS_d2 <-
  gam(dp10m_ELAP1~ s(lunar_phase, bs="cc", k=6) +
        te(derivate_tide, tide, bs=c("tp","tp"),
           k=c(10,10)) + 
        t2(derivate_tide, tide, cluster_hcpc,
           bs=c("tp","tp","re"), k=c(10,10)) + 
        s(SST, k=10, bs="tp") +
        s(h, k=10, 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_mod_cluster_S_d2 <-
  gam(dp10m_ELAP1~ s(lunar_phase, bs="cc", k=6) +
        #te(derivate_tide, tide, bs=c("tp","tp"),
        #   k=c(10,10)) + 
        t2(derivate_tide, tide, cluster_hcpc,
           bs=c("tp","tp","re"), k=c(10,10)) + 
        s(SST, k=10, bs="tp") +
        s(h, k=10, 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_mod_cluster_GI_d2 <-
  gam(dp10m_ELAP1~ s(lunar_phase, bs="cc", k=6) +
        te(derivate_tide, tide, bs=c("tp","tp"),
           k=c(10,10)) + 
        te(derivate_tide, tide, by=cluster_hcpc,
           bs=c("tp","tp"), k=c(10,10)) + 
        s(SST, k=10, bs="tp") +
        s(h, k=10, bs="cc") +
        effort +
        s(depth, k=4) +
        s(distance_BLAP, k=4) + 
        s(cluster_hcpc, bs="re"),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

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

dp10m_mod_cluster_GS_s2 <-
  gam(dp10m_ELAP1~ 
        s(lunar_phase, k=6, bs="cc") + 
        te(derivate_tide, tide, bs=c("tp","tp"),
           k=c(10,10)) + 
        s(SST, k=6, bs="tp") + 
        s(SST, cluster_hcpc, k=10, bs="fs",xt=list(bs="tp")) + 
        s(h, k=10, 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_mod_cluster_S_s2 <-
  gam(dp10m_ELAP1~ 
        s(lunar_phase, k=6, bs="cc") + 
        te(derivate_tide, tide, bs=c("tp","tp"),
           k=c(10,10)) + 
        s(SST, cluster_hcpc, k=10, bs="fs",xt=list(bs="tp")) + 
        s(h, k=10, 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_mod_cluster_I_s2 <-
  gam(dp10m_ELAP1~ 
        s(lunar_phase, k=6, bs="cc") + 
        te(derivate_tide, tide, bs=c("tp","tp"),
           k=c(10,10)) + 
        s(SST, by=cluster_hcpc, k=10, bs="tp", m=1) + 
        s(h, k=10, bs="cc") + 
        effort +
        s(depth, k=4) +
        s(distance_BLAP, k=4) + 
        s(cluster_hcpc, bs="re"),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

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





####
dp10m_mod_cluster_S_l2 <-
  gam(dp10m_ELAP1~ s(lunar_phase, cluster_hcpc, k=6, 
                     bs="fs",xt=list(bs="cc")) +
        te(derivate_tide, tide, bs=c("tp","tp"),
           k=c(10,10)) + 
        s(SST, k=10, bs="tp") +
        s(h, k=10, 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_mod_cluster_I_l2 <-
  gam(dp10m_ELAP1~ s(lunar_phase, by = cluster_hcpc, bs="cc",  
                     m=1, k=6) +
        te(derivate_tide, tide, bs=c("tp","tp"),
           k=c(10,10)) + 
        s(SST, k=10, bs="tp") +
        s(h, k=10, bs="cc") + 
        effort +
        s(depth, k=4) +
        s(distance_BLAP, k=4) + 
        s(cluster_hcpc, bs="re"),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

dp10m_mod_cluster_GS_l2 <-
  gam(dp10m_ELAP1~ s(lunar_phase, bs="cc", k=6) +
        s(lunar_phase, cluster_hcpc, k=6, bs="fs",xt=list(bs="cc")) + 
        te(derivate_tide, tide, bs=c("tp","tp"),
           k=c(10,10)) + 
        s(SST, k=10, bs="tp") +
        s(h, k=10, 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)))

# Note that s(cluster_hcpc, bs="re") has to be explicitly included here, as the h by cluster smoother does not include an intercept
dp10m_mod_cluster_GI_l2 <-
  gam(dp10m_ELAP1~ s(lunar_phase, bs="cc", k=6) +
        s(lunar_phase, by = cluster_hcpc, bs="cc", k=6, m=1) + 
        te(derivate_tide, tide, bs=c("tp","tp"),
           k=c(10,10)) + 
        s(SST, k=10, bs="tp") +
        s(h, k=10, bs="cc") + 
        effort +
        s(depth, k=4) +
        s(distance_BLAP, k=4) + 
        s(cluster_hcpc, bs="re"),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))


summary(final)
summary(dp10m_mod_cluster_GS_h)
summary(dp10m_mod_cluster_S_h)
summary(dp10m_mod_cluster_GI_h)
summary(dp10m_mod_cluster_I_h)
summary(dp10m_mod_cluster_GS_s2)
summary(dp10m_mod_cluster_S_s2)
summary(dp10m_mod_cluster_GI_s2)
summary(dp10m_mod_cluster_I_s2)
summary(dp10m_mod_cluster_GS_d2)
summary(dp10m_mod_cluster_S_d2)
summary(dp10m_mod_cluster_GI_d2)
summary(dp10m_mod_cluster_I_d2)
summary(dp10m_mod_cluster_GS_l2)
summary(dp10m_mod_cluster_S_l2)
summary(dp10m_mod_cluster_GI_l2)
summary(dp10m_mod_cluster_I_l2)

MODELS <- list(dp10m_mod_cluster_GS_h, dp10m_mod_cluster_S_h, dp10m_mod_cluster_GI_h, dp10m_mod_cluster_I_h, dp10m_mod_cluster_GS_s2, dp10m_mod_cluster_S_s2, dp10m_mod_cluster_GI_s2, dp10m_mod_cluster_I_s2, dp10m_mod_cluster_GS_d2, dp10m_mod_cluster_S_d2, dp10m_mod_cluster_GI_d2, dp10m_mod_cluster_I_d2, dp10m_mod_cluster_GS_l2, dp10m_mod_cluster_S_l2, dp10m_mod_cluster_GI_l2, dp10m_mod_cluster_I_l2)
MODELS_names <- c("dp10m_mod_cluster_GS_h", "dp10m_mod_cluster_S_h", "dp10m_mod_cluster_GI_h", "dp10m_mod_cluster_I_h", "dp10m_mod_cluster_GS_s2", "dp10m_mod_cluster_S_s2", "dp10m_mod_cluster_GI_s2", "dp10m_mod_cluster_I_s2", "dp10m_mod_cluster_GS_d2", "dp10m_mod_cluster_S_d2", "dp10m_mod_cluster_GI_d2", "dp10m_mod_cluster_I_d2", "dp10m_mod_cluster_GS_l2", "dp10m_mod_cluster_S_l2", "dp10m_mod_cluster_GI_l2", "dp10m_mod_cluster_I_l2")
AIC_models <- AIC(dp10m_mod_cluster_GS_h, dp10m_mod_cluster_S_h, dp10m_mod_cluster_GI_h, dp10m_mod_cluster_I_h, dp10m_mod_cluster_GS_s2, dp10m_mod_cluster_S_s2, dp10m_mod_cluster_GI_s2, dp10m_mod_cluster_I_s2, dp10m_mod_cluster_GS_d2, dp10m_mod_cluster_S_d2, dp10m_mod_cluster_GI_d2, dp10m_mod_cluster_I_d2, dp10m_mod_cluster_GS_l2, dp10m_mod_cluster_S_l2, dp10m_mod_cluster_GI_l2, dp10m_mod_cluster_I_l2)
bbmle::AICtab(dp10m_mod_cluster_GS_h, dp10m_mod_cluster_S_h, dp10m_mod_cluster_GI_h, dp10m_mod_cluster_I_h, dp10m_mod_cluster_GS_s2, dp10m_mod_cluster_S_s2, dp10m_mod_cluster_GI_s2, dp10m_mod_cluster_I_s2, dp10m_mod_cluster_GS_d2, dp10m_mod_cluster_S_d2, dp10m_mod_cluster_GI_d2, dp10m_mod_cluster_I_d2, dp10m_mod_cluster_GS_l2, dp10m_mod_cluster_S_l2, dp10m_mod_cluster_GI_l2, dp10m_mod_cluster_I_l2,weights=T,  sort = F)

number_invalid_parameters <- invalid_parameters <-  number_parameters <- rep(NA, length(MODELS))
for(k in 1:length(MODELS)){
  number_invalid_parameters[k] <- (length(summary(MODELS[[k]])$s.pv)-table(summary(MODELS[[k]])$s.pv >.05)[1])
  
  number_parameters[k] <- length(summary(MODELS[[k]])$s.pv)
  invalid_parameters[k] <- 
    (paste(rownames((summary(MODELS[[k]])[24])$s.table)[c(1:length(summary(MODELS[[k]])[24]$s.table[,4])) [summary(MODELS[[k]])[24]$s.table[,4] >.05]], collapse=", "))
}
model_inf <- 
  data.frame( model_name = MODELS_names,
              AIC = AIC_models$AIC,
              dAIC = bbmle::AICtab(MODELS, sort=F)$dAIC,
              AIC_weight=round(bbmle::AICtab(MODELS, sort=F, 
                                               weights=T)$weight, 
                                 digits=3),
              number_invalid_parameters =  number_invalid_parameters,
              invalid_parameters = invalid_parameters,
              number_parameters = number_parameters)  
View(model_inf)

write.table(model_inf, "clipboard", sep="\t")

# summarys used to obtain supplemtal tables 2 and 3
summary(dp10m_mod_cluster_I_h)
summary(dp10m_mod_cluster_S_h)