Add spawner-recruitment with "a, b" formulation; Also major refactoring and [Bug]: new features for impact of timevary biology on benchmark calcs

[Rick-Methot-NOAA]

paper by Miller and Brooks advocates for using the "a", "b" formulation of SRR curves, rather than steepness, when there is time varying biology which affects the spawners per recruit calculation.

A related / affected issues is #341 and #625

The fundamental issue is that if biology changes, the SSB/R changes. So, the equilibrium point with no fishing changes even though the SRR is unchanged. This is true for the steepness and the a,b formulation. However, what gets reported in mgmt_quant(1) as SSB_unf does NOT take this into account; it simply takes the benchmark SSB/R and multiplies it by the R0. The correct reporting of the equilibrium benchmark R0,B0 only appears in the SPR/YPR profile.

With changing adult biology, the SS3 paradigm asserts that R0 is unchanged, even though the equilibrium virgin biomass that results from that R0 is changed by the changed biology. With SS3 strongly organized around use of R0 and initial equilibrium and spawner-recruitment, the MSY calcs from SS3 are affected. Additional controls are needed in order to elucidate the relationships and some of the benchmark quantities were found to be in error when there is changing biology. The changes made in this issue provide more control over whether the initial year (aka virgin) SSB/R is used or the contemporary (recent) biology is used in various relevant places. The options provide backward compatibility to the historical approach.

With the a,b approach, R0 is a derived quantity. With the steepness approach, it is a leading parameter. If a,b is used, then the R0, h are calculated using start year biology and copied to their normal location as MGparms 1 and 2.

If R0 or h themselves change, then there is a question of whether the SRR should continue to use start year bio, or bio at that year. The intention before 3.30.34 was to use the updated biology, but unclear whether it was fully implemented. If the a,b parameterization was used then changing biology is not invoked as a or b changes, although the a,b values are influenced by the biology. But with R0,h there is need to specify SSBpR0. Three possibilities when R0 or h changes:

1. keep SSBpR0 from virgin even if biology has changed
2. use SSBpR0 from the year of the R0 or h change. SS3 now tracks this year with a "1" in timevary_SRparm(y)
3. update SSBpR0 each year beginning in the year of the R0, h change.
   I think option 2 is logical and is consistent with how the benchmark calculations would use benchmark biology along with any changes R0, h. Options 1 and 3 will not be implemented at this time.

Task List:

• add SRR function 10 to add a alpha, beta parameterization

• substantially refactor the variables names to clarify SPR and SSB/R usages

• add input in SRR section to allow for specification regarding biology used in SSBpR0, it was always using the benchmark biology for SSBpR0, but should use virgin biology to maintain the same SRR as is used in the time series.

• fix to biology used in SPR profile and global MSY; it was using biology at end of forecast, not the requested benchmark biology; fix by referring to values stored during benchmark at styr-3

• address biology used in calculation of depletion as it affects the control rule inflection by adding HCR_anchor control in forecast.ss

• output the benchmark and forecast controls to report.sso

• add depletion option 6 that will use benchmark SSB_unf as basis

• clarify that the internal depletion values are reported as Bratio

• clarify and report which SSB_unf is used for spawner-recruitment calcs, which for depletion, and which for HCR_anchor

• output SSB for HCR_inflection to mgmt_quant,

• Change approach for invoking Bmsy as the HCR_inflection. Instead of using -1 as HCR_inflection value, now invoke as option 3 for HCR_anchor. Retain backward compatibility, but write ss_new in new approach only.

• add 5 elements to mgmt_quant:
  Recr_MSY_bmarkbio
  Recr_unfished_bmarkbio
  SSB_unfished_bmarkbio // note that this is same as element 1: SSB_unfished
  SSB_HCR_infl
  SSB_virgin_again // note this is redundant with SSB_virgin in the #_Spawn_Bio section. Redundant items may be removed after testing of feature is complete. @iantaylor-NOAA

• address which wt_age_beg is used for benchmark smrybio when there is time-vary biology

• change mgmt_quant(1) to report the equilibrium SSB, not simply the R * SSB/R.

