#install.packages("rjags")
library(rjags)
library(ggplot2)
df <- as.data.frame(read.csv("Data_AO_YFT.csv"))#Loading data
Year <- df$Year
sc <- 10000 #Scaling
Catch <- df$Catch/sc
CPUE <- df$CPUE
f.Year <- length(Year)
CPUE_Year_initial <- 30
CPUE_Year_final <- 69
Model <-
"model {
for(i in 2:f.Year){
D[i]~dlnorm(Dtmp[i],inv.tau2)
tmp[i]<-D[i-1]+r*D[i-1]*PF[i-1]-Catch[i-1]/K
PF[i-1]<- 1-D[i-1]
Dtmp[i]<-log(max(tmp[i],0.0001))
}
for(i in CPUE_Year_initial:CPUE_Year_final){
CPUE.hat[i]<-q*B[i]
muC[i]<-log(CPUE.hat[i])
CPUE[i]~dlnorm (muC[i],inv.sigma2.C)
CPUE.est[i]~dlnorm(muC[i],inv.sigma2.C)
}
##Bratio,Fratio
for(i in 1:f.Year){
B[i]<-K*D[i]
Bratio[i]<-B[i]/BMSY
F[i]<-Catch[i]/B[i]
Fratio[i]<-F[i]/FMSY
}
MSY<-r*K/4
BMSY<-K/2
FMSY<-MSY/BMSY
##prior
r~dunif(0.01,2)
K~dunif(1,500)
Dep <- 1
logDep<-log(Dep)
D[1]~dlnorm(logDep,inv.tau2)
tau~dunif(0.00001,1)
inv.tau2 <- pow(tau,-2)
qL <- pow(10,-5)
qU <- pow(10,-1)
q~dunif(qL,qU)
sigma.C~dunif(0.00001,1)
inv.sigma2.C <- pow(sigma.C,-2)
}"
writeLines(Model, con ="JAGS.bug.txt")
Data bust be sapplied as list(). Initialvalue must be sapplied as list(). *Setting parameters which are sampled by MCMC.
#data
JAGS.data <- list(CPUE=CPUE,Catch=Catch, f.Year = f.Year,CPUE_Year_initial=CPUE_Year_initial,
CPUE_Year_final=CPUE_Year_final)
# Initial values
inits <- list(r = runif(1,0.01,2),K = runif(1,1,500),tau = runif(1,0.00001,1),q = runif(1,10^-5,10^-1),sigma.C = runif(1,0.00001,1))
parameters <- c("r","K","D","B","F","CPUE.est","BMSY","FMSY","MSY","q","tau","sigma.C","Dep")
#MCMC settings
ni <- 500000
nb <- ni/2
nt <- 100
nc <- 3
JAGS.out <- jags.model(
file = "JAGS.bug.txt",
data = JAGS.data,
inits = list(inits, inits, inits),
n.chain = 3)
update(JAGS.out, nb)
res_JAGS <- coda.samples(
JAGS.out,
parameters,
thin = nt, n.iter = ni-nb
)
table1 <- function(sample){
c(median(sample),quantile(sample,c(0.1,0.9),names = FALSE))
}#function
colnames <- c("Median","Q10","Q90","Rhat")
rownames <- c("r","K","initDep","Sig_Pro","Sig_Obs","MSY","BMSY")
#labels(res_JAGS[[1]])[[2]]
K_sample <- c(res_JAGS[[1]][,251],res_JAGS[[2]][,251],res_JAGS[[3]][,251])*sc
r_sample <- c(res_JAGS[[1]][,254],res_JAGS[[2]][,254],res_JAGS[[3]][,254])
D1_sample<- c(res_JAGS[[1]][,180],res_JAGS[[2]][,180],res_JAGS[[3]][,180])
q_sample<- c(res_JAGS[[1]][,253],res_JAGS[[2]][,253],res_JAGS[[3]][,253])
sigma.p_sample<- c(res_JAGS[[1]][,256],res_JAGS[[2]][,256],res_JAGS[[3]][,256])
sigma.o_sample<- c(res_JAGS[[1]][,255],res_JAGS[[2]][,255],res_JAGS[[3]][,255])
MSY_sample<- c(res_JAGS[[1]][,252],res_JAGS[[2]][,252],res_JAGS[[3]][,252])*sc
BMSY_sample<- c(res_JAGS[[1]][,70],res_JAGS[[2]][,70],res_JAGS[[3]][,70])*sc
table1_JAGS <- rbind(table1(r_sample),table1(K_sample),c(1,"NA","NA"),table1(sigma.p_sample),
table1(sigma.o_sample),table1(MSY_sample),table1(BMSY_sample))
Rhats <- c("NA","NA","NA","NA","NA","NA","NA")
table1_JAGS <- cbind(table1_JAGS,Rhats)
dimnames(table1_JAGS) <- list(rownames,colnames)
#write.csv(table1_JAGS,"table1-S_JAGS.csv",row.names = TRUE)
table1_JAGS
## Median Q10 Q90
## r "0.164418818981463" "0.0810278247141653" "0.292026489091953"
## K "3194534.438305" "1928272.62542323" "4584313.24984061"
## initDep "1" "NA" "NA"
## Sig_Pro "0.115521377747541" "0.0881970597866279" "0.141065819696729"
## Sig_Obs "0.0380665528393847" "0.00847793051236359" "0.0736637032697919"
## MSY "125920.426331905" "79338.2149353183" "175734.630401282"
## BMSY "1597267.2191525" "964136.312711617" "2292156.62492031"
## Rhat
## r "NA"
## K "NA"
## initDep "NA"
## Sig_Pro "NA"
## Sig_Obs "NA"
## MSY "NA"
## BMSY "NA"
colnames <- c("Year","Median_Biomass","Q10_Biomass","Q90_Biomass","Median_CPUE","Q10_CPUE","Q90_CPUE","Median_Depletion","Q10_Depletion","Q90_Depletion")
table2 <- function(matrix){
quant01<-function(x){quantile(x,0.1)}
quant09<-function(x){quantile(x,0.9)}
cbind(apply(matrix,2,median),apply(matrix,2,quant01),apply(matrix,2,quant09))
}#Define fuction
#labels(res_JAGS[[1]])[2]
B_matrix <- rbind(res_JAGS[[1]][,1:69],res_JAGS[[2]][,1:69],res_JAGS[[3]][,1:69])*sc
CPUEest_matrix <- rbind(res_JAGS[[1]][,71:110],res_JAGS[[2]][,71:110],res_JAGS[[3]][,71:110])
D_matrix <- rbind(res_JAGS[[1]][,111:179],res_JAGS[[2]][,111:179],res_JAGS[[3]][,111:179])
resCPUEest<- rbind(matrix(NA,dim(table2(B_matrix))[1]-dim(table2(CPUEest_matrix))[1],3),table2(CPUEest_matrix))
table2_JAGS <- cbind(Year,table2(B_matrix),resCPUEest,table2(D_matrix))
table2_JAGS <- rbind(colnames,table2_JAGS)
#write.csv(table2_JAGS,"table2-S_JAGS.csv",row.names = FALSE)
head(table2_JAGS)
## Year
## colnames "Year" "Median_Biomass" "Q10_Biomass" "Q90_Biomass"
## B[1] "1950" "3190117.00660066" "1860414.28708401" "4708175.46375863"
## B[2] "1951" "3168594.83893451" "1849383.73826674" "4782152.68697437"
## B[3] "1952" "3161763.85701685" "1805774.18299977" "4820149.74975922"
## B[4] "1953" "3141856.76505951" "1788510.65811151" "4864082.55010801"
## B[5] "1954" "3149113.40309497" "1789102.92435557" "4929100.62993945"
##
## colnames "Median_CPUE" "Q10_CPUE" "Q90_CPUE" "Median_Depletion"
## B[1] NA NA NA "1.00084200375943"
## B[2] NA NA NA "0.999628179735674"
## B[3] NA NA NA "0.996590477601626"
## B[4] NA NA NA "0.995705455867898"
## B[5] NA NA NA "0.998697141341826"
##
## colnames "Q10_Depletion" "Q90_Depletion"
## B[1] "0.865686705046667" "1.1595816531195"
## B[2] "0.828100104355455" "1.20348122211214"
## B[3] "0.804868643319409" "1.23562497286099"
## B[4] "0.787418210087664" "1.25214038673299"
## B[5] "0.781025605961017" "1.2679635276962"
Median_Biomass <- table2(B_matrix)[,1]
Q10_Biomass<- table2(B_matrix)[,2]
Q90_Biomass<- table2(B_matrix)[,3]
Median_K <-table1(K_sample)[1]
plot(Year, Median_Biomass, type ="n", ylim=c(0, 5000000),
ylab = "Biomass", xlab ="Year")
polygon(c(Year, rev(Year)),
c(Q10_Biomass, rev(Q90_Biomass)),
col = "gray90", border = "gray90")
par(new=T)
plot(Year, Median_Biomass, type ="l", ylim=c(0, 5000000),
ylab = "", xlab ="")
abline(h=median(Median_K), col="pink")
Bratio <- Median_Biomass/table1(BMSY_sample)[1]
Fratio <- (Catch*sc)/table1(MSY_sample)[1]
YearLab <- paste(c(seq(50,99,length=50),00,01,02,03,04,
05,06,07,08,09,seq(10,19,length=10)))
plot(Bratio, Fratio, type = "n", xlim = c(0,3), ylim = c(0,2),
xlab = "B/Bmsy", ylab="U/Umsy")
rect(0, 0, 1, 1, col="yellow")
rect(1, 0, 3, 1, col="green")
rect(0, 1, 1, 2, col="red")
rect(1, 1, 3, 2, col="orange")
par(new=T)
plot(Bratio, Fratio, type = "l", xlim = c(0,3), ylim = c(0,2),
xlab = "", ylab="")
for(i in 1:length(Bratio)){
points(Bratio[i], Fratio[i], cex = 3, pch = 19, col="gray95")
text(Bratio[i], Fratio[i], cex = 1.3, YearLab[i])
}