rm(list=ls())

#### SBT CPUE ####
## Load packages
library(tidyverse)
library(lubridate)
library(splines)
library("lunar")
library(boot)
library(influ)
library(MASS)
library(fields)
library(maps)
library(mapdata)
library(maptools)
library(readxl)
library(dplyr)
library(ggplot2)
library(ggpubr)

### Set up directories
basedir <- "D:/working/CPUE_STD/SBT_publication/public/"
figdir <- paste0(basedir, "figures/")
stddir <- paste0(basedir, "std/")
std_dl_dir <- paste0(stddir, "dl/")
clustdir <- paste0(basedir, "clusters/")
clust_dl_dir <- paste0(clustdir, "dl/")
datadir <- paste0(basedir,"Data/")
Rdir <- paste0(basedir,"RFiles/")

source(paste0(Rdir,"/make_ccsbt_sa.r"))
source(paste0(Rdir,"/dataclean_SBT.r"))
source(paste0(Rdir,"/support_functions_1.r"))
source(paste0(Rdir,"/support_functions_2.r"))
source(paste0(Rdir,"/plot_SBT_map.r"))
source(paste0(Rdir,"/plot.map.r"))
source(paste0(Rdir,"/SBT_areas_map.r"))

## Load data
rawdat <- as.data.frame(read.csv(paste0(datadir,"SBT_KOR_LL_OP_data_20190610.csv"), stringsAsFactors = FALSE))
nm <- c("op_yr","op_mon","VESSEL_CD","VESSEL_NAME_rev","DATE","Lat01","NS","Long01","EW","hooks","floats",
              "alb","bet","bft","blm","bum","mls","oth","sbt","sfa","sha","skj","swo","yft","Total","SBT_DisRelN")
names(rawdat) <- nm

prepdat <- dataprep_KRSBT(rawdat)

## Prepare data
prepdat2 <- make_ccsbt_sa(prepdat,latvar=prepdat$lat5,lonvar=prepdat$lon5)

dat_allyrs <- dataclean_KR_SBT(prepdat2, endyr = 2019)
dat <- dat_allyrs[dat_allyrs$op_yr > 1995,]

save(prepdat,prepdat2,dat_allyrs,dat, file=paste0(datadir, "KRdat.RData"))
load(file=paste0(datadir, "KRdat.RData"))

## adjust variables
dat <- dat %>% mutate(Year = op_yr, month = op_mon)

a0 <- unique(dat[,c("VESSEL_CD", "VESSEL_NAME_rev", "vessid")])
ch.vid <- dat %>% group_by(VESSEL_CD) %>% summarise(cnt = n())
ch.vid <- data.frame(ch.vid)
VESSEL_CD <- sort(unique(dat[,c("VESSEL_CD")]))
kid <- paste("K", seq(1,length(VESSEL_CD),1), sep="")
id <- data.frame(cbind(VESSEL_CD,kid))

a1 <- dat %>% dplyr::select(c(-vessid,-tripidmon)) 
a1 <- merge(a1, id, by='VESSEL_CD', all = T)
dat <- a1 %>% rename(vessid = kid)
dat$tripidmon <- paste(dat$vessid,dat$op_yr,dat$op_mon)

### make plots
setwd(figdir)

## Catch by year
cd2 <- read.csv(paste0(datadir, "KR_TotalCatchByFleet_1991-2018.csv"),header=T,stringsAsFactors=F)
catch.yr <- cd2 %>% mutate(catch = as.numeric(gsub(",", "", catch)))

windows()
ggplot(data = catch.yr, aes(x = year, y = catch)) +
  geom_bar(stat = 'identity') +
  xlab("Year") + ylab("Catch (ton)") +
  scale_x_continuous(limits = c(1990,2019), breaks = seq(1991,2018,5), minor_breaks=seq(1991,2018,5), expand = c(0,0)) +
  scale_y_continuous(limits = c(0,2010), breaks = seq(0,2000,500), expand = c(0,0), labels = c("0", "500", "1,000", "1,500", "2,000")) +
  theme_bw() +
  theme(text = element_text(size=16))
ggsave(filename = "Figure 1.png", dpi = 300)
graphics.off()

## maps of fishing effort
x.lab <- c("0", "50E", "100E", "150E")
y.lab <- c("50S", "30S", "10S")

windows(height=15,width=17); par(mfrow=c(3,2),mar=c(5,5,2,1.5))
bs <- tapply(dat$hooks,list(dat$lonfac,dat$latfac),length)
bs[is.na(bs)] <- 1
for(dec in seq(1995,2015,5)) {
  a <- dat[dat$Year >= dec & dat$Year < dec+5,]
  x <- tapply(a$hooks,list(a$lonfac,a$latfac),sum)
  image(as.numeric(dimnames(bs)[[1]]),as.numeric(dimnames(bs)[[2]]),bs,col="lightblue",xlab="Longitude",ylab="Latitude",xlim=c(-20,190),ylim=c(-55,0),main=paste(dec),cex.lab=1.4,axes = F,cex.main = 1.4)
  image(as.numeric(dimnames(x)[[1]]),as.numeric(dimnames(x)[[2]]),x,add=T,breaks=exp(seq(log(min(x,na.rm=T)),log(max(x*1.01,na.rm=T)),length.out=13)),col=rev(heat.colors(12)))
  contour(as.numeric(dimnames(x)[[1]]),as.numeric(dimnames(x)[[2]]),x,add=T,
          levels=c(10,10^2,10^3,10^4,10^5,10^6,10^7,10^8,10^9,10^10,10^11,10^12))
  axis(1, at=seq(0,150,50), labels = x.lab, cex.axis=1.1)
  axis(2, at=seq(-50,-10,20), labels = y.lab, cex.axis=1.1)
  box()
  maps::map("world",add=T, interior=F,fill=T)
  add_SBT_areas_map(ulwd=1)
}
savePlot(filename="Figure 2", type = "png")
graphics.off()


## monthly effort
len <- max(dat$Year) - min(dat$Year) + 1
a1 <- dat %>% group_by(sarea, month) %>% summarise(mean.hooks = sum(hooks/1000)/len) %>%
  filter(sarea == 2 | sarea == 8 | sarea == 9 | sarea == 14) %>%
  mutate(SA = as.factor(paste("SA", sarea, sep=" ")))

windows()
a1 %>% ggplot() +
  geom_line(aes(month, mean.hooks, linetype = SA), size = 1) +
  geom_point(aes(month, mean.hooks, shape = SA), size = 3) +
  scale_shape_manual(values = c(0,1,2,3)) + 
  xlab("Month") + ylab("Mean anuual efforts(X1,000 hooks)") +
  scale_x_continuous(limits = c(1,12), breaks = seq(1,12,1), minor_breaks=seq(1,12,2)) +
  scale_y_continuous(limits = c(0,550), breaks = seq(0,550,100), expand = c(0,0)) +
  # scale_linetype_discrete(name = "Statistical area") +
  theme_bw() +
  theme(text = element_text(size = 15)) +
  theme(legend.title = element_blank()) +
  theme(legend.position = c(0.913, 0.89), legend.text = element_text(size = 12), legend.background = element_rect(color = "black"))
ggsave(filename = "Figure 3.png", dpi = 300)
graphics.off()

## HBF
a0 <- dat %>% mutate(SA = ifelse(sarea == 8 | sarea == 9, "SA 8-9", "Others")) %>%
  mutate(SA = as.factor(SA))
a1 <- a0 %>% group_by(hbf, SA) %>% summarise(freq = n())
a1 <- data.frame(a1)
a2 <- levels(a1$SA)
label <- c(a2[2], a2[1])

windows()
ggplot(data = a1, aes(x = hbf, y = freq, fill = SA)) +
  geom_bar(stat = 'identity') +
  xlab("Number of hooks between floats (HBF)") + ylab("Frequency") +
  scale_fill_manual(breaks=c("SA 8-9", "Others"), values = c("SA 8-9" = "black", "Others" = "grey")) +
  scale_x_continuous(limits = c(4,26), breaks = seq(5,25,5), minor_breaks=seq(5,25,5), expand = c(0,0)) +
  scale_y_continuous(limits = c(0,15500), breaks = seq(0,15000,3000), expand = c(0,0), labels = c("0", "3,000", "6,000", "9,000", "12,000", "15,000")) +
  theme_bw() +
  theme(text = element_text(size=16)) +
  theme(legend.title = element_blank(), legend.position = c(0.905,0.934), legend.text = element_text(size = 12), legend.background = element_rect(color = "black"))
