# packages
library(lme4)
library(effects)
library(DHARMa)
library(MASS)
library(car)
library(ggplot2)
library(GGally)
library(glmmTMB)
library(pscl)
library(lmtest)
library(splines)
library(visreg)
library(knitr)
library(splines)
library(sjPlot)
library(performance)
library(rstatix)

# dataset: names Stats
Stats$prop <- Stats$No_colonies/3026
Stats$prop <- as.numeric(Stats$prop)
Stats$Year <- as.numeric(Stats$Year)
Stats$Month <- as.factor(Stats$Month)
Stats$PFM <- as.numeric(Stats$PFM)
Stats$Pos_PFM <- as.numeric(Stats$Pos_PFM) 
Stats$Atoll <- as.factor(Stats$Atoll)
Stats$No_colonies <- as.numeric(Stats$No_colonies)
Stats$DailyWindSpeed <- as.numeric(Stats$DailyWindSpeed)
Stats$DailyPpt <- as.numeric(Stats$DailyPpt) 
Stats$TideDepth <- as.numeric(Stats$TideDepth)
Stats$SST30 <- as.numeric(Stats$SST30)
Stats$SST90 <- as.numeric(Stats$SST90)
Stats$Date <- as.Date(as.character(Stats$Date),format = "%Y%m%d")
vif(lm(Pos_PFM ~ DailyWindSpeed + DailyPpt + SST30 + TideDepth, data = Stats))
cor(Stats[, c("DailyWindSpeed", "DailyPpt", "SST30", "SST90", "TideDepth")], use = "complete.obs")

#### pfm / no cols ~ year * atoll ####
Stats$Year2 <- Stats$Year
Stats$Year2 <- as.factor(Stats$Year2)
j <- glm.nb(Pos_PFM ~ Year2*Atoll, data = Stats)
plot(simulateResiduals(j))
summary(j)
plot(allEffects(j))

k <- glm.nb(No_colonies ~ Year2*Atoll, data = Stats)
plot(simulateResiduals(k))
summary(k)
plot(allEffects(k))


#### Histograms of No Colonies ####
hist(statsNM$No_colonies)
hist(statsB$No_colonies)
hist(Stats$No_colonies)

figurehNM <- ggplot(data = statsNM, aes(x = No_colonies)) +
  geom_histogram(binwidth = 50, fill = "gray51",color = "black", boundary = 0) +  
  scale_x_continuous(limits = c(0, 250), breaks = seq(0, 250, by = 50)) + 
  scale_y_continuous(breaks = seq(0, 30, by = 5)) +  
  theme_bw() +
  theme(panel.grid = element_blank(),  
        axis.text = element_text(size = 10),  
        axis.title = element_text(size = 10, face = "bold")) +  
  labs(x = "Number of Colonies", y = "Frequency")  

figurehBaa <- ggplot(data = statsB, aes(x = No_colonies)) +
  geom_histogram(binwidth = 50, fill = "gray51",color = "black", boundary = 0) +  
  scale_x_continuous(limits = c(0, 250), breaks = seq(0, 250, by = 50)) + 
  scale_y_continuous(breaks = seq(0, 35, by = 5)) +  
  theme_bw() +
  theme(panel.grid = element_blank(),  
        axis.text = element_text(size = 10),  
        axis.title = element_text(size = 10, face = "bold")) +  
  labs(x = "Number of Colonies", y = "Frequency")  

histNoCols <- ggplot(data = Stats, aes(x = No_colonies)) +
  geom_histogram(binwidth = 50, fill = "gray51",color = "black", boundary = 0) +  
  scale_x_continuous(limits = c(0, 250), breaks = seq(0, 250, by = 50)) + 
  scale_y_continuous(breaks = seq(0, 70, by = 10)) +  
  theme_bw() +
  theme(panel.grid = element_blank(),  
        axis.text = element_text(size = 10),  
        axis.title = element_text(size = 10, face = "bold")) +  
  labs(x = "Number of Colonies", y = "Frequency")


#### Synchronicity ####
kruskal.test(prop ~ Atoll, data = Stats) # no difference 
kruskal_effsize(prop ~ Atoll, data = Stats)

kruskal.test(prop ~ Year, data = Stats) # difference 
kruskal_effsize(prop ~ Year, data = Stats)

kruskal.test(prop ~ Month, data = Stats) # difference 
kruskal_effsize(prop ~ Month, data = Stats)

cor(Stats$SST30, Stats$SST90)

# using proportion to spawn as response
# proportion: beta regression with random effects
# all data with random effects: atoll 
hist(Stats$prop)
hist(Stats$No_colonies)

m1 <- glmmTMB(prop ~ DailyWindSpeed + DailyPpt + SST30 + TideDepth + Pos_PFM +
                (1|Atoll) + (1| Year), 
              data = Stats, 
              family = beta_family())
