computeFbrps gives negative yield and recruitment

[hgerritsen]

I was running a bootstrap SR on my stock object and some of the runs resulted in reference points that were higher than fcrash and which gave negative yield and recruitment. Presumably the correct behaviour in this case should be zero yield and recruitment.
Here is some code to illustrate the issue:

library(stockassessment)
library(FLCore)
library(FLfse)
library(FLRef)

SAMfit <- fitfromweb('whg.7b-ce-k_WGCSE22_RevRec_2023')
stk <- SAM2FLStock(SAMfit,catch_estimate = T)
`
# one of the iterations in my bootstrap that gives strange results
i <- c(19L, 10L, 6L, 24L, 14L, 2L, 13L, 18L, 22L, 14L, 6L, 1L, 19L, 
            19L, 8L, 6L, 23L, 12L, 6L, 8L, 7L, 11L, 17L, 4L)  
y <- as.FLSR(stk,model='ricker')
rec(y) <- rec(y)[, i]
ssb(y) <- ssb(y)[, i]
sr <- srrTMB(y, spr0 = mean(spr0y(stk)))

brps <- computeFbrps(stock = stk, sr = sr, proxy = 'sprx', f0.1 = TRUE)
brps@refpts
# fmax and spr.30 are higher than fcrash and give negative yield and recr

I think the steepness of the SR is so low that the entire SR curve is below the replacement line:

# following https://www.fao.org/3/v8400e/V8400E02.htm - plot replacement line (fig 6)
fi <- an(fbar(brps))
## SPR vector corresponding to these F's
spri <- an(spr(brps))
## reference points
rpi <- refpts(brps)

par(mfrow = c(1, 2))
## SPR
plot(fi, spri, type = "l", bty = "l", xlab = "Fbar", ylab = "SPR", ylim = c(0, max(spri)))

## SR
pars <- params(brps)
curve(pars["a"] * x * (exp(-pars["b"] * x)),
      from = 0, to = 1e5, bty = "l",
      xlab = "SSB", ylab = "Recruitment", xaxs = "i", yaxs = "i",
      main = "No point of intersection with replacement line")
abline(c(0, 1/spri[idx]), col = "blue")

[Henning-Winker]

Could you send me the stk and the sr as .rdata object please?

save(stk,sr,file="ricker.negative.rdata")
brps <- computeFbrps(stock = stk, sr = sr, proxy = 'sprx', f0.1 = TRUE)

For now, I can just speculate that it is not a bug but a property of the ricker given your data.Because the negative values occur on Fmax and spr30, which are quantities computed by FLBRP, I suspect this has to do with FLBRP 

 brps <- computeFbrps(stock = stk, sr = sr, proxy = 'sprx', f0.1 = TRUE)

Just a comment on the above: if x from the sprx is not specified it is x= 40 by default, i.e. spr40. I do not recall that f0.1 is a command. If you want to compare that in ploteq rather usebrps <- computeFbrps(stock = stk, sr = sr, proxy = c('sprx'.'f0.1'),x=40ploteq(brps)

Note that similar properties occur also for segmented regressions (hockey-stick)This is If the estimated break point b is larger than the SSB corresponding to biomass at at spr30, i.e. b > R0*spr30 (or the biomass at Fmax),This is because the hockey-stick has the property of week depensation for SSB < b, thus you would crash the stock for F that leads to an equilibrium biomass below b = Blim.
I am not sure how Ricker looks likes but it is often spurious and implausible I  would generally advice against unless there is strong biological evidence backed up by independent data.

[hgerritsen]

Hi Henning, thanks for looking into this.
ricker.negative.zip
stock and sr attached (as zip because gitgub doesnt allow .rdata files.)
In this bootstrap iteration it happened to be the ricker but i think in other iterations I also got negative yield for BH.