ggsave(filename = "Figure 4.png", dpi = 300)
graphics.off()

## CPUE by area and yrqtr
regYord <- c(14,2,9,8)
windows(height=12,width=14); par(mfrow=c(2,2),mar=c(5,4.5,2,1))
for (r in regYord) {
  llv <- dat[dat$sarea==r,]
  eff <- tapply(llv$hooks,llv$yrqtr,sum)
  yft <- tapply(llv$yft,llv$yrqtr,sum)
  bet <- tapply(llv$bet,llv$yrqtr,sum)
  alb <- tapply(llv$alb,llv$yrqtr,sum)
  sbt <- tapply(llv$sbt,llv$yrqtr,sum)
  yft <- 100*yft/eff
  bet <- 100*bet/eff
  alb <- 100*alb/eff
  sbt <- 100*sbt/eff
  if(is.na(r)) plot(1:2,1:2,type="n",axes=F,xlab="",ylab="") else {
    maxy <- max(c(sbt,bet,yft,alb))
    plot(names(yft),yft+1e-5,ylim=c(1e-5,maxy),xlab="Year",ylab="Catch per hundred hooks",main=paste("SA",r),log="y",xlim=range(dat$yrqtr), cex.lab = 1.2)
    points(names(bet),bet+1e-5,col=2,pch=2)
    points(names(alb),alb+1e-5,col=3,pch=3)
    points(names(sbt),sbt+1e-5,col=4,pch=4)
    if(r==14) legend(x=2015, y=0.007,legend=c("SBT","BET","YFT","ALB"),col=c(4,2,1,3),pch=c(4,2,1,3), cex = 0.9)
  }
}
savePlot(filename="Figure S1", type="png")
graphics.off()

##Proportion sets with zero catches by area and yrqtr
regYord <- c(14,2,9,8)
windows(height=12,width=14); par(mfrow=c(2,2),mar=c(5,4.5,2,1))
for (r in regYord) {
  llv <- dat[dat$sarea==r,]
  ay <- tapply(llv$yft==0,llv$yrqtr,sum)
  a <- tapply(llv$hooks,llv$yrqtr,length)
  ab <- tapply(llv$bet==0,llv$yrqtr,sum)
  aalb <- tapply(llv$alb==0,llv$yrqtr,sum)
  asbt <- tapply(llv$sbt==0,llv$yrqtr,sum)
  ay <- 100*ay/a
  ab <- 100*ab/a
  aalb <- 100*aalb/a
  asbt <- 100*asbt/a
  if(is.na(r)) plot(1:2,1:2,type="n",axes=F,xlab="",ylab="") else {
    maxy <- max(c(ay,ab,asbt,aalb))
    plot(names(ay),ay,xlim=range(dat$yrqtr),ylim=c(0,maxy),xlab="Year",ylab="Proportion of zero catches",main=paste("SA",r),las=1, cex.lab=1.2)
    points(names(ab),ab,col=2,pch=2)
    points(names(aalb),aalb,col=3,pch=3)
    points(names(asbt),asbt,col=4,pch=4)
    if(r==14) legend(x=2015, y=80,legend=c("SBT","BET","YFT","ALB"),col=c(4,2,1,3),pch=c(4,2,1,3), cex=0.9)
  }
}
savePlot(filename="Figure S2",type="png")
graphics.off()

## proportion of SBT & ALB
a5y <- aggregate(cbind(bet,yft,alb,swo,blm,bum,sbt,sfa,mls,sha,oth,Total,hooks) ~ lonx + lat + eval(5*floor((Year)/5)),data=dat,FUN=sum)

am <- aggregate(cbind(bet,yft,alb,swo,blm,bum,sbt,sfa,mls,sha,oth,Total,hooks) ~ lonx + lat + month,data=dat[dat$Year > 2005,],FUN=sum)

names(a5y)[3] <- "decade"

brk=seq(0,1,.1)

windows(height=15,width=17); par(mfrow=c(3,2),mar=c(5,5,2,1.5))
for(d in seq(1995,2015,5)) plot_catchmapx(indat=a5y,vbl=a5y$sbt/(a5y$Total),dcd=d,latlim=c(-50,-30),lonlim=c(-20,140),brk2=seq(0,1,.05),x.lab=c("20W","20E","60E","100E","140E"),y.lab=c("50S","40S","30S"))
image.plot(as.numeric(rownames(am)),as.numeric(colnames(am)),am,add=T,col=rev(heat.colors(length(brk)-1)),breaks=brk,legend.only=T,smallplot=c(.95,.98,.32,.87))
savePlot(filename = "Figure 5A", type = "png")

windows(height=15,width=17); par(mfrow=c(3,2),mar=c(5,5,2,1.5))
for(d in seq(1995,2015,5)) plot_catchmapx(indat=a5y,vbl=a5y$alb/(a5y$Total),dcd=d,latlim=c(-50,-30),lonlim=c(-20,140),brk2=seq(0,1,.05),x.lab=c("20W","20E","60E","100E","140E"),y.lab=c("50S","40S","30S"))
image.plot(as.numeric(rownames(am)),as.numeric(colnames(am)),am,add=T,col=rev(heat.colors(length(brk)-1)),breaks=brk,legend.only=T,smallplot=c(.95,.98,.32,.87))
savePlot(filename = "Figure 5B", type = "png")

graphics.off()

## Time distribution of vessels
vy <- list()
windows(height=7,width=12); par(mfrow=c(1,2),mar=c(4,4,2,1))
llv <- dat
vess <- as.numeric(unique(as.factor(llv$vessid)))
minyr <- maxyr <- vess
vvv <- llv[order(llv$vessid,llv$Year),]
minyr <- vvv[match(vess,as.numeric(vvv$vessid)),"Year"]
vvv <- llv[order(llv$vessid,-llv$Year),]
maxyr <- vvv[match(vess,as.numeric(vvv$vessid)),"Year"]
vessyrs <- data.frame(vess=vess,minyr=as.numeric(minyr),maxyr=as.numeric(maxyr),stringsAsFactors=F)
vessyrs <- vessyrs[order(-floor(vessyrs$minyr),-floor(vessyrs$maxyr)),]
vy[[1]] <- vessyrs

plot(1:length(vess),1:length(vess),xlim=c(1995,2019),type="n",xlab="Year",ylab="Vessel",main="Participation by arrival", cex.main=1.5, cex.lab = 1.4, cex.axis=1.3,las=1)
for (i in 1:length(vess)) {
  lines(c(floor(vessyrs[i,]$minyr),floor(vessyrs[i,]$maxyr)),c(i,i),lwd=2,type="b")
}

vessyrs <- vy[[1]]
vessyrs <- vessyrs[order(-floor(vessyrs$maxyr)),]
llv <- dat
vess <- as.numeric(unique(as.factor(llv$vessid)))

plot(1:length(vess),1:length(vess),xlim=c(1995,2019),type="n",xlab="Year",ylab="Vessel",main="Participation by departure", cex.main=1.5, cex.lab = 1.4, cex.axis=1.3, las=1)
for (i in 1:length(vess)) {
  lines(c(floor(vessyrs[i,]$minyr),floor(vessyrs[i,]$maxyr)),c(i,i),lwd=2,type="b")
}
savePlot(filename = "Figure S3", type="png")
graphics.off()


## Number of cells by area
lu <- function(x) {
  length(unique(x))
}

areas <- c(1,2,7,8,9,11,13,14,15)
fromage <- dat[dat$sarea %in% areas & dat$Year > 1995,]
a <- with(fromage,tapply(latlong,list(sarea,Year),lu));a[is.na(a)] <- 0
an <- with(fromage,tapply(latlong,list(sarea,Year),length));an[is.na(an)] <- 0
b <- with(fromage,tapply(latlong,Year,lu))
bn <- with(fromage,tapply(latlong,list(Year),length))
yrs <- as.numeric(colnames(a))

