rm(list=ls())
if (!require("glmmTMB")) install.packages("glmmTMB"); library("glmmTMB")
if (!require("DHARMa")) install.packages("DHARMa"); library("DHARMa")
if (!require("ggeffects")) install.packages("ggeffects"); library("ggeffects")
if (!require("ggplot2")) install.packages("ggplot2"); library("ggplot2")
if (!require("performance")) install.packages("performance"); library("performance")
if (!require("gridExtra")) install.packages("gridExtra"); library("gridExtra")
if (!require("sjPlot")) install.packages("sjPlot"); library("sjPlot")
if (!require("marginaleffects")) install.packages("marginaleffects"); library("marginaleffects")


main <- read.csv("maindf.csv")


fishsimpglm <- glmmTMB(fish_simpson ~ PC1 + PC2+ HC + Protection, 
                       data = main,
                       family = beta_family(link = "logit")) # no random effects because of singularity problem
summary(fishsimpglm) # marginal sig for PC2
fishsimpsummary <- avg_slopes(fishsimpglm) # to get average marginal effects

check_collinearity(fishsimpglm)
fishsimpsim<- simulateResiduals(fittedModel = fishsimpglm)
testResiduals(fishsimpsim)
performance::r2(fishsimpglm)

#fish simpson plot
fishsimp <-ggpredict(fishsimpglm, terms = "PC2") 
fishsimpplot <- ggplot(fishsimp, aes(x=x, y = predicted))+
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
              alpha = 0.1, colour = NA) +
  geom_line(linewidth = 1)+
  geom_point(data = main,
             aes(x = PC2, y = fish_simpson),
             inherit.aes = FALSE,
             alpha = 0.4,
             colour = ifelse(main$Location == "Koh Rong Marine National Park", "blue", "darkgrey"),
             size = 2.5) +
  labs(x = "PC2 scores",
       y = "Predicted fish diversity (Simpson's 1-D)", 
       tag = "A")+
  theme_classic(base_size = 14) +
  theme(
    legend.position = c(0.85, 0.2),
    legend.background = element_rect(fill = "white")
  )


#### fish abundance 
fishabunglm <- glmmTMB(fish_total_abun ~ PC1 + PC2 + HC + Protection + (1|Location), 
                       data = main,
                       family = nbinom2()) # pc1, hc sig

summary(fishabunglm)
r2(fishabunglm)
check_collinearity(fishabunglm)

#model diagnostics
fishabunsim<- simulateResiduals(fittedModel = fishabunglm)
testResiduals(fishabunsim)

#fish abundance plot
quantile(main$HC)[2:4]
fishabundance <-ggpredict(fishabunglm, terms = c("PC1", "HC[0.308, 0.435, 0.589]"))
plot(fishabundance)
fishabunplot <- ggplot(fishabundance, aes(x=x, y = predicted, colour = group))+
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group),
              alpha = 0.1, colour = NA) +
  geom_line(linewidth = 1)+
  geom_point(data = main,
             aes(x = PC1, y = fish_total_abun),
             inherit.aes = FALSE,
             alpha = 0.4,
             colour = ifelse(main$Location == "Koh Rong Marine National Park", "blue", "darkgrey"),
             size = 2.5) +
  labs(x = "PC1 scores",
       y = "Predicted total fish abundance", 
       colour = "Hard coral cover", 
       fill = "Hard coral cover",
       tag = "A")+
  theme_classic(base_size = 14) +
  theme(
    legend.position = c(0.85, 0.8),
    legend.background = element_rect(fill = "white")
  )


### fish biomass 

fishbiomassglm <- glmmTMB(log(total_site_biomass) ~ PC1 + PC2 + HC + Protection + (1|Location), 
                       data = main,
                       family = gaussian)  
summary(fishbiomassglm)  #hard coral sig only
fishbiomasssim<- simulateResiduals(fittedModel = fishbiomassglm)
testResiduals(fishbiomasssim)
r2(fishbiomassglm)
check_collinearity(fishbiomassglm)

#fish biomass plot
fishbiomass <- ggpredict(fishbiomassglm, terms = "HC")
plot(fishbiomass)

