spawner-recruitment with "a, b" formulation and 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

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
• 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
• do testing with EWAA model (hake)
• 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
• add recruitment at MSY to the mgmt_quant
• clarify and report which SSB is used for spawner-recruitment calcs, which for depletion, and which for HCR_anchor
• add SSB for HCR_anchor to mgmt_quant,

[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

[iantaylor-NOAA]

Trying again to tag @timjmiller and @liz-brooks for any tips on initial values and ranges for the alpha beta parameterization of the Beverton-Holt SRR. See messages and attachment from @Rick-Methot-NOAA in this github thread.

[Rick-Methot-NOAA]

Thanks Ian. In my first trials, I can get the same answer (haven't done any trials with time-varying biology yet) as with the R0,h approach, but convergence path is ragged. (haven't done any trials with time-varying biology yet). ADMB makes some bad jumps early in the search even though the initial values for parms were not far off. Definitely worse than performance with R0,h.
I haven't done any trials with time-varying biology yet to replicate the slippery slope findings.
What I think we need is a parameterization with the leading parameters being Rbar and parm2 where parm2 provides curvature. This will be tough because expected curvature depends on Rbar/R0, which we don't know.

[Rick-Methot-NOAA]

something got broken in MSY calcs

[timjmiller]

in wham initial value for ln(a) = 0. For B-H initial value for ln(b) = 0 and for Ricker ln(b) = -10

maybe a better guess would be by determining a and b from some initial guess for R0 (e.g. exp(10)) and h (e.g., 0.7) and unfished ssb/r (averaged over the time series if necessary).

You probably already have these but eqs for getting a and b are here for BH and here for Ricker

[Rick-Methot-NOAA]

Thanks Tim. I got off-track for awhile until realizing that SS3 already had some alpha, beta code based on R = aS/(b+S) such that your a,b and mine did not match. I have now aligned SS3 a,b to be same as yours.

[liz-brooks]

are you all set on this based on tim's reply? i let this slip too far down the queue, but wasn't sure if you "close" issues on this list serve or not. cheers liz

…

On Mon, Feb 13, 2023 at 9:09 PM Richard Methot ***@***.***> wrote: Thanks Tim. I got off-track for awhile until realizing that SS3 already had some alpha, beta code based on R = aS/(b+S) such that your a,b and mine did not match. I have now aligned SS3 a,b to be same as yours. — Reply to this email directly, view it on GitHub <#191 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACFL5LIZANJAW5BASJVP3JLWXLSMDANCNFSM5ADSFN2A> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Liz Brooks, PhD Operations Research Analyst Population Dynamics Branch NOAA/NMFS Northeast Fisheries Science Center 166 Water Street phone: 508.495.2238 Woods Hole, MA 02543 fax: 508.495.2393 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

[iantaylor-NOAA]

NWFSC and SWFSC folks discussed this new Brooks (2024) paper on recruitment https://doi.org/10.1016/j.fishres.2023.106896 which renewed interest in availability of the a-b (alpha-beta) parameterization of the Beverton-Holt stock-recruit relationship.

@Rick-Methot-NOAA, do you think this is feasible to get into the next release and is there any support you need to do so?

[Rick-Methot-NOAA]

We need careful consideration of implications for depletion calculations. They already are somewhat wrong when there is time-varying growth.

[Rick-Methot-NOAA]

Tim,
The SS3 implementation of the R0h approach does not update the unfished SSB/R when there is time-varying biology (except in the rare case where there is a regime shift in R0). So, it is equivalent in performance to the a,b which implicitly finds values for those correlated a,b parameters that find the same shape and elevation for the curve. I'll work this comparison into the demonstration I'm developing.

[timjmiller]

Rick, That makes sense that for a given year of biology you get equivalent results for either a,b or R0,steepness parameterization. You can translate from one to the other for expected recruitment. I'm pretty sure this is what ASAP does, too. The issue though is that one needs to keep track of which year of biology is used to define R0 (and S0) and steepness. For example, if one is using more recent biology (that is different) to calculate SSB/R and Y/R for management reference points, you would not want to use the estimated steepness and R0 in those calculations because they are functions of different biology parameters. You would need to translate the earlier R0 and steepness into a,b and then calculate the new R0 and steepness with the current biology to be comparable to other equilibria based on recent biology. Using the a,b parameterization makes the bookkeeping much easier, because the post-recruit biology only comes into the equation when calculating reference points or any other equilibrium points and it is easier to make sure the time-specific biology is consistent between different equilibrium attributes. e.g., unfished SSB and R and those at Fmsy. It also makes it easier to understand where the variation in equilibrium points is coming from (pre-recruit biology or post-recruit biology). A regime shift in R0 could be due to changes in post-recruit biology and/or the a,b parameters. If the shift is due to the former, we can calculate this change in R0 over time from annual calculations of R0 (and other equilibria) with the changing biology but constant a,b. I think one would want to be careful about whether the regime shift is in a and/or b parameters vs. post-recruit biology. Tim

…

On Tue, May 21, 2024 at 5:40 PM Richard Methot ***@***.***> wrote: Tim, The SS3 implementation of the R0h approach does not update the unfished SSB/R when there is time-varying biology (except in the rare case where there is a regime shift in R0). So, it is equivalent in performance to the a,b which implicitly finds values for those correlated a,b parameters that find the same shape and elevation for the curve. I'll work this comparison into the demonstration I'm developing. — Reply to this email directly, view it on GitHub <#191 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AEIGN7GEAFPBCI4LFOQ2ZXTZDO5N3AVCNFSM5ADSFN2KU5DIOJSWCZC7NNSXTN2JONZXKZKDN5WW2ZLOOQ5TEMJSGM2DQNJXHAZQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- Timothy J. Miller, PhD (he, him, his) Research Fishery Biologist NOAA, Northeast Fisheries Science Center Woods Hole, MA 508-495-2365

[liz-brooks]

i think there is confusion here between (1) what is assumed mathematically in deriving R0h and (2) how the SRR is coded in SS. for (1), you must define a replacement line to translate from ab to R0h. if none of the biological parameters vary over time, then there is a unique replacement line with slope = 1/spr0. if any of those biological parameters vary over time (M, maturity, "fecundity"), then there is not a unique replacement line. in that case, the choice of biological parameters impacts the calculation of spr0, and this in turn will influence where the replacement line intersects the SRR defined by ab and therefore produces different R0h estimates depending on which biological parameters you pointed to for the replacement line calculation. since the ab parameterization does not depend on the replacement line, fitting that parameterization will return the same estimates of ab even if biological parameters are time varying. for (2), in your previous message you said "The SS3 implementation of the R0h approach does not update the unfished SSB/R when there is time-varying biology..." it sounds like you're saying that if biological parameters are time varying and you ignore that, then there is no difference in the R0h estimates. but maybe i'm misinterpreting you? i think the point we have been trying to make is that it DOES make a difference. pointing to year 1 for the biological parameters is arbitrary, and if you pointed to a different year, you'd get a different estimate of R0h. your time series estimates of recruits and spawners would probably look very much the same (mosly driven by the data). but if you are bothering to estimate the SRR, then presumably you want to use them to derive reference points and make projections -- both of those would be impacted by having time-dependent R0h parameters. the replacement line assumption seems to be a behind the scenes assumption in the code that perhaps some users don't recognize (and maybe can't specify or change or explore). thus, if you move to the ab parameterization, there is no assumption that needs to be made about which biological parameters to use in fitting the SRR. then reference point specification and projections are separate steps where the analyst is cognizant (and transparent) about which biological parameters are assumed to derive the reference points (again, they define the calculation of spr0). cheers liz

…

On Tue, May 21, 2024 at 5:40 PM Richard Methot ***@***.***> wrote: Tim, The SS3 implementation of the R0h approach does not update the unfished SSB/R when there is time-varying biology (except in the rare case where there is a regime shift in R0). So, it is equivalent in performance to the a,b which implicitly finds values for those correlated a,b parameters that find the same shape and elevation for the curve. I'll work this comparison into the demonstration I'm developing. — Reply to this email directly, view it on GitHub <#191 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACFL5LLNHCDBW73MFB3EW5DZDO5N3AVCNFSM5ADSFN2KU5DIOJSWCZC7NNSXTN2JONZXKZKDN5WW2ZLOOQ5TEMJSGM2DQNJXHAZQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Liz Brooks, PhD Operations Research Analyst Population Dynamics Branch NOAA/NMFS Northeast Fisheries Science Center 166 Water Street phone: 508.495.2238 Woods Hole, MA 02543 fax: 508.495.2393 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

[Rick-Methot-NOAA]

Let me work out the examples before responding further. R0 and the R0,h approach are very deep in the DNA of SS3. Major changes to SS3 are not anticipated; instead these issues should be aimed at how best to develop the FIMS approach to use of spawner-recruitment.
Thanks Tim for noting that SS3's use of start year biology for calculating the SPR0 produces an equivalent curve as if a,b was used.
The issue is the MSY calculation in SS3 which allows user to specify a range of years for biology, but it does not update R0,h to be consistent with that biology. That inconsistency probably will be dealt with by calculating an a,b that is equivalent to the original R0,h; then using that a,b in the reference point calculation.
Both R0h and a,b produce recruits from SSB(y) where the SSB(y) is based on that year's biology. If the biology is quasi-randomly time-varying it's effect on the production of recruits from the spawner-recruitment curve should be the same whether that curve is parameterized as a,b or R0h. It's the same curve. Regime shifts are a different matter and do need a recalibration of the curve.

[msupernaw]

It may be useful looking at this formulation using the FIMS_statistical_computing_investigations tool.

here is an example analyzing the double logistic function.

here is the generated report.

[msupernaw]

I took the liberty of running the Beverton-Holt variation in the functional analysis tool. Assuming that I coded it correctly, it would appear that beta has a higher level of derivative stochasticity, which could make convergence for some quasi-newton minimizers more difficult. If you're interested, the code is here and the results are here. You may want to verify the fixed values for phi and S. Hope this helps with your investigation.

[Rick-Methot-NOAA]

Thanks for chiming in Matthew. I don't think I can use that for what I'm doing in SS3, but it seems very useful for eventual development in FIMS.

[Rick-Methot-NOAA]

I now have a comparison between SS3 using R0h and SS3 using AB. The comparison was run using a situation with time-varying growth such that fish body size increases substantially in later portion of the time series and into the projection.

The R0h approach gets expected recruitment from:

      NewRecruits = (4. * steepness * Recr_virgin_adj * SSB_curr_adj) /
          (SSB_virgin_adj * (1. - steepness) + (5. * steepness - 1.) * SSB_curr_adj);

where Recr_virgin_adj and SSB_virgin_adj are from the biology at start of the time series and are not (in this example) allowed to be updated with changing conditions.
The AB approach uses:
NewRecruits = (alpha * SSB_curr_adj) / (1.0 + beta * SSB_curr_adj);
So, the values of alpha and beta subsume to effect of virgin biology and steepness.

The results of the two runs were identical to rounding error. Same recruitment time series, same fit to data, same MSY, same projection.

The MSY calculations in both runs were set to use start year biology, a user control in SS3. This caused MSY to be much less than projected catch because the projected catch used the increased body weight, but the MSY did not. This difference needs to be made more apparent to the SS3 user that currently is done.

Switching the MSY calculation to use end year biology had an effect on Fmsy and MSY and BMSY and now is consistent with the level of MSY from the projection at FMSY. It still was the same between the Roh run and the AB run.

As I continue working on this, my emphasis will be on clearly altering users to what changes might be occurring due to time-varying biology with R0h. I'll keep the AB option, but not emphasize any need for users to switch away from R0h. Equivalent results can be obtained.

[liz-brooks]

hi rick, thanks for following up with these comparisons. your results are what we'd expect when you point to either the beginning or end of the time series for your biological parameters (though i don't think it would work if a prior is specified on h). is the user option to only point to a single year of biological parameters, or would it be possible to specify a range of years (e.g., to reflect average prevailing conditions)? a further thought, if you just fit average recruitment with deviations, are SPR calculations pointing to the first year by default or a user-specified year (or range of years) that is consistent with projections? out of curiosity, what were the estimates from those two cases for R0, h, Fmsy, Bmsy, MSY as well as terminal year estimates of F/Fmsy, B/Bmsy, and B/B0? i would expect that differences in the h estimate due to year selected for biological parameters will carry through to perceptions about sustainable stock levels and current stock status and depletion. alerting users to that specification sounds like a good plan. do you want to leave it up to the user to make individual runs exploring the consequence of the year they choose, or do you want to provide a table and/or plot in the results that automatically makes all of those translations of R0h so the user can see the range of values corresponding to biological parameters for each model year? i.e., since you can get back to ab from your R0h specification, and from ab you can then solve for the annual R0h (this is what is done and is reported in ASAP). cheers liz

…

On Fri, May 24, 2024 at 12:34 PM Richard Methot ***@***.***> wrote: I now have a comparison between SS3 using R0h and SS3 using AB. The comparison was run using a situation with time-varying growth such that fish body size increases substantially in later portion of the time series and into the projection. The R0h approach gets expected recruitment from: NewRecruits = (4. * steepness * Recr_virgin_adj * SSB_curr_adj) / (SSB_virgin_adj * (1. - steepness) + (5. * steepness - 1.) * SSB_curr_adj); where Recr_virgin_adj and SSB_virgin_adj are from the biology at start of the time series and are not (in this example) allowed to be updated with changing conditions. The AB approach uses: NewRecruits = (alpha * SSB_curr_adj) / (1.0 + beta * SSB_curr_adj); So, the values of alpha and beta subsume to effect of virgin biology and steepness. The results of the two runs were identical to rounding error. Same recruitment time series, same fit to data, same MSY, same projection. The MSY calculations in both runs were set to use start year biology, a user control in SS3. This caused MSY to be much less than projected catch because the projected catch used the increased body weight, but the MSY did not. This difference needs to be made more apparent to the SS3 user that currently is done. Switching the MSY calculation to use end year biology had an effect on Fmsy and MSY and BMSY and now is consistent with the level of MSY from the projection at FMSY. It still was the same between the Roh run and the AB run. As I continue working on this, my emphasis will be on clearly altering users to what changes might be occurring due to time-varying biology with R0h. I'll keep the AB option, but not emphasize any need for users to switch away from R0h. Equivalent results can be obtained. — Reply to this email directly, view it on GitHub <#191 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACFL5LLXXV74B4XLVB77HPTZD5TYVAVCNFSM5ADSFN2KU5DIOJSWCZC7NNSXTN2JONZXKZKDN5WW2ZLOOQ5TEMJSHE4TKOBWHA3Q> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Liz Brooks, PhD Operations Research Analyst Population Dynamics Branch NOAA/NMFS Northeast Fisheries Science Center 166 Water Street phone: 508.495.2238 Woods Hole, MA 02543 fax: 508.495.2393 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

[Rick-Methot-NOAA]

Yes, I have built into the reports the translation between R0h and AB.
R0 is extremely hardwired into the DNA of synthesis. So when it uses AB, it calculates R0 and h and stores them in the expected location so that R0 can be used to start off the population and serve as basis for some depletion options. the associated h is not used.
SS3 has long had specification of a range of years over which to average the biology when calculating reference points. Its shortcoming is that it is averaging the results of the biology (mean of weight-at-age, rather than weight-at-age from mean of changing M and growth parameters). That is being pursued in another issue. Changing to work from averaged parameters will also allow for inclusion of density-dependent growth and M which can happen now in the time series and projections but not in the reference point calculations.
Alerting users to the impact on depletion calculations is one of the biggest consequences. SS3 already has a fully developed dynamic B0 which probably will get recommended as the way to go whenever there is time-varying biology. Dynamic B0 follows changing biology as well as changing e(recruits); even density-dependent growth gets accounted for.
I do not anticipate calculating in SS3 an "h" value relevant to each year's biology. So when there is time-varying biology and the user asks for non-start year biology to be used in MSY there is no effect on the R0 and h because their values were estimated to be consistent with the time series of recruits produced by the assessment; same as happens with AB.
Liz - I resist putting biological interpretation on the parameters AB, although you and others have pursued that line of investigation. I see the spawner-recruitment function merely as a tool to calculate the degree of decline in R that is expected as a stock transitions from an unfished state to a fully fished state near Bmsy. R0h satisfies that need; better still with a 3rd parameter to allow uncertainty in the Bmsy/B0 ratio
So, while it is possible to translate the R0h parameters into AB, there is no need to do so when using R0h to describe the central tendency of R as SSB declines. It is just an alternative configuration of this 2 parameter functional form. Where we need to be careful is in stating explicitly the M and growth conditions under which the translation is occurring.

[Rick-Methot-NOAA]

@chantelwetzel-noaa @shcaba
I am in the process of greatly expanding the output associated with the SPAWN_RECR output table. In particular, I am emphasizing aspects affected by time-varying biology and time-varying spawner_recruitment parameters. I am still working on examples to show relationship between how it looks with the R0,h formulation vs the alpha, beta formulation. For now, I am just asking to see if this output makes sense. Much work needed in r4ss to read this revision.

`SPAWN_RECRUIT report:19
SR_Function: 3

parm parm_label value phase
1 SR_LN(R0) 8.81439 1 #_is_time_vary,_so_SRR_updates_base_SPR_annually
2 SR_BH_steep 0.632662 4
3 SR_sigmaR 0.6 -4
4 SR_regime 0 -4
5 SR_autocorr 0 -99

sigmaR: 0.6
SPR_virgin: 6.77161 #_uses_biology_from_start_year: 1971
Ln(R0): 8.81439
R0: 6730.37
steepness: 0.632662
Ln(alpha)_derived: 0.0172082 alpha 1.01736
Ln(beta)_derived: -8.95401 beta 0.000129218

Quantities for MSY and other benchmark calculations
Benchmark_years: 1_beg_bio 2_end_bio 3_beg_selex 4_end_selex 5_beg_relF 6_end_relF 7_beg_recr_dist 8_end_recr_dist 9_beg_SRparm 10_end_SRparm
Benchmark_years: 2001 2001 2001 2001 2001 2001 1971 2001 1971 2001
First_year_with_timevary_MG: 1990
#_range_of_years_is_averaged,_so_reduces_standard_error_of_result;_do_this_only_when_timevarying_makes_necessary: 7 8
#_range_of_years_is_averaged,_so_reduces_standard_error_of_result;_do_this_only_when_timevarying_makes_necessary: 9 10
SPR_unfished_benchmark: 13.5849 #_based_on_averaging_biology_over_benchmark_year_range
Bmsy/Bzero: 1.14525 # using styr bio for Bzero
Bmsy/Bunf: 0.570868 # using MSY's averaged bio for Bunf

RecDev_method: 2
sum_recdev: 0.708989
recr_logL: -8.61409
1971 2001 main_recdev:start_end
1900 1900 2001 2002 1 breakpoints_for_bias_adjustment_ramp
ERA N RMSE RMSE^2/sigmaR^2 mean_BiasAdj est_rho Durbin-Watson
main 31 0.409543 0.465903 1 -0.136275 2.25483
early 0 0 0 0

Initial_equilibrium: 0 # 0/1_to_use_spawner-recruitment_in_initial_equ_recruitment_calculation

Yr SpawnBio exp_recr with_regime bias_adjusted pred_recr dev biasadjuster era mature_bio mature_num raw_dev SPR0_curr h_curr R0_curr P1 P2 P3 P4
S/Rcurve 45575.5 6730.37
Virg 45575.5 6730.37 6730.37 6730.37 6730.37 - 0 Virg 45575.6 14930.4 0.0
Init 45575.5 6730.37 6730.37 6730.37 6730.37 - 0 Init 45575.6 14930.4 0
1971 45575.5 6730.37 6730.37 5621.68 5984.17 0.0624876 1 Main 45575.6 14930.4 0.0624876 6.77161 0.632662 6730.37 8.81439 0.632662 0.6 0
1972 45575.3 6730.37 6730.37 5621.68 3786.63 -0.395153 1 Main 45575.3 14928.5 -0.395153 6.77161 0.632662 6730.37 8.81439 0.632662 0.6 0
1973 45494.4 6728.63 6728.63 5620.23 7035.02 0.224529 1 Main 45494.4 14899.5 0.224529 6.77161 0.632662 6730.37 8.81439 0.632662 0.6 0`

[chantelwetzel-noaa]

I think this expanded output is relatively clear and includes all the needed information. Thanks!

[Rick-Methot-NOAA]

Below is the promised demonstration. My conclusion is that the R0h approach is different from the alpha-beta approach only when invoking time-varying SPR0, which is never done in assessments I have seen. Even then, the difference is small and not wrong, just different. Report below:
Gen data using AB approach and with time-vary growth (DATAGEN_TV). TVG has Linf increasing from ~72cm early in time series, to 80 cm in 1990 and beyond throughout the projection, which is 80 years to allow growth and recruitment to approach stable values following the change in growth

Estimate the model with AB approach, which reports the equivalent R0 and h (steep) and uses that R0 in starting the population.
Result:
logL 1367.46
SR_LN(R0) 8.86372 0 Fix 0.209419
SR_BH_steep_derived 0.492508 0 Fix 0.365634
SR_BH_ln(alpha) -0.650818 0 Act 0.39153
SR_BH_ln(beta) -9.81242 0 Act 0.168276

Estimation the model with R0 and h.
Result:
logL = 1367.46 and everything about the run is same as the AB run.
SR_LN(R0) 8.86372 0 Act 0.209419
SR_BH_steep 0.492508 0 Act 0.365634

In the Spawn_Recr table of report.sso, SS3 presents the annual R0 and h implies by the AB parameters and that year's biology. These values are reported, but are not used in calculations. The values gradually increases to ln(R0) = 11.7566 and steep = 0.60523. Same output for the AB run and the R0h run.

Add to the estimation run the possibility that a spawn_recr parameter is time varying with a block change in 1985. This is before the growth change. Set the block offset to be 0 so that R0 does not change, but the code logic sees that it could change. This invokes code logic that updates SPR0 annually for use in the SRR calculations. FOr the AB run, do the block on the beta parameter; again with no change realized.

With Bmark_yr set to 1971:
#_Mgmt_Quantity
SSB_unfished 50174.3 0
Totbio_unfished 124677 0
SmryBio_unfished 124187 0
Recr_unfished 6644.32 0
SSB_Btgt 17159.6 0
SPR_Btgt 0.488397 0
annF_Btgt 0.0857971 0
Dead_Catch_Btgt 2786.1 0
SSB_SPR 11455.2 0
annF_SPR 0.118586 0
Dead_Catch_SPR 2594.79 0
SSB_MSY 17243 0
SPR_MSY 0.489689 0
annF_MSY 0.0853933 0
Dead_Catch_MSY 2786.13 0
Ret_Catch_MSY 2786.13 0
B_MSY/SSB_unfished 0.343662 0

and long-term projection, which always uses time-varying biology, shows:
SSB_2079 31493.1 0
SSB_2080 31499.3 0
SSB_2081 31505 0
Recr_2079 5869.68 0
Recr_2080 5870.04 0
Recr_2081 5870.37 0
F_2079 0.0853933 0
F_2080 0.0853933 0
F_2081 0.0853933 0
OFLCatch_2079 4312.7 0
OFLCatch_2080 4313.52 0
OFLCatch_2081 4314.29 0

Note that depletion (Bratio) is using Bmsy as denominator:
Bratio_2079 1.82643 0
Bratio_2080 1.82679 0
Bratio_2081 1.82712 0

Now change Bmark_yr to 2001, a year for which biology is still changing:
SRR parms are unaffected:
SR_LN(R0) 8.80152 0 Fix 0.207197
SR_BH_steep_derived 0.529114 0 Fix 0.411393
SR_BH_ln(alpha) -0.518856 0 Act 0.413524
SR_BH_ln(beta) -9.57203 0 Act 0.186767

#_Mgmt_Quantity shows bigger biomass and catch, but note that biology is not in equilibrium
SSB_unfished 75560 0
Totbio_unfished 147662 0
SmryBio_unfished 147172 0
Recr_unfished 6644.32 0
SSB_Btgt 25841.5 0
SPR_Btgt 0.413648 0
annF_Btgt 0.102823 0
Dead_Catch_Btgt 4318.95 0
SSB_SPR 24515.1 0
annF_SPR 0.108106 0
Dead_Catch_SPR 4315.26 0
SSB_MSY 25656.1 0
SPR_MSY 0.411741 0
annF_MSY 0.103543 0
Dead_Catch_MSY 4319.04 0
Ret_Catch_MSY 4319.04 0
B_MSY/SSB_unfished 0.339547 0

long-term projections are larger that MSY because growth continues to increase. Here is bodywt at age for some relevant years:
START YEAR THRU yr before Linf increase
1989: 0.0737219 0.192683 0.369545 0.593713 0.851289 1.12842 1.41297 1.69524 1.96798 2.22616 2.46661 2.68765 2.8887 3.06997 3.23223 3.3766 3.50441 3.61708 3.71605 3.80272 3.87842 3.9444 4.0018 4.05166 4.0949 4.13237 4.1648 4.19285 4.2171 4.23804 4.25611 4.27171 4.28516 4.29677 4.30677 4.31539 4.32281 4.32921 4.33472 4.33947 4.3464

last year of timeseries and basis for MSY caLCS
2001: 0.0737219 0.220287 0.452707 0.758356 1.11794 1.51115 1.91966 2.32846 2.72613 3.10458 3.45852 3.78493 4.08247 4.2682 4.47167 4.6518 4.81045 4.94959 5.07121 5.1772 5.26935 5.34931 5.4186 5.47857 5.5304 5.57518 5.61384 5.64719 5.67595 5.70075 5.72212 5.74054 5.7564 5.77007 5.78184 5.79197 5.8007 5.80821 5.81468 5.82025 5.84089

10th year of long projection
2010: 0.0737219 0.220287 0.452707 0.758356 1.11794 1.51115 1.91966 2.32846 2.72613 3.10458 3.45852 3.78493 4.08247 4.35103 4.59138 4.80493 4.99348 5.15912 5.304 5.43029 5.54008 5.63531 5.69266 5.75383 5.80665 5.85223 5.89155 5.92544 5.95465 5.97981 6.00149 6.02016 6.03624 6.05008 6.062 6.07226 6.08109 6.0887 6.09524 6.10088 6.1218

last year of long projection
2081 0.0737219 0.220287 0.452707 0.758356 1.11794 1.51115 1.91966 2.32846 2.72613 3.10458 3.45852 3.78493 4.08247 4.35103 4.59138 4.80493 4.99348 5.15912 5.304 5.43029 5.54008 5.63531 5.71776 5.78907 5.85066 5.90383 5.94968 5.98922 6.0233 6.05266 6.07795 6.09973 6.11848 6.13463 6.14853 6.1605 6.17081 6.17967 6.18731 6.19388 6.21399

The difference between the MSY result and the long-term projection result is an indication of the outstanding need for the capability to use equilibrium growth in the MSY calculations.

Note that I remembered to turn off the control rule so projections are at FMSY, not reduced for ABC

SSB_2079 26191.9 0
SSB_2080 26196.5 0
SSB_2081 26200.7 0

Recr_2079 5519.84 0
Recr_2080 5520.18 0
Recr_2081 5520.49 0

F_2079 0.103543 0
F_2080 0.103543 0
F_2081 0.103543 0

Bratio_2079 1.02088 0
Bratio_2080 1.02106 0
Bratio_2081 1.02122 0

OFLCatch_2079 4379.24 0
OFLCatch_2080 4379.99 0
OFLCatch_2081 4380.68 0

note that MSY catch is 4319, so slightly smaller because more body growth in the projection.

This section of the forecast-report.sso output has a bit more information, especially the recruitment level. I have created a new issue for adding recruit to the standard management quantity report.:

find_Fmsy_to_maximize_dead_catch
SPR@MSY 0.411741
Fmult 0.103543
ann_F 0.103543
Exploit(Catch/Bsmry) 0.0717375
Recruits@MSY 5479.32
SPBmsy 25656.1 4.68236
SPBmsy/SPB_virgin 0.511341
SPBmsy/SPB_unfished 0.339547
MSY_for_optimize 4319.04 0.788245
MSY_encountered 4319.04 0.788245
MSY_dead 4319.04 0.788245
MSY_retain 4319.04 0.788245
MSY_revenue 4319.04
MSY_cost 0
MSY_profit 4319.04
Biomass_Smry 60206.2 10.9879

The time_vary SR_parm flag causes SS3 to use the current year's biology to recalc SPR0 and use that SPR0 in executing the R = F(SSB) code. This is not routine!!!
This flag is triggered when R0 or h is timevarying, or alpha or beta when using that formulation.
With time-varying SPR0, the estimate of h is 0.55 versus a value of 0.49 when SPR0 is held constant at the start year value (normal approach). The logL between the two model runs is similar (1367.23 vs 1367.46), even though nearly all parameters show some shifting. This demonstrates a small dependence of h on the value of SPR0.
Both of the above runs, used styr biology for MSY calcs.

With the AB approach, and not using timevary beta, the logL is 1367.46.
with AB approach and with beta having timevary flag triggered, the same result is obtained because the alpha-beta parameterization does not depend on SPR0.

Note that time_vary SR_parm flag is nearly never used in assessments. All assessments used the start year biology's SPR0 when executing the spawner-recruitment relationship, so produce identical results as when the AB approach is used. Note also that both the R0h and the AB approaches use the current year's fecundity when calculating the spawning biomass; there is no difference in that regard.

[Rick-Methot-NOAA]

@timjmiller @liz-brooks
Here is experiment to show the effect of a change in fecundity units.
Base run has fecundity = 1.0, so reproductive output is same as mature female biomass.
Experiment run sets fecundity to 0.1
all parameters get updated in the model run.

Base:
SPR0_(virgin): 7.44189 #_uses_biology_from_start_year: 1971
Ln(alpha): -0.510467 alpha 0.600215
Ln(beta): -9.62761 beta 6.58841e-05
steepness_derived: 0.527563
ln(R0)_derived: 8.7637

Experiment:
SPR0_(virgin): 0.744189 #_uses_biology_from_start_year: 1971
Ln(alpha): 1.79203 alpha 6.00162
Ln(beta): -7.32514 beta 0.000658766
steepness_derived: 0.527541
ln(R0)_derived: 8.7637

Conclusion:
SPR0 changes 10x as expected.
parameters alpha and beta get much different estimated values
model time series results and -logL were identical
h and R0 derived from the estimated alpha and beta were identical; e.g. unaffected by the rescaling of SPR0.

[Rick-Methot-NOAA]

update_SRR_AB.zip

[liz-brooks]

hi rick, alpha and beta both changed by a factor of 10 (10 times larger for Experiment), while spr0 got 10 times smaller for Experiment. the fact that h and R0 are unchanged is expected because alpha or beta always multiplies spr0 when you translate from the alpha beta parameterization to the h R0 parameterization: h = alpha*spr0 / (4+alpha*spr0) R0 = (alpha*spr0 -1)/beta*spr0 the thought experiment in section 3 of miller and brooks walks through what happens if you have a simple scalar change within a single assessment time series (i.e., not scaling the whole time series and comparing two separate runs).

…

On Wed, Jun 26, 2024 at 5:52 PM Richard Methot ***@***.***> wrote: @timjmiller <https://github.com/timjmiller> @liz-brooks <https://github.com/liz-brooks> Here is experiment to show the effect of a change in fecundity units. Base run has fecundity = 1.0, so reproductive output is same as mature female biomass. Experiment run sets fecundity to 0.1 all parameters get updated in the model run. Base: SPR0_(virgin): 7.44189 #_uses_biology_from_start_year: 1971 Ln(alpha): -0.510467 alpha 0.600215 Ln(beta): -9.62761 beta 6.58841e-05 steepness_derived: 0.527563 ln(R0)_derived: 8.7637 Experiment: SPR0_(virgin): 0.744189 #_uses_biology_from_start_year: 1971 Ln(alpha): 1.79203 alpha 6.00162 Ln(beta): -7.32514 beta 0.000658766 steepness_derived: 0.527541 ln(R0)_derived: 8.7637 Conclusion: SPR0 changes 10x as expected. parameters alpha and beta get much different estimated values model time series results and -logL were identical h and R0 derived from the estimated alpha and beta were identical; e.g. unaffected by the rescaling of SPR0. — Reply to this email directly, view it on GitHub <#191 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACFL5LN6RKHHZNCQ2EW4OTLZJMZ3JAVCNFSM5ADSFN2KU5DIOJSWCZC7NNSXTN2JONZXKZKDN5WW2ZLOOQ5TEMJZGI3DQNJXGUZA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Liz Brooks, PhD Operations Research Analyst Population Dynamics Branch NOAA/NMFS Northeast Fisheries Science Center 166 Water Street phone: 508.495.2238 Woods Hole, MA 02543 fax: 508.495.2393 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

[Rick-Methot-NOAA]

After much, much staring at the code. I find that the call to Equil_Spawn_Recr_Fxn will not always send the correct SSB0 and R0.
@iantaylor-NOAA
Notes below based on code as found in 3.30.22.1

Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work(2), SR_parm_work(3), SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR

FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const prevariable& SRparm2, const prevariable& SRparm3,
const prevariable& SSB_virgin, const prevariable& Recr_virgin, const prevariable& SPR_temp)

