###Read Me###

#The fishing zones (1-17; Figure 2) consist of 1° x 1° quadrants that were used to identify the areas where red grouper fishing takes place and where fishing effort is concentrated in the southeastern Gulf of Mexico. These zones have been proposed by INAPESCA (Instituto Nacional de Pesca), and have been used to effectively delimit the areas of interest in the study.

#Gear abbreviations refer to the following meanings:Bicycle (B),longline (LL), Dinghy and shortline (DS), Combined (C).

###Data### 

(Sif<-read.csv("Semi-industrial fleet.csv"))
names(Sif)

#Class of data

Sif$Year <- factor(Sif$Year)
Sif$Zone<- factor(Sif$Zone,levels=c("1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","17","1-A","1-B","2-A"))
Sif$Month<- factor(Sif$Month,levels=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"))
Sif$Crew.size<-as.factor(Sif$Crew.size)
Sif$Fishing.gear<- factor(Sif$Fishing.gear,levels=c("LL","DS","B","CFG"))
Sif$Navigation.days<- factor(Sif$Navigation.days,levels=c("2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","17","18","19","20","21","22","23","24","25","26","27","28","29","30","31","32","33"))
   
#CPUE distribution
library(fitdistrplus) 

#Supplemental Figure S1 A and B
descdist((Sif$CPUE), discrete = FALSE)
descdist(log(Sif$CPUE), boot = 1000)



### GVIF for collinearity ###
library(car)
 
mod<- lm(CPUE ~ Year + Zone + Navigation.days + Sea.surface.temperature + Crew.size + Depth + Month + Fishing.gear, data = Sif)
(GVIF<-vif(mod))
       
    
###GAM###
library(mgcv)
options(na.action=na.omit)

gamSif<-gam(log(CPUE) ~ Year + Zone + Month +  Fishing.gear+ Crew.size +Navigation.days + s(Sea.surface.temperature)  + s(Depth), family =gaussian,data =Sif)
summary(gamSif) 


### Multimodel inference ###
library(MuMIn)
options(na.action = "na.fail")

ddSif<- dredge(gamSif,extra=c("adjR^2","R^2","deviance", F = function(x)summary(x)$fstatistic[[1]]))
head(ddSif)

#deviance explained 
mod0Sif<-gam(log(CPUE) ~ 1, family =gaussian,data=Sif)
(ddSif$Devianza.explained<-((deviance(mod0Sif)-ddSif$deviance)/deviance(mod0Sif))*100 )
head(ddSif)
 
bestmodelSif<-get.models(ddSif,1)[[1]]
(anova.gam(bestmodelSif))

#Histogram  supplemenal Figure S1 C
library(moments)

res <- data.frame(res = residuals(bestmodelSif))
hist(res$res, col = "gray", border = "gray5", prob = TRUE, ylim = c(0, 1), las = 1, main = "", xlab = "Residuals (GAM)", ylab = "Density",xlim=c(-4,4))
lines(density(res$res, adjust=2), lwd = 2) 
box(bty = "L")  
text(-3, 0.9, paste("Skewness =", round(skewness(res$res), 3)))  
text(-3, 0.85, paste("Kurtosis =", round(kurtosis(res$res), 3)))


## GAM figure ##
library(visreg) 
library(ggplot2)
library(ggpubr)
library(gridExtra)

   
v1 <- visreg(bestmodelSif,"Year", line=list(col="black"), fill = list(fill ="lightsteelblue1"),ylab = "", partial = FALSE, gg = TRUE,scale="response") + xlab("Year")+ylim(3.5,5) + theme_classic() + theme(axis.text.x = element_text(colour = "black", size = 8, angle = 90,vjust = 0.5, hjust=1), axis.text.y = element_text(colour = "black", size = 10))

v2 <- visreg(bestmodelSif,"Fishing.gear", line=list(col="black"),fill = list(fill ="lightsteelblue1"), ylab = "", partial = FALSE, gg = TRUE,scale="response") + xlab("Fishing gear") +ylim(3,5) + theme_classic() + theme(axis.text.x = element_text(colour = "black", size = 8,vjust = 0.5, hjust=1), axis.text.y = element_text(colour = "black", size = 10)) 

v3 <- visreg(bestmodelSif,"Month", line=list(col="black"),fill = list(fill ="lightsteelblue1"), ylab = "", partial = FALSE, gg = TRUE,scale="response") + xlab("Month") +ylim(4,4.8) +  theme_classic() + theme(axis.text.x = element_text(colour = "black", size = 8, angle = 90,vjust = 0.5, hjust=1), axis.text.y = element_text(colour = "black", size = 10))

v4 <- visreg(bestmodelSif,"Navigation.days", line=list(col="black"), fill = list(fill ="lightsteelblue1"),ylab = "", partial = FALSE, gg = TRUE,scale="response") + xlab("Navigation days") + ylim(3,8.5)+ theme_classic() + theme(axis.text.x = element_text(colour = "black", size = 7, angle = 90,vjust = 0.5, hjust=1), axis.text.y = element_text(colour = "black", size = 10))
 
v5 <- visreg(bestmodelSif,"Depth", line=list(col="black"),fill = list(fill ="lightsteelblue1"), ylab = "", partial = FALSE, gg = TRUE,scale="response") + xlab("Depth (m)") +ylim(4,4.8) +  scale_x_continuous(breaks=seq(from=0,to=300,by=50))+theme_classic() + theme(axis.text.x = element_text(colour = "black", size = 8), axis.text.y = element_text(colour = "black", size = 10))

v6 <- visreg(bestmodelSif,"Sea.surface.temperature", line=list(col="black"),fill = list(fill ="lightsteelblue1"), ylab = "", partial = FALSE, gg = TRUE,scale="response") + xlab("Sea surface temperature (°C)")  + ylim(3.5,5) +theme_classic() + theme(axis.text.x = element_text(colour = "black", size = 8), axis.text.y = element_text(colour = "black", size = 10))+  scale_x_continuous(breaks=seq(from=20,to=32,by=2))

v7 <- visreg(bestmodelSif,"Crew.size", line=list(col="black"),fill = list(fill ="lightsteelblue1"), ylab = "", partial = FALSE, gg = TRUE,scale="response") + xlab("Crew size") +ylim(2.5,7) + theme_classic() + theme(axis.text.x = element_text(colour = "black", size = 8,angle=90,vjust = 0.5, hjust=1), axis.text.y = element_text(colour = "black", size = 10))

v8 <- visreg(bestmodelSif,"Zone", line=list(col="black"),fill = list(fill ="lightsteelblue1"), ylab = "", partial = FALSE, gg = TRUE,scale="response") + xlab("Zone") + ylim(3.5,5) +theme_classic() + theme(axis.text.x = element_text(colour = "black", size = 8, angle = 90,vjust = 0.5, hjust=1), axis.text.y = element_text(colour = "black", size = 10))
  
#Figure X
figure <- ggarrange(v6, v5, v1, v3, v8, v4, v7, v2, ncol = 2, nrow = 4)  
annotate_figure(figure, left = text_grob("Partial for ln (CPUE)", color = "black", rot = 90, size = 12))


###Fitted CPUE###

FittedCPUE<-exp(predict(bestmodelSif,type="response"))

   
####Comparation nominal and fitted CPUE###
     
(Bw<-read.csv("Sif.Wilcoxon.csv"))
Bw$Group<-factor(Bw$Group, levels=c("Nominal CPUE","Fitted CPUE"))
Bw$Month<- factor(Bw$Month,levels=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"))

library(ggplot2)
library(ggpubr)

#Year    

lab_coords<-c(rep(350,23))
    
YearSif<- ggboxplot(Bw, x = "Year", y = "Value",fill = "Group",outlier.shape = NA)  +
  scale_fill_manual(name="",values=c("#5773CC","lightsteelblue1" ),guide = FALSE ) + 
    coord_cartesian(ylim = c(0, 420))+
       theme_classic()+
theme(legend.title = element_text( size = 20))+ 
theme(axis.title=element_text(size=18),axis.text.x = element_text(face="bold", color="black",size=16),axis.text.y = element_text(face="bold",colour="black", size=12))+
theme(legend.text = element_text( size=14, face="bold"))+ theme(legend.position="bottom")+labs(x="Year", y=expression(paste("CPUE (",   italic("kg EFD"^-1*")")))) 

YearSif1<-YearSif + stat_compare_means(aes(group = Group), label = "p.signif", label.y = lab_coords,size=7)

#Month      

lab_coords<-c(rep(350,23))
     
MonthSif<- ggboxplot(Bw, x = "Month", y = "Value",fill = "Group",outlier.shape = NA)+
  scale_fill_manual(name="",values=c("#5773CC","lightsteelblue1" )) +
  coord_cartesian(ylim = c(0, 420))+ 
              theme_classic() +
theme(legend.title = element_text( size = 20))+ 
theme(axis.title=element_text(size=18),axis.text.x = element_text(face="bold", color="black",size=16),axis.text.y = element_text(face="bold",colour="black", size=12))+
theme(legend.text = element_text( size=14, face="bold"))+ theme(legend.position="bottom")+labs(x="Month", y=expression(paste("CPUE (",   italic("kg EFD"^-1*")")))) 

MonthSif1<-MonthSif+ stat_compare_means(aes(group = Group), label = "p.signif", label.y = lab_coords,size=7)

#Figure 
library(ggpubr)
ggarrange(YearSif1,MonthSif1, ncol = 1, nrow = 2)