From 0911ff6ef5ef6aebf3edd8475afb40287bd6fb9d Mon Sep 17 00:00:00 2001 From: Eric Fletcher <37851243+efletcherPIFSC@users.noreply.github.com> Date: Fri, 30 Jan 2026 15:48:08 -1000 Subject: [PATCH] Removed jabbaz.R Need to assess how the functionality of jabbaz.R and how central the script is to jabba before removing it, An alternative method is to have "relative" paths, rework source data. Related to jabbamodel/JABBA#13. --- jabbaz.R | 223 ------------------------------------------------------- 1 file changed, 223 deletions(-) delete mode 100644 jabbaz.R diff --git a/jabbaz.R b/jabbaz.R deleted file mode 100644 index f8bba7c..0000000 --- a/jabbaz.R +++ /dev/null @@ -1,223 +0,0 @@ - - -setwd("C:/Work/Research/Laurie/JabbaZ") -load("C:/Work/Research/Laurie/JabbaZ/om.75.RData",verbose=T) -library(FLCore) -library(FLBRP) - -library(JABBA) -plot(om) - -stk = iter(om,1) -plot(stk) -stk = window(stk,start=70,end=115) -# priors -fmsy = an(refpts(eq)["msy","harvest"]) -bmsy = an(refpts(eq)["msy","ssb"]) -b0 = an(refpts(eq)["virgin","ssb"]) -msy = an(refpts(eq)["msy","yield"]) -shape = bmsy/b0 - - -# Always convert Z to rate z = 1-exp(-Z) -zobs = function(flq,frange = "missing"){ - out=NULL - CN = flq - if(missing(frange)){ - Amin = which(apply(CN,1,median)==max(apply(CN,1,median))) - Amax = nrow(CN)-1 - frange= Amin:Amax - } - - for(t in 1:ncol(CN)){ - Ct= CN[,t] - df = as.data.frame(log(Ct[frange,])) - Z = -coef(lm(data~age,df))[[2]] - df = df[1,] - df$data = 1-exp(-Z) - df$age = "all" - out = rbind(out,df) - } - return(out) -} - -years = an(dimnames(stk)$year) - -# Prepare df -# Important don't use too many age-class because of lag -# If I choose ages 5:15 the default lag 1 works - for 5:20 I need lag2+ -set.seed(1234) -df = as.data.frame(FLQuants(Catch=catch(stk), - Index= vb(stk)*ar1rlnorm(0,years,sdlog=0.25), - Z1=zobs(catch.n(stk),frange=5:15)*ar1rlnorm(0,years,sdlog=0.1), - Z2=zobs(catch.n(stk),frange=5:15)*ar1rlnorm(0,years,sdlog=0.2) - )) - -om = NULL -om$stock = as.data.frame(ssb(stkfw)/bmsy) -om$harvest = as.data.frame(fbar(stkfw)/fmsy) -om$catch = as.data.frame(catch(stkfw)) - -om$fmsy = fmsy -om$bmsy = bmsy -om$msy = c(refpts(eq)["msy","yield"]) - -catch = Catch -index = Index -catch.n = as.data.frame(catch.n(stk)) -catch.n = catch.n[catch.n$year>=85,] - -save(catch,index,catch.n,om,file="data/pol.sim.rdata") - - -ggplot(df[df$qname%in%c("Z1","Z2"),],aes(year,data,color=qname))+geom_line() - - -inp = reshape2::dcast(df,year~qname,value.var="data",fun.aggregate=sum) -Index = inp[,c("year","Index")] -Index[Index==0] = NA -Index$Index[1:15] = NA # exclude observations -Index$Index = Index$Index/mean(Index$Index,na.rm=T) -Catch = inp[,c("year","Catch")] -Z= inp[,c("year","Z1","Z2")] -Z.se = Z -Z.se[,2:3] = 0.1 - -Z[Z==0] = NA -Z[1:15,-1] = NA # exclude observations - -# Fit with Catch + Index: Simple Fox with r = Fmsy -jbinput= build_jabba(catch=Catch,cpue=Index,auxiliary=Z, - auxiliary.se = Z.se, - model.type = "Fox",scenario="Index", - r.prior = c(fmsy,0.2), - auxiliary.type = "z", - auxiliary.lag = 4, - shape.CV = 0.2, - auxiliary.obsE = 0.01,auxiliary.sigma = T,verbose=F, - psi.prior=c(0.7,0.3)) - -fI= fit_jabba(jbinput,quickmcmc = T,verbose=F) -jabba=fI -jbplot_cpuefits(fI,index=1:2) -jbplot_runstest(fI,index=1:2) -jbplot_residuals(fI) - -prj1 = fw_jabba(fI,quant="F",type="ratio",imp.yr =2,stochastic = T,imp.values=1,nyears = 15,AR1=T ) - -jbplot_ensemble(prj1) - -# Compare: Black dashed line is TRUE pop dyn from OM -stkfw = iter(om,1) -jbpar(mfrow=c(2,2)) -jbplot_ensemble(prj1,subplots = 1,add=T) -lines(data~year,data=as.data.frame(ssb(stkfw)/bmsy),lwd=2,lty=2) -jbplot_ensemble(prj1,subplots = 2,add=T) -lines(data~year,data=as.data.frame(fbar(stkfw)/fmsy),lwd=2,lty=2) -jbplot_ensemble(prj1,subplots = 5,add=T) -jbplot_ensemble(prj1,subplots = 6,add=T) -lines(data~year,data=as.data.frame(catch(stkfw)),lwd=2,lty=2) - - - -# Fit Catch + Index + Z, testing with lag = 4 -jbIZ.lag4 = build_jabba(catch=Catch,cpue=Index,auxiliary = Z, - model.type = "Fox",scenario="Ind+Z.lag4", - r.prior = c(fmsy,0.3), - auxiliary.type = "z", - psi.prior=c(0.7,0.3), - auxiliary.obsE = 0.1, - auxiliary.sigma=TRUE, - auxiliary.lag = 4, # Lag4 - qA.cv=0.1, - verbose=F) - - -# Laurie - how to objectively check lag of F to Z response in stock? -error = ar1rlnorm(0,years,sdlog=0.) -Z35 =(zobs(stk,frange=5:35)*error) -Z20 =(zobs(stk,frange=5:20)*error) -Z15 =(zobs(stk,frange=5:15)*error) -Z10 =(zobs(stk,frange=5:10)*error) -Zom = apply(1-exp(-z(stk)[2:29]),2,median) -invI= 1/(vb(stk)*error) -invI = invI/mean(invI)*mean(Z35) -plot(invI,Z35,Z10,Zom) - - -# Fit Catch + Index + Z, testing with lag = 4 -jbIZ.lag4 = build_jabba(catch=Catch,cpue=Index,auxiliary = Z, - model.type = "Fox",scenario="Ind+Z.lag4", - r.prior = c(fmsy,0.3), - auxiliary.type = "z", - psi.prior=c(0.7,0.3), - auxiliary.obsE = 0.2, - auxiliary.lag = 4, # Lag4 - verbose=F) - -fIZ= fit_jabba(jbIZ.lag4,quickmcmc = T,verbose=F) - -# Fit only relative H ~ qZ, H is effort proxy with lag 4 -jbIE.lag4 = build_jabba(catch=Catch,cpue=Index,auxiliary = Z, - model.type = "Fox",scenario="Ind+Eff.lag4", - r.prior = c(fmsy,0.3), - auxiliary.type = "effort", - auxiliary.obsE = 0.2, - auxiliary.lag = 4, - verbose=F, - psi.prior=c(0.7,0.3)) - -fIE= fit_jabba(jbIE.lag4,quickmcmc = T,verbose=F) - - -# Compare: Black dashed line is TRUE pop dyn from OM -jbpar(mfrow=c(2,2)) -jbplot_ensemble(list(fI,fIZ,fIE),subplots = 1,add=T) -lines(data~year,data=as.data.frame(ssb(stk)/bmsy),lwd=2,lty=2) -jbplot_ensemble(list(fI,fIZ,fIE),subplots = 2,add=T) -lines(data~year,data=as.data.frame(fbar(stk)/fmsy),lwd=2,lty=2) -jbplot_ensemble(list(fI,fIZ,fIE),subplots = 5,add=T) -jbplot_ensemble(list(fI,fIZ,fIE),subplots = 6,add=T) - -# Check fits -jbpar(mfrow=c(2,2)) -jbplot_cpuefits(fIZ,add=T) -jbplot_cpuefits(fIE,add=T) - - - - -# Compare Lags -lags = 1:4 -jbs = lapply(lags,function(x){ - jb = build_jabba(catch=Catch,cpue=Index,auxiliary = Z, - model.type = "Fox", - scenario=paste0("Ind+Z.lag",x), - r.prior = c(fmsy,0.3), - auxiliary.type = "z", - psi.prior=c(0.7,0.3), - auxiliary.obsE = 0.2, - auxiliary.lag = x, # Lag2 - verbose=F) - fit_jabba(jb,quickmcmc = T,verbose=F) -}) - -jbpar(mfrow=c(2,2)) -jbplot_ensemble(jbs,subplots = 1,add=T) -lines(data~year,data=as.data.frame(ssb(stk)/bmsy),lwd=2,lty=2) -jbplot_ensemble(jbs,subplots = 2,add=T) -lines(data~year,data=as.data.frame(fbar(stk)/fmsy),lwd=2,lty=2) -jbplot_ensemble(jbs,subplots = 5,add=T) -jbplot_ensemble(jbs,subplots = 6,add=T) - -jbpar(mfrow=c(2,2)) -for(i in 1:4){ - jbplot_runstest(jbs[[i]],index=2,add=T) - legend("topleft",paste("Lag",i),bty="n") -} - -jbpar(mfrow=c(2,2)) -for(i in 1:4){ - jbplot_cpuefits(jbs[[i]],index=2,add=T) - legend("topleft",paste("Lag",i),bty="n") -}