Problem with above pairing is that Equil_Spawn_Recr_Fxn needs to access the real virgin values unless the SR_update_SPR0 flag is set due to time-varying
spawner-recruitment parameters. If only biology is time-varying, the the spawn-recruit equations should use the virgin values.
So, the call to Equil_Spawn_Recr_Fxn should have a switch such that the correct SPR0 will be used.

For the Spawn_Recr function which is called during time series and forecast, the pairing is:
Recruits = Spawn_Recr(SSB_use, R0_use, SSB_current); // calls to function Spawn_Recr

FUNCTION dvariable Spawn_Recr(const prevariable& SSB_virgin_adj, const prevariable& Recr_virgin_adj, const prevariable& SSB_current)

where SSB_use and R0_use are set to virgin or are updated depending on the SR_update_SPR0.
So, if only biology is time-varying, the Spawn_Recr correctly uses the virgin values
Repeat: if SR_update_SPR0 is set due to time-vary SRparms, then SSB_use will use those parameters and the updated biology.

To align the terminology:
For the calling routines,
SSB_use for the call to Spawn_Recr and SSB_unf for the call to Equil_Spawn_Recr_Fxn serve the same function. They should be merged to one name.

For the receiving functions,
SSB_virgin and SSB_virgin_adj serve the same role for Equil_Spawn_Recr_Fxn and Spawn_Recr respectively.