fromage$latlong1 <- paste(fromage$lat,fromage$lon)
a.2 <- with(fromage,tapply(latlong1,list(sarea,Year),lu));a.2[is.na(a.2)] <- 0
an.2 <- with(fromage,tapply(latlong1,list(sarea,Year),length));an.2[is.na(an.2)] <- 0
b.2 <- with(fromage,tapply(latlong1,Year,lu))
bn.2 <- with(fromage,tapply(latlong1,list(Year),length))
yrs.2 <- as.numeric(colnames(a.2))

fromage$llsarea <- as.factor(paste(fromage$latlong,fromage$sarea))
fromage$ll1sarea <- as.factor(paste(fromage$latlong1,fromage$sarea))
a.3 <- aggregate(ll1sarea ~ latlong+sarea+Year,lu,data=fromage)
b.3 <- tapply(a.3$ll1sarea, list(a.3$sarea, a.3$Year), mean, na.rm = TRUE)
b.3[is.na(b.3)] <- 0
b1.3 <- t(t(b.3)/apply(b.3,2,sum))
b3.3 <- tapply(a.3$ll1sarea, list(a.3$Year), mean, na.rm = TRUE)

library(reshape)
a0 <- data.frame(t(a)); names(a0) <- areas;  a0$yr <- row.names(a0)
a1 <- melt(a0, id=c("yr")); a1$yr <- as.numeric(a1$yr); names(a1) <- c("Year","sarea","no")
a1 <- a1 %>% arrange(desc(sarea))
mean.ops <- bn/b; b1 <- data.frame(mean.ops)
b1$Year <- row.names(b1); b1$Year <- as.numeric(b1$Year)

a0.2 <- data.frame(t(a.2)); names(a0.2) <- areas;  a0.2$yr <- row.names(a0.2)
a1.2 <- melt(a0.2, id=c("yr")); a1.2$yr <- as.numeric(a1.2$yr); names(a1.2) <- c("Year","sarea","no")
mean.ops.2 <- bn.2/b.2; b1.2 <- data.frame(mean.ops.2)
b1.2$Year <- row.names(b1.2); b1.2$Year <- as.numeric(b1.2$Year)

a0.3 <- data.frame(t(b1.3)); names(a0.3) <- areas;  a0.3$yr <- row.names(a0.3)
a1.3 <- melt(a0.3, id=c("yr")); a1.3$yr <- as.numeric(a1.3$yr); names(a1.3) <- c("Year","sarea","no")
b1.3 <- data.frame(b3.3); names(b1.3) <- c("mean.fished"); b1.3$Year <- as.numeric(rownames(b1.3))

detach(package:reshape)

palette <- c('#002955','#074ca1','dodgerblue','darkgreen','chartreuse3','darkolivegreen1','#7c0022','#c70125','#ff6600')

p1 <- ggplot() +
  geom_bar(data = a1, aes(x=Year, y=no, fill = sarea), stat = 'identity') +
  xlab("Year") + ylab("Number of 5x5 cells") +
  scale_fill_manual(values = palette) +
  scale_x_continuous(limits=c(1995,2019), breaks=seq(1996,2018,3), minor_breaks=seq(1996,2018,3), expand=c(0,0)) +
  theme_bw() + 
  theme(legend.title = element_blank(), legend.position = 'none', legend.text = element_text(size = 10), legend.background = element_rect(color = "black"), legend.direction = "horizontal", legend.key.size = unit(0.3, "cm")) +
  theme(text = element_text(size=11), axis.text = element_text(size=10), axis.text.y = element_text(angle=90,hjust=0.5)) +
  theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank()) +
  geom_line(data = b1, aes(x=Year, y=mean.ops*0.5), size = 1) +
  scale_y_continuous(limits = c(0,81), breaks = seq(0,80,20), expand = c(0,0), sec.axis = sec_axis(~.*2, name = "Mean number of \noperations per cell")) +
  theme(axis.title.y.right=element_text(angle=90), axis.text.y.right=element_text(hjust=0.5)) +
  guides(fill = guide_legend(nrow=1))

p2 <- ggplot() +
  geom_bar(data = a1.2, aes(x=Year, y=no, fill = sarea), stat = 'identity') +
  xlab("Year") + ylab("Number of 1x1 cells") +
  scale_fill_manual(values = palette) +
  scale_x_continuous(limits=c(1995,2019), breaks=seq(1996,2018,3), minor_breaks=seq(1996,2018,3), expand=c(0,0)) +
  theme_bw() + 
  theme(legend.title = element_blank(), legend.position = 'none', legend.text = element_text(size = 10), legend.background = element_rect(color = "black"), legend.direction = "horizontal", legend.key.size = unit(0.3, "cm")) +
  theme(text = element_text(size=11), axis.text = element_text(size=10), axis.text.y = element_text(angle=90,hjust=0.5)) + 
  theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank()) +
  geom_line(data = b1.2, aes(x=Year, y=mean.ops.2*25), size = 1) +
  scale_y_continuous(limits = c(0,510), breaks = seq(0,500,100), expand = c(0,0), sec.axis = sec_axis(~.*0.04, name = "Mean number of \noperations per cell")) +
  theme(axis.title.y.right=element_text(angle=90), axis.text.y.right=element_text(hjust=0.5)) +
  guides(fill = guide_legend(nrow=1))

p3 <- ggplot() +
  geom_bar(data = a1.3, aes(x=Year, y=no, fill = sarea), stat = 'identity') +
  xlab("Year") + ylab("Proportions of 5x5 cells fished") +
  scale_fill_manual(values = palette) +
  scale_x_continuous(limits=c(1995,2019), breaks=seq(1996,2018,3), minor_breaks=seq(1996,2018,3), expand=c(0,0)) +
  theme_bw() +
  theme(legend.title = element_blank(), legend.position = 'bottom', legend.text = element_text(size = 10), legend.background = element_rect(color = "black"), legend.direction = "horizontal", legend.key.size = unit(0.3, "cm")) +
  theme(text = element_text(size=11), axis.text = element_text(size=10), axis.text.y = element_text(angle=90,hjust=0.5)) +
  theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank()) +
  geom_line(data = b1.3, aes(x=Year, y=mean.fished*1/8), size = 1) +
  scale_y_continuous(limits = c(0,1.01), breaks = seq(0,1,0.2), expand = c(0,0), sec.axis = sec_axis(~.*8, name = "Mean number 1x1 cells \nfished per 5x5 cell", breaks = seq(0,8,2))) +
  theme(axis.title.y.right=element_text(angle=90), axis.text.y.right=element_text(hjust=0.5)) +
  guides(fill = guide_legend(nrow=1))

windows()
ggarrange(p1, p2, p3, nrow = 3, common.legend = T, legend = "bottom") +
  theme(plot.margin = unit(c(0.1,0,0,0), "cm"))
ggsave(filename = "Figure 6.png", dpi=300)

graphics.off()


#### Standardization ####
### Data selection
setwd(stddir)

dat <- dat %>% dplyr::rename(ALB=alb, BET=bet, BFT=bft, BLM=blm, BUM=bum, MLS=mls, OTH=oth, SBT=sbt, SFA=sfa, SHA=sha, SKJ=skj, SWO=swo, YFT=yft)

# run through with each of these
minqtrs=1;runsp="SBT";maxqtrs=300; domonth=T
fmla.opbin <- make_formula_SBT(runsp,modtype="deltabin",addboat=F,dohbf=F,addmonth=domonth,addmoon=T)
fmla.oppos <- make_formula_SBT(runsp,modtype="deltapos",addboat=F,dohbf=F,addmonth=domonth,addmoon=T)
fmla.boatbin <- make_formula_SBT(runsp,modtype="deltabin",addboat=T,dohbf=F,addmonth=domonth,addmoon=T)
fmla.boatpos <- make_formula_SBT(runsp,modtype="deltapos",addboat=T,dohbf=F,addmonth=domonth,addmoon=T)
fmla.oplogn <- make_formula_SBT(runsp,modtype="logn",addboat=F,dohbf=F,addmonth=domonth,addmoon=T)
fmla.boatlogn <- make_formula_SBT(runsp,modtype="logn",addboat=T,dohbf=F,addmonth=domonth,addmoon=T)
fmla.opnegb <- make_formula_SBT(runsp,modtype="negbin",addboat=F,dohbf=F,addmonth=domonth,addmoon=T)
fmla.boatnegb <- make_formula_SBT(runsp,modtype="negbin",addboat=T,dohbf=F,addmonth=domonth,addmoon=T)

