###Read Me###

#Gear abbreviations refer to the following meanings:Short longline (SL), Headlines (H), Scuba diving (SD), Jimba and handlines (JH), Gillnets and handlines (GH).


###Data### 

(Ssf<-read.csv("Small-scale fleet.csv"))
names(Ssf)

#Class of data
   
Ssf$Year <- factor(Ssf$Year)
Ssf$Month<- factor(Ssf$Month,levels=c("Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"))
Ssf$Fishing.gear<- factor(Ssf$Fishing.gear,levels=c("SL","H","SD","JH","GH"))
Ssf$Place.of.origin<- factor(Ssf$Place.of.origin, levels=c("Celestun","Sisal","Chuburna","Progreso","TelchacXXPuerto","SanXXCrisanto","Dzilam deXXBravo","SanXXFelipe","RiaXXLagartos","LasXXColoradas","El Cuyo"))
levels(Ssf$Place.of.origin) <- gsub("XX", "\n", levels(Ssf$Place.of.origin))
Ssf$Crew.size<-factor(Ssf$Crew.size)

#CPUE distribution
library(fitdistrplus) 

#Supplemental Figure S2
descdist((Ssf$CPUE), boot = 1000, discrete = FALSE) 
 
 
### GVIF for collinearity ###
library(car)

mod<- lm(CPUE ~ Year + Place.of.origin + Month +  Fishing.gear+ Crew.size + Sea.surface.temperature + Depth,data =Ssf)
(GVIFScfleet<-vif(mod))
       
       
###GAM###
library(mgcv)
options(na.action=na.omit)

gamSsf<-gam(CPUE ~ Year + Place.of.origin + Month +  Fishing.gear + Crew.size + s(Sea.surface.temperature)  + s(Depth), family =Gamma(link=inverse),data =Ssf)
summary(gamSsf)


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

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

#deviance explained 
mod0Ssf<-gam(CPUE ~ 1, family =Gamma(link=inverse),data =Ssf)
(ddSsf$Devianza.explained<-((deviance(mod0Ssf)-ddSsf$deviance)/deviance(mod0Ssf))*100 )

bestmodelSsf<-get.models(ddSsf,1)[[1]]
(anova.gam(bestmodelSsf))
     

## GAM figure ##
library(visreg) 
library(ggplot2)
library(ggpubr)
library(gridExtra)
 
v1 <- visreg(bestmodelSsf,"Year", line=list(col="black"),fill = list(fill ="lightgoldenrod"),ylab = "", partial = F, gg = TRUE,scale="response") + xlab("Year")+ylim(0.5,2.5) + theme_classic() + theme(axis.text.x = element_text(colour = "black", size = 8), axis.text.y = element_text(colour = "black", size = 10))
                                                                     
v2 <- visreg(bestmodelSsf,"Fishing.gear", line=list(col="black"), fill = list(fill ="lightgoldenrod"),ylab = "", partial = FALSE, gg = TRUE,scale="response") + xlab("Fishing gear") +ylim(0,3.2) + 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)) 
                                                               
v3 <- visreg(bestmodelSsf,"Month", line=list(col="black"),fill = list(fill ="lightgoldenrod"), ylab = "", partial = FALSE, gg = TRUE,scale="response") + xlab("Month") +ylim(0.5,2.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))
   
v4 <- visreg(bestmodelSsf,"Place.of.origin", line=list(col="black"),fill = list(fill ="lightgoldenrod"), ylab = "", partial = FALSE, gg = TRUE,scale="response") + xlab("Place of origin") + ylim(0.5,3)+ 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))
                               
v5 <- visreg(bestmodelSsf,"Depth", line=list(col="black"),fill = list(fill ="lightgoldenrod"), ylab = "", partial = FALSE, gg = TRUE,scale="response") + xlab("Depth (m)") +ylim(0.5,2.5) +  scale_x_continuous(breaks=seq(from=0,to=60,by=15))+theme_classic() + theme(axis.text.x = element_text(colour = "black", size = 8), axis.text.y = element_text(colour = "black", size = 10))

v6 <- visreg(bestmodelSsf,"Sea.surface.temperature", line=list(col="black"),fill = list(fill ="lightgoldenrod"), ylab = "", partial = FALSE, gg = TRUE,scale="response") + xlab("Sea surface temperature (?C)")  + ylim(0.5,2.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=22,to=32,by=2))

v7 <- visreg(bestmodelSsf,"Crew.size", line=list(col="black"),fill = list(fill ="lightgoldenrod"), ylab = "", partial = FALSE, gg = TRUE,scale="response") + xlab(" Crew size") +ylim(0.5,2.5) + theme_classic() + theme(axis.text.x = element_text(colour = "black", size = 8), axis.text.y = element_text(colour = "black", size = 10))

#Figure X
figure <- ggarrange(v6, v5, v1, v3,v4,v7, v2, ncol = 2, nrow = 4)  
annotate_figure(figure, left = text_grob("Partial for CPUE", color = "black", rot = 90, size = 12))


###Fitted CPUE###

FittedCPUE<-predict(bestmodelSsf,type="response")

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

library(ggplot2)
library(ggpubr)

#Year

lab_coords<-c(rep(8,23))
     
YearSsf<- ggboxplot(Bw, x = "Year", y = "Value",fill = "Group",outlier.shape = NA)  +
  scale_fill_manual(name="",values=c("goldenrod4","lightgoldenrod"),guide = FALSE ) + coord_cartesian(ylim = c(0, 10))+
     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=12),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 EFH"^-1*")")))) 

YearSsf1<-YearSsf + stat_compare_means(aes(group = Group), label = "p.signif", label.y = lab_coords,size=6)

#Month

lab_coords<-c(rep(12.5,23))

MonthSsf<- ggboxplot(Bw, x = "Month", y = "Value",fill = "Group",outlier.shape = NA)  +  scale_fill_manual(name="",values=c("goldenrod4","lightgoldenrod"))+coord_cartesian(ylim = c(0, 15))+ 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=12),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 EFH"^-1*")")))) 
   
MonthSsf1<-MonthSsf+ stat_compare_means(aes(group = Group), label = "p.signif", label.y = lab_coords,size=6)

#Figure 
library(ggpubr)
ggarrange(YearSsf1,MonthSsf1, ncol = 1, nrow = 2)