[iantaylor-NOAA]

@Rick-Methot-NOAA, I haven't followed the progress on the stock-recruit revisions very closely, but I think your explanation makes sense. However makes sense to you to sort all this out will be fine with me.

[Rick-Methot-NOAA]

I determined that we need fuller control over when the SSBpR0 is updated and used. Therefore, I have re-purposed an used placeholder in the Spawner-recruitment input to provide this control:

1 #  SR_update_SSBpR0
#  0 - OK, but only if no timevary biology or SR parm
#  1 - best: update SSBpR0 for benchmark and for time series only if SRparm R0 or h (not regime) is set to have time-varying property
#  2 - incorrect (old, incorrect SS3 approach):  always update SSBpR0 for benchmark's use of spawner-recruitment, but only for the time series if there is a timevary SR parm
#  3 - option:  do not update SSBpR0 (do keep start year SSBpR0), even if R0 or h is set to have time-varying property

[Rick-Methot-NOAA]

add control for basis of HCR inflection:

1 # Control rule method (0: none; 1: ramp does catch=f(SSB), buffer on F; 2: ramp does F=f(SSB), buffer on F; 3: ramp does catch=f(SSB), buffer on catch; 4: ramp does F=f(SSB), buffer on catch) 
# values for top, bottom and buffer exist, but not used when Policy=0
0.4 # Control rule inflection for constant F (as frac of Bzero, e.g. 0.40); must be > control rule cutoff, or set to -1 to use Bmsy/SSB_unf 
# Also see HCR_anchor below
0.1 # Control rule cutoff for no F (as frac of Bzero, e.g. 0.10) 
0.75 # Buffer:  enter Control rule target as fraction of Flimit (e.g. 0.75), negative value invokes list of [year, scalar] with filling from year to YrMax 
#
3 #_N forecast loops (1=OFL only; 2=ABC; 3=get F from forecast ABC catch with allocations applied)
3 # First forecast loop with stochastic recruitment
1 # Forecast base recruitment:  0= spawn_recr; 1=mult*spawn_recr_fxn; 2=mult*VirginRecr; 3=deprecated; 4=mult*mean_over_yr_range
# for option 4, set phase for fore_recr_devs to -1 in control to get constant mean in MCMC, else devs will be applied
1 # multiplier on base recruitment 
2 # HCR_anchor: 0 or 1 uses unfished benchmark SSB (old hardwired approach), 2 = virgin SSB