# Data for each region
dat8 <- dat[dat$month %in% 7:12 & dat$sarea==8 & dat$hbf %in% 9:12,]; runreg=8
dat9 <- dat[dat$month %in% 3:10 & dat$sarea==9 & dat$hbf %in% 9:12,]; runreg=9

# first run with all the data
a8 <- select_data_KRSBT(dat8,runreg=8,minqtrs=minqtrs,runsp=runsp,mt="deltabin",addcl=NA,addpca=NA,minperyr=50,minpervess=50,minperll=100)
a9 <- select_data_KRSBT(dat9,runreg=9,minqtrs=minqtrs,runsp=runsp,mt="deltabin",addcl=NA,addpca=NA,minperyr=50,minpervess=50,minperll=100)
mn8 <- with(a8,0.1* mean(SBT/hooks))
mn9 <- with(a9,0.1* mean(SBT/hooks))

for(r in c(8, 9)) {
  runreg <- r
  a <- get(paste0("a",r))
  mn <- get(paste0("mn",r))
  
  # plot nominal vs stdized, with lognormal
  nom <- aggregate(cbind(SBT,hooks) ~ yrqtr,data=a,sum)
  with(nom,plot(unfactor(yrqtr),SBT/hooks))
  windows()
  nom <- aggregate(cbind(SBT,hooks) ~ Year,data=a,sum)
  with(nom,plot(unfactor(Year),SBT/hooks,ylim=c(0,max(SBT/hooks))))
  table(a$hbf,a$Year)

  glmdat.bin <- a
  glmdat.pos <- glmdat.bin[glmdat.bin[,3]>0,]
  wtt.oppos.area   <- mk_wts_Year(glmdat.pos,wttype="area")
  wtt.opbin.equal  <- mk_wts_Year(glmdat.bin,wttype="equal")
  wtt.oplogn.area  <- mk_wts_Year(glmdat.bin,wttype="area")
  fname <- paste0("Run_KR_",runsp,"_R",runreg,"eq ",minqtrs,"-",maxqtrs,"_qtrs")
  abin <- length(unique(glmdat.bin$Year))
  apos <- length(unique(glmdat.pos$Year))

  # Lognormal conzstant
  model.oplogn.area   <- glm(fmla.oplogn,data=glmdat.bin,weights=wtt.oplogn.area,family="gaussian");gc()
  model.opboatlogn.area   <- glm(fmla.boatlogn,data=glmdat.bin,weights=wtt.oplogn.area,family="gaussian");gc()
  coefs.oplogn.area <- get.coefs(model.oplogn.area,abin)
  coefs.opboatlogn.area <- get.coefs(model.opboatlogn.area,abin)
  summarize_and_store(mod=model.oplogn.area,dat=glmdat.bin,fname=paste0("op",runreg),modtype="logn")
  summarize_and_store(mod=model.opboatlogn.area,dat=glmdat.bin,fname=paste0("boat",runreg),modtype="logn")

  # Negative binomial
  model.opnegb.area   <- glm.nb(fmla.opnegb,data=glmdat.bin,weights=wtt.oplogn.area);gc()
  model.opboatnegb.area   <- glm.nb(fmla.boatnegb,data=glmdat.bin,weights=wtt.oplogn.area);gc()
  summarize_and_store(mod=model.opnegb.area,dat=glmdat.bin,fname=paste0("op",runreg),modtype="negb")
  summarize_and_store(mod=model.opboatnegb.area,dat=glmdat.bin,fname=paste0("boat",runreg),modtype="negb")
  
  # Binomial
  model.opbin.equal  <- glm(fmla.opbin,data=glmdat.bin,family="binomial");gc()
  coefs.opbin.equal <- get.bin.coefs(model=model.opbin.equal,abin,glmdat.bin)
  model.opboatbin.equal  <- glm(fmla.boatbin,data=glmdat.bin,weights=wtt.opbin.equal,family="binomial");gc()
  summarize_and_store(mod=model.opbin.equal,dat=glmdat.bin,fname=paste0("op",runreg),modtype="deltabin")
  summarize_and_store(mod=model.opboatbin.equal,dat=glmdat.bin,fname=paste0("boat",runreg),modtype="deltabin")
  coefs.opboatbin.equal <- get.bin.coefs(model.opboatbin.equal,abin,glmdat.bin)
  plot_effects_SBT_yr(model=model.opboatbin.equal,indat=glmdat.bin,dohbf=F,domonth=T,domoon=F,mt="deltabin"); 
  savePlot(paste(fname,"bin_effects.png",sep=""),type="png")
  
  # Lognormal positive
  model.oppos.area   <- glm(fmla.oppos,data=glmdat.pos,weights=wtt.oppos.area,family="gaussian");gc()
  coefs.oppos.area <- get.coefs(model.oppos.area,apos)
  model.opboatpos.area   <- glm(fmla.boatpos,data=glmdat.pos,weights=wtt.oppos.area,family="gaussian");gc()
  coefs.opboatpos.area <- get.coefs(model.opboatpos.area,apos)
  windows(); par(mfrow=c(2,1),mar=c(4,4,3,1))                    # plot diagnostics
  plotdiags(model.opboatpos.area$residuals,ti="Model including vessel")
  savePlot(paste(fname,"_diags.png",sep=""),type="png")
  
  plot_effects_SBT_yr(model=model.opboatpos.area,indat=glmdat.pos,dohbf=F,domonth=domonth,mt="deltapos"); 
  savePlot(paste(fname,"pos_effects.png",sep=""),type="png")
  
  summarize_and_store(mod=model.oppos.area,dat=glmdat.pos,fname=paste0("op",runreg),modtype="deltapos")
  summarize_and_store(mod=model.opboatpos.area,dat=glmdat.pos,fname=paste0("boat",runreg),modtype="deltapos")
  
  setwd(std_dl_dir)
  summarize_and_store_dl(modb=model.opbin.equal, modp=model.oppos.area, dat=glmdat.bin, datpos=glmdat.pos, fname=paste0("dl_op",runreg), modlab="dellog", dohbf = FALSE, addcl = FALSE, addmon = domonth, keepd = TRUE)
  summarize_and_store_dl(modb=model.opboatbin.equal, modp=model.opboatpos.area, dat=glmdat.bin, datpos=glmdat.pos, fname=paste0("dl_opboat",runreg), modlab="dellog", dohbf = FALSE, addcl = FALSE, addmon = domonth, keepd = TRUE)

  setwd(stddir)

  # Influence plots (lognormal constant) 
  infYear <- Influence$new(model.opboatlogn.area,data=glmdat.bin)
  infYear$calc()
  graphics.off()
  infYear$summary
  windows()
  savePlot(paste0("area",runreg,"_vessid_Influence_plots"),type="png")
  infYear$cdiPlot('latlong')
  savePlot(paste0("area",runreg,"_ll5_Influence_plots"),type="png")
  infYear$cdiPlot("ns(hooks, 5)")
  savePlot(paste0("area",runreg,"_hooks_Influence plots"),type="png")
  infYear$cdiPlot("ns(month, df = 3)")
  savePlot(paste0("area",runreg,"_month_Influence plots"),type="png")
  infYear$cdiPlot("ns(moon, df = 4)")
  savePlot(paste0("area",runreg,"_moon_Influence plots"),type="png")
  graphics.off()

  setwd(std_dl_dir)
  
  # Influence plots (delta-positive)
  infYear <- Influence$new(model.opboatpos.area,data=glmdat.pos)
  infYear$calc()
  graphics.off()
  infYear$summary
  windows()
  infYear$cdiPlot('vessid')
  savePlot(paste0("dl_area",runreg,"_vessid_Influence_plots"),type="png")
  infYear$cdiPlot('latlong')
  savePlot(paste0("dl_area",runreg,"_ll5_Influence_plots"),type="png")
  infYear$cdiPlot("ns(hooks, 5)")
  savePlot(paste0("dl_area",runreg,"_hooks_Influence plots"),type="png")
  infYear$cdiPlot("ns(month, df = 3)")
  savePlot(paste0("dl_area",runreg,"_month_Influence plots"),type="png")
  infYear$cdiPlot("ns(moon, df = 4)")
  savePlot(paste0("dl_area",runreg,"_moon_Influence plots"),type="png")
  graphics.off()
  
  # Influence plots (delta-binomial)
  infYear <- Influence$new(model.opboatbin.equal,data=glmdat.bin)
  infYear$calc()
  graphics.off()
  infYear$summary
  windows()
  infYear$cdiPlot('vessid')
  savePlot(paste0("dl_area",runreg,"_vessid_Influence_plots_bin"),type="png")
  infYear$cdiPlot('latlong')
  savePlot(paste0("dl_area",runreg,"_ll5_Influence_plots_bin"),type="png")
  infYear$cdiPlot("ns(hooks, 5)")
  savePlot(paste0("dl_area",runreg,"_hooks_Influence plots_bin"),type="png")
  infYear$cdiPlot("ns(month, df = 3)")
  savePlot(paste0("dl_area",runreg,"_month_Influence plots_bin"),type="png")
  infYear$cdiPlot("ns(moon, df = 4)")
  savePlot(paste0("dl_area",runreg,"_moon_Influence plots_bin"),type="png")
  graphics.off()
  
  setwd(stddir)
  
  # Explore effects of other species  
  a <- model.opboatlogn.area$data
  res <- model.opboatlogn.area$residuals
  str(a)
  boxplot(res ~ a$YFT)
  boxplot(res ~ a$BET)
  library(mgcv)
  mody <- gam(res ~ s(a$BET, k = 5) + s(a$YFT, k = 30))
  windows(width = 10, height = 8)
  plot.gam(mody, pages = 1, xlim = c(0, 50))
  mody <- gam(res ~ te(log(a$BET+1), log(a$YFT+1), k = c(10,30)))
  windows(width = 10, height = 8)
  plot.gam(mody, pages = 1, scheme = 2)
}
graphics.off()


