SRR

New srrTMB algorithm

library(FLCore)
library(FLBRP)
library(ggplotFL)

require(FLasher)

require(FLCandy)

library(FLSRTMB)

library(plyr)
library(dplyr) 
library(reshape) 

library(akima)

library(scales)

library(ggplot2); theme_set(theme_bw())

setwd("~/Desktop/active/srr/Rmd")

Modifications

1 Call function with parameters/priors for steepness, \(R_0\), sigma R and sigma \(R_{AR1}\)

2 Allow \(SPR_0\) to vary by year

3 Fit using h, \(B_0 = SPR_0*R_0\)

4 Return time varying \(\alpha\) & \(\beta\)

5 Allow use of constricted likelihood, by estimating sigma R and sigma RAR1 iteratively

6 Allow \(R_0\) to be estimated in blocks

Then

• Fit with and without time varying parameters and compare SRRs using AIC
• Use time varying SRRs to estimate time varying reference points
• Use MASE to compare ability to predict reference points
• Run hindcast of ICES HCR

Single Stock

Data

load("/home/laurie/Desktop/active/biomassrebuild/data/Updated_stks_n81_R0_updated2023_V2.RData")

stksel = c("pil.27.8c9a",
           "pra.27.3a4a",
           "her.27.3a47d",
           "ple.27.420",
           "cod.27.6a",
           "mon.27.78abd",
           "mac.27.nea",
           "whb.27.1-91214",
           "bss.27.8ab",
           "sol.27.4",
           "cod.27.5a",
           "cod.27.22-24",
           "hke.27.3a46-8abd",
           "sol.27.7a")[-c(2,5,9,10,11)]

stksel = unique(c(stksel,c("ank.27.78abd","her.27.3031","hke.27.3a46-8abd","hke.27.8c9a","her.27.25-2932")[-3]))

stks=ICESStocks[stksel]
bm=ldply(stks, function(x) as.data.frame(x@benchmark))
es=ldply(stks, function(x) as.numeric(x@eqsim))
es=transform(es,Catchequi     =an(V2),
                BMSY          =an(V3),
                B0            =an(V4),            
                FmsyMedianC   =an(V5),
                FmsyMedianL   =an(V6),
                F5percRiskBlim=an(V7), 
                FlimEqsim     =an(V8), 
                R0            =an(V9))[,-c(2:10)]
rf=merge(bm,es)
fl=ldply(stks, function(x) attributes(x)$fishlife)
rf=merge(rf,fl)

spr0=ldply(stks,function(x) data.frame(spr0=c(mean(spr0Yr(x)))))
srrs=llply(stks, function(x) as.FLSR(x,model="bevholt"))

Fits

x    =stks[["sol.27.7a"]]
spr0 =mean(spr0y(x))
r0   =c(an(x@eqsim["R0"]),0.2)
srs  =FLSRs(
  s.est=srrTMB(as.FLSR(x,model=bevholtSV),spr0=spr0,r0.pr=r0),
  s0.6 =srrTMB(as.FLSR(x,model=bevholtSV),spr0=spr0,r0.pr=r0,s.est=FALSE,s=0.6),
  s0.9 =srrTMB(as.FLSR(x,model=bevholtSV),spr0=spr0,r0.pr=r0,s.est=FALSE,s=0.9),
  s.d  =srrTMB(as.FLSR(x,model=bevholtDa),spr0=spr0,r0.pr=r0,                   ld=1,d=1.5,d.logitsd=1),
  s.d2 =srrTMB(as.FLSR(x,model=bevholtDa),spr0=spr0,r0.pr=r0,                   ld=1,d=2.0,d.est=FALSE))

plot(srs)

Diagnostics

AIC

ldply(srs,AIC)

    .id       V1
1 s.est 106.4978
2  s0.6 107.0030
3  s0.9 119.7675
4   s.d 101.1197
5  s.d2 101.1081

\(SPR_0\)

rodFn=FLife:::rodFn

spr0=FLQuants(llply(stks["sol.27.7a"], function(x) FLCandy:::spr0Yr(x)))
rgm =ldply(spr0, FLife:::rod)

ggplot(transform(as.data.frame(spr0,drop=TRUE),.id=qname))+
      geom_point(aes(year,data),col="grey25")+
      geom_line( aes(year,data),col="grey75")+
      geom_polygon(aes(year,data,group=regime),
          fill="lavender",col="blue",
          lwd=.25,alpha=.5,
          data=rgm)+
      facet_wrap(.id~.,scale="free",ncol=2)+
      xlab("Year")+ylab("SPR0")+
      theme_bw(16)+
      theme(axis.text.x=element_text(angle=45),legend.position="bottom")

Figure 1 SPR0

Residuals

rsd=FLQuants(llply(srs, function(x) residuals(x)))
rgm=ldply(rsd, FLife:::rod)

ggplot(transform(as.data.frame(rsd,drop=TRUE),.id=qname))+
      geom_point(aes(year,data),col="grey25")+
      geom_line( aes(year,data),col="grey75")+
      geom_polygon(aes(year,data,group=regime),
          fill="lavender",col="blue",
          lwd=.25,alpha=.5,
          data=rgm)+
      facet_wrap(.id~.,scale="free_y",ncol=3)+
      scale_y_continuous(labels=percent)+
      xlab("Year")+ylab("SPR0")+
      theme_bw(16)+
      theme(axis.text.x=element_text(angle=45),legend.position="bottom")

Figure 2 Recruitment residuals.

ACF

par(mfrow=c(3,2))
l_ply(srs, function(x) acf(residuals(x)))

Figure 3 ACF

MASE

Reference points

rfs=ldply("sol.27.7a", function(.id) 
      cbind(.id=.id,nonStationarity(stks[[.id]],srs[[1]])))
rf2=ddply(rfs, .(.id,refpt,quant), transform, data=data/mean(data,na.rm=T))

ggplot(subset(rf2,quant=="ssb"&refpt%in%"msy"))+
    geom_point(aes(year,data),col="grey25")+
      geom_line( aes(year,data),col="grey75")+
      scale_y_continuous(labels=percent)+
      facet_wrap(.id~.,scale="free_x",ncol=3)+
      xlab("Year")+ylab("MSY")+
      theme_bw(16)+
      theme(axis.text.x=element_text(angle=45),
            legend.position="bottom")

Figure 4 MSY

ggplot(subset(rf2,quant=="ssb"&refpt%in%"virgin"))+
    geom_point(aes(year,data),col="grey25")+
      geom_line( aes(year,data),col="grey75")+
      facet_wrap(.id~.,scale="free_x",ncol=3)+
      scale_y_continuous(labels=percent)+
      xlab("Year")+ylab("Virgin")+
      theme_bw(16)+
      theme(axis.text.x=element_text(angle=45),
            legend.position="bottom")

Figure 5 Virgin

ggplot(cast(subset(rfs,quant=="ssb"&refpt%in%c("virgin","msy")),
          .id+year~refpt,value="data"))+
      geom_point(aes(year,msy/virgin),col="grey25")+
      geom_line( aes(year,msy/virgin),col="grey75")+
      facet_wrap(.id~.,scale="free_x",ncol=3)+
      xlab("Year")+ylab("BMSY/BVirgin")+
      theme_bw(16)+
      theme(axis.text.x=element_text(angle=45),
            legend.position="bottom")

Figure 6 Ratio of \(B_{MSY}/B_0\)

Multiple Stocks