plot(simulateResiduals(m1))
summary(m1)
plot(allEffects(m1))
r2(m1)
check_singularity(m1) # TRUE, remove atoll as random,  var < 0.001; not contributing to the model

m1 <- glmmTMB(prop ~ DailyWindSpeed + DailyPpt + SST30 + TideDepth + Pos_PFM, 
              data = Stats, 
              family = beta_family())
plot(simulateResiduals(m1))
summary(m1)
plot(allEffects(m1))
tab_model(m1)

# by atoll
statsNM <- subset(Stats, Atoll == "NorthMale")
statsB <- subset(Stats, Atoll == "Baa")

NMm1 <- glmmTMB(prop ~ DailyWindSpeed + DailyPpt + SST30 + TideDepth + Pos_PFM, 
                data = statsNM, 
                family = beta_family())
plot(simulateResiduals(NMm1))
summary(NMm1)
plot(allEffects(NMm1))
tab_model(NMm1)

Baam1 <- glmmTMB(prop ~ DailyWindSpeed + DailyPpt + SST30 + TideDepth + Pos_PFM, 
                 data = statsB, 
                 family = beta_family())
plot(simulateResiduals(Baam1)) # issues with residuals
testQuantiles(Baam1)
summary(Baam1)
plot(allEffects(Baam1))
tab_model(Baam1)

AIC(Baam1)

#### Synchroncity: species specific where N > 30 (N = days spawning recorded) ####
# humilis
humprop <- glmmTMB(prop ~ DailyWindSpeed + DailyPpt + SST30 + TideDepth + Pos_PFM, 
                   data = statsHum, 
                   family = beta_family())
plot(simulateResiduals(humprop))
check_singularity(humprop) # singular with Atoll & Year as random: removed (var < 0.001))
summary(humprop)
plot(allEffects(humprop))
tab_model(humprop)


# secale 
secprop <- glmmTMB(prop ~ DailyWindSpeed + DailyPpt + SST30 + TideDepth + Pos_PFM, 
                   data = statsSec, 
                   family = beta_family())
plot(simulateResiduals(secprop))
check_singularity(secprop) # singular with Atoll & Year as random: removed (var < 0.001))
summary(secprop) 
plot(allEffects(secprop))
tab_model(secprop)


# tenuis 
hist(statsTen$prop)
hist(statsTen$No_colonies)
tenprop <- glmmTMB(prop ~ DailyWindSpeed + DailyPpt + SST30 + TideDepth + Pos_PFM, 
                   data = statsTen, 
                   family = beta_family())
plot(simulateResiduals(tenprop))
check_singularity(tenprop) # singular with Atoll & Year as random: removed (var < 0.001)
summary(tenprop)
plot(allEffects(tenprop))
tab_model(tenprop)



#### Timing ####
# all data 
time1 <- glmmTMB(Pos_PFM ~ DailyWindSpeed*Atoll + DailyPpt*Atoll + SST30*Atoll + TideDepth*Atoll,
                 data = Stats, family = nbinom1(link = "log"))
plot(simulateResiduals(time1))
summary(time1)ci
check_singularity(time1)
plot(allEffects(time1))
tab_model(time1)
# removed Atoll & Month as random effect -> variance so close to - causing singularity issues with model
library(ggeffects)
ggpredict(time1, terms = c("SST30", "Atoll")) %>%
  plot()
ggpredict(time1, terms = c("TideDepth", "Atoll")) %>%
  plot()

plot(statsNM$Pos_PFM, statsNM$SST30)

# by species
# secale
statsSec$obs <- 1:36
timesec <- glmmTMB(Pos_PFM ~ DailyWindSpeed*Atoll + DailyPpt*Atoll + SST30*Atoll + TideDepth*Atoll,
                   data = statsSec, family = nbinom1(link = "log"))
plot(simulateResiduals(timesec)) # issues with residuals
check_singularity(timesec) # FALSE
summary(timesec)
plot(allEffects(timesec))
tab_model(timesec)

glmmTMB::diagnose(timesec) # use poisson

timesec2 <- glmmTMB(Pos_PFM ~ DailyWindSpeed*Atoll + DailyPpt*Atoll + SST30*Atoll + TideDepth*Atoll,
                    data = statsSec, family = poisson())
plot(simulateResiduals(timesec2)) # better
summary(timesec2)
check_singularity(timesec2) # FALSE
plot(allEffects(timesec2))
tab_model(timesec2)

timesec3 <- glmmTMB(Pos_PFM ~ DailyWindSpeed*Atoll + DailyPpt + SST30 + TideDepth*Atoll,
                    data = statsSec, family = poisson())
plot(simulateResiduals(timesec3)) 
summary(timesec3)
check_singularity(timesec3) # FALSE
plot(allEffects(timesec3))
tab_model(timesec3)