### Output standardization result
# Lognormal constatn
windows(width=11,height=7);par(mfrow=c(1,2))
logn.reslist <- list()
for (r in c(8,9)) {
  # Load lognormal results & model from file
  fn <- paste0(stddir,"/boat",r,"_logn_predictions.RData")
  fn2 <- paste0(stddir,"/boat",r,"_logn_model.RData")
  load(fn);  load(fn2)
  # Nominal CPUE / slected data
  nom <- with(mod$data,tapply(SBT/hooks,Year,mean))
  nom <- nom/mean(nom,na.rm=TRUE)
  nom <- data.frame(yy=as.numeric(names(nom)),ind=nom)
  mn <- with(mod$data,0.1* mean(SBT/hooks))
  # Simple lognormal model
  mod1 <- glm(log(SBT/hooks + mn) ~ Year,data=mod$data)
  newdat <- data.frame(Year = sort(unique(mod$data$Year)))
  newdat$ind2 <- exp(predict(mod1,newdata=newdat,type="response")) - mn
  yr <- as.numeric(as.character(nd$newdat$Year))
  yrs <- seq(min(yr),max(yr))
  # Lognormal model result from file
  bb <- data.frame(yr=yr,a=exp(nd$predresp$fit)-mn)
  bb$ul <- exp(nd$predresp$fit + 2 * nd$predterms$se.fit[,1]) -mn
  bb$ll <- exp(nd$predresp$fit - 2 * nd$predterms$se.fit[,1]) -mn
  bb$cv <- nd$predterms$se.fit[,1]
  bb[,2:4] <- bb[,2:4] / mean(bb$a,na.rm=T)
  b <- bb[match(yrs,yr),]
  logn.reslist[[paste0("reg",r)]] <- b
  write.csv(b, file = paste0("indices_logn_reg",r,".csv"))
  # Plot 
  plot(yrs,b$a,type="b",xlab="Year",ylab="Relative CPUE",ylim=c(0,3),lty=2,main=paste("Area",r), xlim = c(1995, 2018))
  lines(nom$yy,nom$ind,type="b",pch=16)
  nd2 <- data.frame(yr = yrs)
  nd2$pr <- newdat$ind[match(nd2$yr, newdat$Year)]
  segments(yrs,b$a,yrs,b$ul)
  segments(yrs,b$a,yrs,b$ll)
  legend("topleft", legend = c("standardized", "nominal"), col = c(1,1), pch = c(1, 16), lty = c(2,1))
  print(b)
}
graphics.off()

# Delta lognormal
windows(width=11,height=7);par(mfrow=c(1,2))
for (r in c(8,9)) {
  ndpos <- pcoefs <- NULL
  load(paste0(std_dl_dir, "/dl_opboat", r, "_pos_dellog_predictions.RData"))
  load(paste0(std_dl_dir, "/dl_opboat", r, "_dellog_indices.RData"))
  xx <- data.frame(yq=as.numeric(as.character(ndpos$newdat$Year)))
  xx$pr <- pcoefs
  xx$ln.cv <- ndpos$predterms$se.fit[,1]
  xx$ll <- exp(log(pcoefs) - 1.96*ndpos$predterms$se.fit[,1])
  xx$ul <- exp(log(pcoefs) + 1.96*ndpos$predterms$se.fit[,1])
  
  a <- data.frame(yr=seq(min(xx$yq,na.rm=T),max(xx$yq,na.rm=T),1))
  a <- cbind(yr=a$yr,xx[match(a$yr,xx$yq),2:5])
  a[,c(2,4:5)] <- a[,c(2,4:5)]/mean(a[,2],na.rm=TRUE)
  
  est <- a$pr/mean(a$pr,na.rm=T)
  upr <- a$ul/mean(a$pr,na.rm=T)
  lwr <- a$ll/mean(a$pr,na.rm=T)
  write.csv(a,file=paste0(std_dl_dir,"/dl_opboat", r, "_dellog_indices.csv"))
  
  plot(a$yr,est,xlab="Year",ylab="Relative CPUE",main=paste0(" R",r),type="l",ylim=c(0,4), xlim = c(1995, 2018))
  findna <- complete.cases(a)
  for(i in 1:(length(findna) - 1)) {
    if(findna[i] & findna[i+1]) {
      j <- i + 1
      with(a,polygon(c(yr[i:j],yr[j:i]),c(lwr[i:j],upr[j:i]),col = "grey75", border = FALSE))
    }
  }
  
  matlines(a[,1],cbind(est,lwr,upr),
           lwd=c(2,1,1),
           lty=1,
           col=c("black","red","red"))
  points(a$yr,est,type = "b")
}
graphics.off()

library(car)
for (r in c(8,9)) {
  fn <- paste0(stddir,"/boat",r,"_logn_model.RData")
  load(fn)
  tst <- drop1(mod)
  tst$dAIC <- tst$AIC - tst$AIC[1]
  sink(paste0(stddir, "Table 1_logn_R",r,"_drop1.txt"))
  print(tst)
  sink()
  }

for (r in c(8,9)) {
  fn <- paste0(std_dl_dir,"/dl_opboat",r,"_dellog moddelta.RData")
  load(fn)
  tst <- drop1(modp)
  tst$dAIC <- tst$AIC - tst$AIC[1]
  sink(paste0(std_dl_dir, "Table 1_dellog_R",r,"_drop1.txt"))
  print(tst)
  sink()
}

for (r in c(8,9)) {
  fn <- paste0(stddir,"/boat",r,"_deltabin_model.RData")
  load(fn)
  tst <- drop1(mod)
  tst$dAIC <- tst$AIC - tst$AIC[1]
  sink(paste0(std_dl_dir, "Table 1_dellog_R",r,"_drop1_bin.txt"))
  print(tst)
  sink()
}

for (r in c(8,9)) {
  fn <- paste0(stddir,"/boat",r,"_deltapos_model.RData")
  load(fn)
  tst <- drop1(mod)
  tst$dAIC <- tst$AIC - tst$AIC[1]
  sink(paste0(std_dl_dir, "Table 1_dellog_R",r,"_drop1_pos.txt"))
  print(tst)
  sink()
}