• create options for setting Btgttgt and Btgttgt2 in benchmark so can use SSB_bench or SSB_virgin; current default is SSB_bench. Solution is to expand depletion basis options to include new (6) to use SSB_bench as basis. Depletion basis already included option (1) SSB_virgin. Now the choice made for depletion basis will also be used for Btarget and Btarget2 basis.

• rename Do_Equil_Calc to SSBpR_Calc.

• remove extra application of regime in the initial equilibrium (R1) calculations when SRR is used.

• SS3 lacks a entity that records the first year of timevary SRR parameters, which is the first year for which changed biology could be incorporated into the updated SRR. Create internal quantity: timevary_SRparm(y) and timevary_bio_4SRR. The later is now read from starter.ss with 0 for legacy approach and 1 to invoke the new revised code. The following warnings are produced if legacy is used where it should not be:
  Warning 10 : There is timevary biology and the legacy approach to benchmark calculations is being used; user should be aware of possible impacts to benchmark results
  Warning 11 : There is timevary biology, so the flag for timevary_bio_4SRR must be set to 0 (old default) or 1 (new improved) in starter.ss

• Forecast uses SRparm_work as continued from timeseries; it does not use the averaged values used in benchmark. Make this more explicit by creating SRparm_bench as a specific version of SRpark_work to be used for benchmark calcs, and by passing a version of SRparm into Spawn_Recr rather than using the global variable.

• The spawn_recr curve reported in report.sso uses styr biology and virgin SSBpR. It seems useful to also show the SR curve using benchmark quantities and perhaps the SRcurve as it exists in last year of the forecast.

• correct the label in starter.ss_new for SPR_reporting. Should be: SPR_report_basis: 0=skip; 1=(1-SPR)/(1-SPR_tgt); 2=(1-SPR)/(1-SPR_MSY); 3=(1-SPR)/(1-SPR_Btarget); 4=1-SPR. Also add option 5=SPR to get the raw %SPR values in timeseries to match the units reported in benchmark.

• Add SPR_reporting option 5 to the manual. Other entries in manual are correct.

Task List for after coding this is completed:

• check consistency of results with info in spawn_recr curve and spr/ypr profile
• do testing with EWAA model (hake)
• Get SS3 power users to test this out and provide feedback
• Update user manual
• Create a short one-pager that people can bring to the council that indicates that there was not a bug but incomplete implementation in SS3

[k-doering-NOAA]

Miller and Brooks 2021

[Rick-Methot-NOAA]

Brainstorming:
set alpha = parm(3) and beta = parm(4) as leading, estimated parameters
then R0 = parm(1) and steepness = parm (2) can be derived and reported per existing outputs

with unexploited spawning biomass per recruit = ϕ0
then

and

But this could only be used to report h and R0 conditioned on there not being time-varying biology. If biology is time-varying, for a given a,b there will be a different implied steepness and R0 for each year. This could easily be reported in the SPR_SERIES table because that table already calculates SSBzero from annual biology. Note that annual biology is the result of the time series of annual biology parameters, it is not the equilibrium result of annual biology parameters. That should be an option.

[iantaylor-NOAA]

Sounds good to me.
Would it be worth considering the reverse case where the alpha and beta are available as reported as derived quantities when you use R0 and h as estimated parameters? That could presumably be added in the future if someone wants it.

[Rick-Methot-NOAA]

got it working and produces same result as R0, h approach.
implemented parameters as ln(alpha) and ln(beta)
code calcs derived R0 and steepness and saves them to SRparm(1) and SRparm(2)
by copying results to SRparm(2), existing priors on steepness should still be usable.
Update: I ran a test with a strong prior on derived steepness and it was successful in causing a changed estimates of the a, b parameters.

[iantaylor-NOAA]

This is excellent, thank you.

[Rick-Methot-NOAA]

attached is spreadsheet in which I investigate plausible ranges for the alpha and beta parameters.
@timjmiller/wham Do you have any tips on parameterizing the alpha beta?
test_AB_range.xlsx