AIC(timesec2) 
AIC(timesec3) # best


# tenuis
statsTen$obs <- 1:41
timeten <- glmmTMB(Pos_PFM ~ DailyWindSpeed*Atoll + DailyPpt*Atoll + SST30*Atoll + TideDepth*Atoll,
                   data = statsTen, family = nbinom1(link = "log"))
plot(simulateResiduals(timeten))
summary(timeten)
check_singularity(timeten)
plot(allEffects(timeten))
tab_model(timeten)

timeten2 <- glmmTMB(Pos_PFM ~ DailyWindSpeed*Atoll + DailyPpt + TideDepth + SST30,
                    data = statsTen, family = nbinom1(link = "log"))
plot(simulateResiduals(timeten2))
summary(timeten2)
check_singularity(timeten2)
plot(allEffects(timeten2))
tab_model(timeten2)

AIC(timeten)
AIC(timeten2) # had to remove SST30 because quantile deviations detected in residuals vs predictor plots

timeten3 <- glm.nb(Pos_PFM ~ DailyWindSpeed*Atoll + DailyPpt*Atoll + SST30*Atoll + TideDepth*Atoll,
                   data = statsTen)
plot(simulateResiduals(timeten3))
summary(timeten3)
check_singularity(timeten3)
plot(allEffects(timeten3))
tab_model(timeten3)

timeten4 <- glm.nb(Pos_PFM ~ DailyWindSpeed*Atoll + DailyPpt + SST30 + TideDepth,
                   data = statsTen)
plot(simulateResiduals(timeten4))
summary(timeten4)
check_singularity(timeten4)
plot(allEffects(timeten4))
tab_model(timeten4)

AIC(timeten3)
AIC(timeten4) # use

timeten5 <- glm.nb(Pos_PFM ~ DailyWindSpeed + DailyPpt + SST30 + TideDepth,
                   data = statsTen)
plot(simulateResiduals(timeten5))
summary(timeten5)
check_singularity(timeten5)
plot(allEffects(timeten5))
tab_model(timeten5)

AIC(timeten4)
AIC(timeten5) #

# humilis
timeHum <- glmmTMB(Pos_PFM ~ DailyWindSpeed*Atoll + DailyPpt*Atoll + SST30*Atoll + TideDepth*Atoll,
                   data = statsHum, family = nbinom1(link = "log"))
plot(simulateResiduals(timeHum)) # dispersio test p = 0.024
summary(timeHum)
check_singularity(timeHum)
plot(allEffects(timeHum))
tab_model(timeHum)

timeHum2 <- glmmTMB(Pos_PFM ~ DailyWindSpeed*Atoll + DailyPpt + SST30 + TideDepth,
                    data = statsHum, family = nbinom1(link = "log")) # warning re converging
plot(simulateResiduals(timeHum2)) 
summary(timeHum2)
check_singularity(timeHum2)
plot(allEffects(timeHum2))
tab_model(timeHum2)

glmmTMB::diagnose(timeHum2) # use poisson
statsHum$obs <- 1:30

timehum3 <- glmmTMB(Pos_PFM ~ DailyWindSpeed*Atoll + DailyPpt*Atoll + SST30*Atoll + TideDepth*Atoll,
                    data = statsHum, family = poisson())
plot(simulateResiduals(timehum3)) # dispersion issue! 
summary(timehum3)
check_singularity(timehum3)
plot(allEffects(timehum3))
tab_model(timehum3)

timeHum4<- glmmTMB(Pos_PFM ~ DailyWindSpeed*Atoll + DailyPpt*Atoll + SST30*Atoll + TideDepth*Atoll +
                     (1|obs),
                   data = statsHum, family = poisson()) 
plot(simulateResiduals(timeHum4)) 
summary(timeHum4)
check_singularity(timeHum4)
plot(allEffects(timeHum4))
tab_model(timeHum4)

AIC(timehum3) # use: drop interactions p>0.1
AIC(timeHum4)

timehum5 <- glmmTMB(Pos_PFM ~ DailyWindSpeed*Atoll + DailyPpt + SST30 + TideDepth,
                    data = statsHum, family = poisson())
plot(simulateResiduals(timehum5)) 
summary(timehum5)
check_singularity(timehum5)
plot(allEffects(timehum5))
tab_model(timehum5)

AIC(timehum3)
AIC(timehum5) # use 


# final models - put here to access 
summary(time1) # all data
plot(allEffects(time1))
tab_model(time1)

summary(timesec3) # secale
plot(allEffects(timesec3))
tab_model(timesec3)

summary(timeten5) # tenuis
plot(allEffects(timeten5))
tab_model(timeten5)

summary(timehum5) # humilis
plot(allEffects(timehum5))
tab_model(timehum5)