#### Cluster analysis
## Clustering
library(randomForest)
library(influ)
library("nFactors")
library(data.table)
library(tidyverse)
library(cluster)
library(fastcluster)
library(splines)
library(boot)
library(beanplot)
library(maps)
library(lubridate)
library(mapdata)
library(maptools)
library(mgcv)
library(NbClust)

setwd(clustdir)

allsp <- c("ALB","BET","BLM","BUM","MLS","OTH","SBT","SFA","SHA","SWO","YFT")

allabs <- c("Year","month","VESSEL_CD","VESSEL_NAME_rev","vessid","DATE","Lat01","NS","Long01","EW","hooks","floats","ALB","BET","BLM","BUM","MLS","OTH","SBT","SFA","SHA","SKJ","SWO","YFT","Total","SBT_DisRelN","dmy","hbf","moon","lat","lon","lonx","lat5","lon5","lonx5","sarea","yrqtr","latlong","tripidmon")

dat <- data.frame(dat)

nclY=c(1,2,3,4,5,6,7,3,3)
flag="KR"
cvn <- c("yrqtr","latlong","hooks","hbf","vessid","Total","lat","lon","lonx","lat5","lon5","lonx5","moon","Year","month")

regtype="sarea"

for(r in 8:9) {
  fnh <- paste(flag,regtype,r,sep="_")
  dataset <- clust_PCA_run(r=r,ddd=dat,allsp=allsp,allabs=allabs,regtype=regtype,ncl=nclY[r],plotPCA=F,clustid="tripidmon",allclust=F,flag=flag,fnhead=fnh,covarnames=cvn, plotfile = "tiff")
  dataset$sarea <- dataset$reg
  assign(paste0("datcl",r),dataset)
  save(dataset,file=paste0(fnh,".RData"))
}

# Standardize
minqtrs=1;runsp="SBT";maxqtrs=300; domonth=T
fmla.opbin <- make_formula_SBT(runsp,modtype="deltabin",addboat=F,dohbf=F,addmonth=domonth,addmoon=T, addcl = "hcltrp")
fmla.oppos <- make_formula_SBT(runsp,modtype="deltapos",addboat=F,dohbf=F,addmonth=domonth,addmoon=T, addcl = "hcltrp")
fmla.boatbin <- make_formula_SBT(runsp,modtype="deltabin",addboat=T,dohbf=F,addmonth=domonth,addmoon=T, addcl = "hcltrp")
fmla.boatpos <- make_formula_SBT(runsp,modtype="deltapos",addboat=T,dohbf=F,addmonth=domonth,addmoon=T, addcl = "hcltrp")
fmla.oplogn <- make_formula_SBT(runsp,modtype="logn",addboat=F,dohbf=F,addmonth=domonth,addmoon=T, addcl = "hcltrp")
fmla.boatlogn <- make_formula_SBT(runsp,modtype="logn",addboat=T,dohbf=F,addmonth=domonth,addmoon=T, addcl = "hcltrp")
fmla.opnegb <- make_formula_SBT(runsp,modtype="negbin",addboat=F,dohbf=F,addmonth=domonth,addmoon=T, addcl = "hcltrp")
fmla.boatnegb <- make_formula_SBT(runsp,modtype="negbin",addboat=T,dohbf=F,addmonth=domonth,addmoon=T, addcl = "hcltrp")

a8 <- select_data_KRSBT(datcl8,runreg=8,minqtrs=minqtrs,runsp=runsp,mt="deltabin",addcl="hcltrp",addpca=NA,minperyr=50,minpervess=50,minperll=100)
a9 <- select_data_KRSBT(datcl9,runreg=9,minqtrs=minqtrs,runsp=runsp,mt="deltabin",addcl="hcltrp",addpca=NA,minperyr=50,minpervess=50,minperll=100)
mn8 <- with(a8,0.1* mean(SBT/hooks))
mn9 <- with(a9,0.1* mean(SBT/hooks))