fishbiomassplot <- ggplot(fishbiomass, aes(x=x, y = log(predicted))) + 
  geom_ribbon(aes(ymin = log(conf.low), ymax = log(conf.high)),
              alpha = 0.1, colour = NA) +
  geom_line(linewidth = 1)+
  geom_point(data = main,
             aes(x = HC, y = log(total_site_biomass)),
             inherit.aes = FALSE,
             alpha = 0.4,
             colour = ifelse(main$Location == "Koh Rong Marine National Park", "blue", "darkgrey"),
             size = 2.5) +
  labs(x = "Hard coral cover",
       y = "Predicted log site-wide fish biomass (g)")+
  theme_classic(base_size = 14)

# invert simpson diversity glmm
invertsimpglm <- glmmTMB(invert_simpson ~ PC1 + PC2 + HC + Protection, 
                       data = main,
                       family = beta_family(link = "logit"))
summary(invertsimpglm) #strong pc1 effect - increases
check_collinearity(invertsimpglm)
r2(invertsimpglm)
invertsimpsim<- simulateResiduals(fittedModel = invertsimpglm)
testResiduals(invertsimpsim) # this suggests that we have significant outliers
invertsimpout <- outliers(invertsimpsim, return = "logical") # identifies outliers
invertsimpdf <- main[!invertsimpout,] # remove outlier rows
invertsimpglm_nooutliers <- glmmTMB(invert_simpson ~ PC1 + PC2 + HC + Protection, 
                         data = invertsimpdf,
                         family = beta_family(link = "logit"))
check_collinearity(invertsimpglm_nooutliers)
summary(invertsimpglm_nooutliers)
avg_slopes(invertsimpglm_nooutliers)
r2(invertsimpglm_nooutliers)
invertsimpsim<- simulateResiduals(fittedModel = invertsimpglm_nooutliers)
testResiduals(invertsimpsim) # this looks OK now

invertsimp <-ggpredict(invertsimpglm_nooutliers, terms = "PC1")
plot(invertsimp)
invertsimpplot <- ggplot(invertsimp, aes(x=x, y = predicted))+
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
              alpha = 0.1, colour = NA) +
  geom_line(linewidth = 1)+
  geom_point(data = main,
             aes(x = PC1, y = invert_simpson),
             inherit.aes = FALSE,
             alpha = 0.4,
             colour = ifelse(main$Location == "Koh Rong Marine National Park", "blue", "darkgrey"),
             size = 2.5) +
  labs(x = "PC1 scores",
       y = "Predicted invertebrate diversity (Simpson's 1-D)",
       tag = "B")+
  theme_classic(base_size = 14)

# invert abundance glmm
invertabunglm <- glmmTMB(invert_total_abun ~ PC1 + PC2 + HC + Protection + (1|Location), 
                       data = main,
                       family = nbinom2()) 

summary(invertabunglm) #pc1 - decreases. nbinom2 same effect
r2(invertabunglm)

invertabunglmsim<- simulateResiduals(fittedModel = invertabunglm)
resinvertabun <- plotResiduals(invertabunglmsim)
testResiduals(invertabunglmsim)
invertabun<-ggpredict(invertabunglm, terms = "PC1")
plot(invertabun)

invertabunplot <- ggplot(invertabun, aes(x=x, y = predicted))+
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
              alpha = 0.1, colour = NA) +
  geom_line(linewidth = 1)+
  geom_point(data = main,
             aes(x = PC1, y = invert_total_abun),
             inherit.aes = FALSE,
             alpha = 0.4,
             colour = ifelse(main$Location == "Koh Rong Marine National Park", "blue", "darkgrey"),
             size = 2.5) +
  labs(x = "PC1 scores",
       y = "Predicted total invertebrate abundance",
       tag = "B")+
  theme_classic(base_size = 14)

grid.arrange(fishsimpplot, invertsimpplot, nrow = 1)
grid.arrange(fishabunplot, invertabunplot, nrow = 1)

tab_model(fishsimpglm, invertsimpglm_nooutliers, transform = NULL, show.stat = T) # dont transform to keep it on logit scale here. we also use AME for interpretation  
tab_model(fishabunglm, invertabunglm, transform = NULL, show.stat = T) #original coeffs
tab_model(fishabunglm, invertabunglm, show.stat = T)# using IRR for interpretation
tab_model(fishbiomassglm,show.stat = T)