for(r in c(8,9)) {
  runreg <- r
  a <- get(paste0("a",r))
  mn <- get(paste0("mn",r))
  # plot nominal vs stdized, with lognormal
  nom <- aggregate(cbind(SBT,hooks) ~ yrqtr,data=a,sum)
  with(nom,plot(unfactor(yrqtr),SBT/hooks))
  windows()
  nom <- aggregate(cbind(SBT,hooks) ~ Year,data=a,sum)
  with(nom,plot(unfactor(Year),SBT/hooks,ylim=c(0,max(SBT/hooks))))
  table(a$hbf,a$Year)

  glmdat.bin <- a
  glmdat.pos <- glmdat.bin[glmdat.bin[,3]>0,]
  wtt.oppos.area   <- mk_wts_Year(glmdat.pos,wttype="area")
  wtt.opbin.equal  <- mk_wts_Year(glmdat.bin,wttype="equal")
  wtt.oplogn.area  <- mk_wts_Year(glmdat.bin,wttype="area")
  fname <- paste0("Run_KR_clust_",runsp,"_R",runreg,"eq ",minqtrs,"-",maxqtrs,"_qtrs")
  abin <- length(unique(glmdat.bin$Year))
  apos <- length(unique(glmdat.pos$Year))

  # Lognormal constant
  model.oplogn.area   <- glm(fmla.oplogn,data=glmdat.bin,weights=wtt.oplogn.area,family="gaussian");gc()
  model.opboatlogn.area   <- glm(fmla.boatlogn,data=glmdat.bin,weights=wtt.oplogn.area,family="gaussian");gc()
  coefs.oplogn.area <- get.coefs(model.oplogn.area,abin)
  coefs.opboatlogn.area <- get.coefs(model.opboatlogn.area,abin)
  summarize_and_store(mod=model.oplogn.area,dat=glmdat.bin,fname=paste0("clust_op",runreg),modtype="logn", addcl= T, plotfile = "tiff")
  summarize_and_store(mod=model.opboatlogn.area,dat=glmdat.bin,fname=paste0("clust_boat",runreg),modtype="logn", addcl = T, plotfile = "tiff")

  # Negative binomial
  model.opnegb.area   <- glm.nb(fmla.opnegb,data=glmdat.bin,weights=wtt.oplogn.area);gc()
  model.opboatnegb.area   <- glm.nb(fmla.boatnegb,data=glmdat.bin,weights=wtt.oplogn.area);gc()
  summarize_and_store(mod=model.opnegb.area,dat=glmdat.bin,fname=paste0("clust_op",runreg),modtype="negb", addcl = T, plotfile = "tiff")
  summarize_and_store(mod=model.opboatnegb.area,dat=glmdat.bin,fname=paste0("clust_boat",runreg),modtype="negb", addcl = T, plotfile = "tiff")
  
  # Binomial
  model.opbin.equal  <- glm(fmla.opbin,data=glmdat.bin,family="binomial");gc()
  coefs.opbin.equal <- get.bin.coefs(model=model.opbin.equal,abin,glmdat.bin)
  model.opboatbin.equal  <- glm(fmla.boatbin,data=glmdat.bin,weights=wtt.opbin.equal,family="binomial");gc()
  summarize_and_store(mod=model.opbin.equal,dat=glmdat.bin,fname=paste0("clust_op",runreg),modtype="deltabin", addcl = T)
  summarize_and_store(mod=model.opboatbin.equal,dat=glmdat.bin,fname=paste0("clust_boat",runreg),modtype="deltabin", addcl = T)
  coefs.opboatbin.equal <- get.bin.coefs(model.opboatbin.equal,abin,glmdat.bin)
  plot_effects_SBT_yr(model=model.opboatbin.equal,indat=glmdat.bin,dohbf=F,domonth=T,domoon=F,mt="deltabin", addcl = T); 
  savePlot(paste(fname,"bin_effects.png",sep=""),type="png")
  
  # Lognormal positive
  model.oppos.area   <- glm(fmla.oppos,data=glmdat.pos,weights=wtt.oppos.area,family="gaussian");gc()
  coefs.oppos.area <- get.coefs(model.oppos.area,apos)
  model.opboatpos.area   <- glm(fmla.boatpos,data=glmdat.pos,weights=wtt.oppos.area,family="gaussian");gc()
  coefs.opboatpos.area <- get.coefs(model.opboatpos.area,apos)
  windows(); par(mfrow=c(2,1),mar=c(4,4,3,1))                           # plot diagnostics
  plotdiags(model.opboatpos.area$residuals,ti="Model including vessel")
  savePlot(paste(fname,"_diags.png",sep=""),type="png")
  savePlot(paste(fname,"_diags.tif",sep=""),type="tiff")
  plot_effects_SBT_yr(model=model.opboatpos.area,indat=glmdat.pos,dohbf=F,domonth=domonth,mt="deltapos", addcl = T); 
  savePlot(paste(fname,"pos_effects.png",sep=""),type="png")
  mn <- with(glmdat.bin,0.1* mean(get(runsp)/hooks))
  summarize_and_store(mod=model.oppos.area,dat=glmdat.pos,fname=paste0("clust_op",runreg),modtype="deltapos", addcl = T)
  summarize_and_store(mod=model.opboatpos.area,dat=glmdat.pos,fname=paste0("clust_boat",runreg),modtype="deltapos", addcl = T)
  
  setwd(clust_dl_dir)
  
  summarize_and_store_dl(modb=model.opbin.equal, modp=model.oppos.area, dat=glmdat.bin, datpos=glmdat.pos, fname=paste0("clust_op",runreg), modlab="dellog", dohbf = FALSE, addcl = TRUE, addmon = domonth, keepd = TRUE)
  summarize_and_store_dl(modb=model.opboatbin.equal, modp=model.opboatpos.area, dat=glmdat.bin, datpos=glmdat.pos, fname=paste0("clust_opboat",runreg), modlab="dellog", dohbf = FALSE, addcl = TRUE, addmon = domonth, keepd = TRUE)

  setwd(clustdir)
  
  # Influence plots (lognormal constant)
  infYear <- Influence$new(model.opboatlogn.area,data=glmdat.bin)
  infYear$calc()
  graphics.off()
  infYear$summary
  windows()
  infYear$cdiPlot('vessid')
  savePlot(paste0("area",runreg,"_vessid_Influence_plots"),type="png")
  infYear$cdiPlot('latlong')
  savePlot(paste0("area",runreg,"_ll5_Influence_plots"),type="png")
  infYear$cdiPlot("ns(hooks, 5)")
  savePlot(paste0("area",runreg,"_hooks_Influence plots"),type="png")
  infYear$cdiPlot("ns(month, df = 3)")
  savePlot(paste0("area",runreg,"_month_Influence plots"),type="png")
  infYear$cdiPlot("ns(moon, df = 4)")
  savePlot(paste0("area",runreg,"_moon_Influence plots"),type="png")
  infYear$cdiPlot("clust")
  savePlot(paste0("area",runreg,"_clust_Influence plots"),type="png")
  graphics.off()
  
  setwd(clust_dl_dir)
  
  # Influence plots (delta-positive)
  infYear <- Influence$new(model.opboatpos.area,data=glmdat.pos)
  infYear$calc()
  graphics.off()
  infYear$summary
  windows()
  infYear$cdiPlot('vessid')
  savePlot(paste0("dl_area",runreg,"_vessid_Influence_plots"),type="png")
  infYear$cdiPlot('latlong')
  savePlot(paste0("dl_area",runreg,"_ll5_Influence_plots"),type="png")
  infYear$cdiPlot("ns(hooks, 5)")
  savePlot(paste0("dl_area",runreg,"_hooks_Influence plots"),type="png")
  infYear$cdiPlot("ns(month, df = 3)")
  savePlot(paste0("dl_area",runreg,"_month_Influence plots"),type="png")
  infYear$cdiPlot("ns(moon, df = 4)")
  savePlot(paste0("dl_area",runreg,"_moon_Influence plots"),type="png")
  infYear$cdiPlot("clust")
  savePlot(paste0("dl_area",runreg,"_clust_Influence plots"),type="png")
  graphics.off()
  
  # Influence plots (delta-binomial)
  infYear <- Influence$new(model.opboatbin.equal,data=glmdat.bin)
  infYear$calc()
  graphics.off()
  infYear$summary
  windows()
  infYear$cdiPlot('vessid')
  savePlot(paste0("dl_area",runreg,"_vessid_Influence_plots_bin"),type="png")
  infYear$cdiPlot('latlong')
  savePlot(paste0("dl_area",runreg,"_ll5_Influence_plots_bin"),type="png")
  infYear$cdiPlot("ns(hooks, 5)")
  savePlot(paste0("dl_area",runreg,"_hooks_Influence plots_bin"),type="png")
  infYear$cdiPlot("ns(month, df = 3)")
  savePlot(paste0("dl_area",runreg,"_month_Influence plots_bin"),type="png")
  infYear$cdiPlot("ns(moon, df = 4)")
  savePlot(paste0("dl_area",runreg,"_moon_Influence plots_bin"),type="png")
  infYear$cdiPlot("clust")
  savePlot(paste0("dl_area",runreg,"_clust_Influence plots_bin"),type="png")
  graphics.off()
  
  setwd(clustdir)
  
}

### Output standardization result
# Lognormal constant
windows(width=11,height=7);par(mfrow=c(1,2))
logn.reslist <- list()
for (r in c(8,9)) {
  # Load lognormal results & model from file
  fn <- paste0(clustdir,"/clust_boat",r,"_logn_predictions.RData")
  fn2 <- paste0(clustdir,"/clust_boat",r,"_logn_model.RData")
  load(fn);  load(fn2)
  # Nominal CPUE
  nom <- with(mod$data,tapply(SBT/hooks,Year,mean))
  nom <- nom/mean(nom,na.rm=TRUE)
  nom <- data.frame(yy=as.numeric(names(nom)),ind=nom)
  mn <- with(mod$data,0.1* mean(SBT/hooks))
  # Simple lognormal model
  mod1 <- glm(log(SBT/hooks + mn) ~ Year,data=mod$data)
  newdat <- data.frame(Year = sort(unique(mod$data$Year)))
  newdat$ind2 <- exp(predict(mod1,newdata=newdat,type="response")) - mn
  yr <- as.numeric(as.character(nd$newdat$Year))
  yrs <- seq(min(yr),max(yr))
  # Lognormal model result from file
  bb <- data.frame(yr=yr,a=exp(nd$predresp$fit)-mn)
  bb$ul <- exp(nd$predresp$fit + 2 * nd$predterms$se.fit[,1]) -mn
  bb$ll <- exp(nd$predresp$fit - 2 * nd$predterms$se.fit[,1]) -mn
  bb$cv <- nd$predterms$se.fit[,1]
  bb[,2:4] <- bb[,2:4] / mean(bb$a,na.rm=T)
  b <- bb[match(yrs,yr),]
  logn.reslist[[paste0("reg",r)]] <- b
  write.csv(b, file = paste0("indices_logn_reg",r,".csv"))
  # Plot 
  plot(yrs,b$a,type="b",xlab="Year",ylab="Relative CPUE",ylim=c(0,3),lty=2,main=paste("Area",r), xlim = c(1995, 2018))
  lines(nom$yy,nom$ind,type="b",pch=16)
  nd2 <- data.frame(yr = yrs)
  nd2$pr <- newdat$ind[match(nd2$yr, newdat$Year)]
  segments(yrs,b$a,yrs,b$ul)
  segments(yrs,b$a,yrs,b$ll)
  legend("topleft", legend = c("standardized", "nominal"), col = c(1,1), pch = c(1, 16), lty = c(2,1))
  print(b)
}
graphics.off()

# Delta lognormal
windows(width=11,height=7);par(mfrow=c(1,2))
for (r in c(8,9)) {
  ndpos <- pcoefs <- NULL
  load(paste0(clust_dl_dir, "/clust_opboat", r, "_pos_dellog_predictions.RData"))
  load(paste0(clust_dl_dir, "/clust_opboat", r, "_dellog_indices.RData"))
  xx <- data.frame(yq=as.numeric(as.character(ndpos$newdat$Year)))
  xx$pr <- pcoefs
  xx$ln.cv <- ndpos$predterms$se.fit[,1]
  xx$ll <- exp(log(pcoefs) - 1.96*ndpos$predterms$se.fit[,1])
  xx$ul <- exp(log(pcoefs) + 1.96*ndpos$predterms$se.fit[,1])
  
  a <- data.frame(yr=seq(min(xx$yq,na.rm=T),max(xx$yq,na.rm=T),1))
  a <- cbind(yr=a$yr,xx[match(a$yr,xx$yq),2:5])
  a[,c(2,4:5)] <- a[,c(2,4:5)]/mean(a[,2],na.rm=TRUE)
  
  est <- a$pr/mean(a$pr,na.rm=T)
  upr <- a$ul/mean(a$pr,na.rm=T)
  lwr <- a$ll/mean(a$pr,na.rm=T)
  write.csv(a,file=paste0(clust_dl_dir,"/clust_opboat", r, "_dellog_indices.csv"))
  
  plot(a$yr,est,xlab="Year",ylab="Relative CPUE",main=paste0(" R",r),type="l",ylim=c(0,4), xlim = c(1995, 2018))
  findna <- complete.cases(a)
  for(i in 1:(length(findna) - 1)) {
    if(findna[i] & findna[i+1]) {
      j <- i + 1
      with(a,polygon(c(yr[i:j],yr[j:i]),c(lwr[i:j],upr[j:i]),col = "grey75", border = FALSE))
    }
  }
  
  matlines(a[,1],cbind(est,lwr,upr),
           lwd=c(2,1,1),
           lty=1,
           col=c("black","red","red"))
  points(a$yr,est,type = "b")
}
graphics.off()

library(car)
for (r in c(8,9)) {
  fn <- paste0(clustdir,"/clust_boat",r,"_logn_model.RData")
  load(fn)
  tst <- drop1(mod)
  tst$dAIC <- tst$AIC - tst$AIC[1]
  sink(paste0(clustdir, "Table 1_clust_logn_R",r,"_drop1.txt"))
  print(tst)
  sink()
}

for (r in c(8,9)) {
  fn <- paste0(clust_dl_dir,"/clust_opboat",r,"_dellog moddelta.RData")
  load(fn)
  tst <- drop1(modp)
  tst$dAIC <- tst$AIC - tst$AIC[1]
  sink(paste0(clust_dl_dir, "Table 1_clust_dellog_R",r,"_drop1.txt"))
  print(tst)
  sink()
}

for (r in c(8,9)) {
  fn <- paste0(clustdir,"/clust_boat",r,"_deltabin_model.RData")
  load(fn)
  tst <- drop1(mod)
  tst$dAIC <- tst$AIC - tst$AIC[1]
  sink(paste0(clust_dl_dir, "Table 1_clust_dellog_R",r,"_drop1_bin.txt"))
  print(tst)
  sink()
}

for (r in c(8,9)) {
  fn <- paste0(clustdir,"/clust_boat",r,"_deltapos_model.RData")
  load(fn)
  tst <- drop1(mod)
  tst$dAIC <- tst$AIC - tst$AIC[1]
  sink(paste0(clust_dl_dir, "Table 1_clust_dellog_R",r,"_drop1_pos.txt"))
  print(tst)
  sink()
}


## comparison plots of standardized CPUEs
d8 <- dat %>% dplyr::filter(sarea == 8)
d8 <- d8 %>% dplyr::filter(Year != 2003 & Year != 2005)   # omit too low sample size
d9 <- dat %>% dplyr::filter(sarea == 9)
dat1 <- data.frame(rbind(d8, d9))

# Lognormal constant
windows(width=11,height=5);par(mfrow=c(1,2), mar=c(4.5,5,3,1))
for (r in c(8,9)) {
  fn <- paste0(clustdir,"/clust_boat",r,"_logn_predictions.RData")
  load(fn)
  bstd <- read.csv(file = paste0(stddir,"indices_logn_reg",r,".csv"))
  # Plot 
  yr <- as.numeric(as.character(nd$newdat$Year))
  yrs <- seq(min(yr),max(yr))
  a0 <- dat1[dat1$sarea == r, ]
  nom <- tapply(a0$SBT/a0$hooks, a0$Year, mean)
  nom1 <- data.frame(yr = as.numeric(names(nom)), a = nom / mean(nom))
  nom2 <- data.frame(yrs = yrs, a = nom1$a[match(yrs, nom1$yr)])
  bb <- data.frame(yr=yr,a=exp(nd$predresp$fit))
  bb$ul <- exp(nd$predresp$fit + 2 * nd$predterms$se.fit[,1])
  bb$ll <- exp(nd$predresp$fit - 2 * nd$predterms$se.fit[,1])
  bb$cv <- nd$predterms$se.fit[,1]
  bb[,2:4] <- bb[,2:4] / mean(bb$a,na.rm=T)
  b <- bb[match(yrs,yr),]
  plot(yrs,b$a,type="b",xlab="Year",ylab="Relative CPUE",ylim=c(0,4.1),pch = 0, lty=3,main=paste("Statistical area",r), xlim=c(1995, 2018),cex.main=1.4,cex.axis=1.2,cex.lab=1.3,las=1)
  lines(yrs,bstd$a,type="b",lty=2, pch = 2)
  arrows(yrs,b$a,yrs,b$ul, lty=2, angle = 90, length = .03)
  arrows(yrs,b$a,yrs,b$ll, lty=2, angle = 90, length = .03)
  arrows(yrs,bstd$a,yrs,bstd$ul, lty=1, angle = 90, length = .03)
  arrows(yrs,bstd$a,yrs,bstd$ll, lty=1, angle = 90, length = .03)
  lines(nom2$yr,nom2$a,type="b",pch = 1, lty = 1)
  print(b)
}
legend("topleft", legend = c("Nominal", "Data selection", "Clustering"), pch = c(1,2,0), lty = c(1,2,3), cex=1.2)
savePlot(filename=paste0(figdir, "Figure 11A"),type="png")
graphics.off()

# Delta lognormal
windows(width=11,height=5);par(mfrow=c(1,2), mar=c(4.5,5,3,1))
for (r in c(8, 9)) {
  dlclu <- read.csv(file=paste0(clust_dl_dir,"/clust_opboat", r, "_dellog_indices.csv"))
  dlstd <- read.csv(file=paste0(std_dl_dir,"/dl_opboat", r, "_dellog_indices.csv"))
  # Plot 
  yr <- as.numeric(as.character(nd$newdat$Year))
  yrs <- seq(min(yr),max(yr))
  a0 <- dat1[dat1$sarea == r, ]
  nom <- tapply(a0$SBT/a0$hooks, a0$Year, mean)
  nom1 <- data.frame(yr = as.numeric(names(nom)), a = nom / mean(nom))
  nom2 <- data.frame(yrs = yrs, a = nom1$a[match(yrs, nom1$yr)])
  plot(dlclu$yr,dlclu$pr,type="b",xlab="Year",ylab="Relative CPUE",ylim=c(0,4.1),pch = 0, lty=3,main=paste("Statistical area",r), xlim = c(1995, 2018),cex.main=1.4,cex.axis=1.2,cex.lab=1.3,las=1)
  arrows(dlclu$yr,dlclu$pr,dlclu$yr,dlclu$ll, lty=1, angle = 90, length = .03)
  arrows(dlclu$yr,dlclu$pr,dlclu$yr,dlclu$ul, lty=1, angle = 90, length = .03)
  lines(dlstd$yr,dlstd$pr,type="b",lty=2, pch = 2, lwd = 1)
  arrows(dlstd$yr,dlstd$pr,dlstd$yr,dlstd$ll, lty=1, angle = 90, length = .03)
  arrows(dlstd$yr,dlstd$pr,dlstd$yr,dlstd$ul, lty=1, angle = 90, length = .03)
  lines(nom2$yr,nom2$a,type="b",pch = 1, lty = 1)
  print(b)
}
legend("topleft", legend = c("Nominal", "Data selection", "Clustering"), pch = c(1,2,0), lty = c(1,2,3), lwd = c(1,1,1), cex=1.2)
savePlot(filename=paste0(figdir, "Figure 11B"), type="png")
graphics.off()