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:

• [x] add SRR function 10 to add a alpha, beta parameterization
• [x] substantially refactor the variables names to clarify SPR and SSB/R usages
• [x] add input in SRR section to allow for specification regarding biology used in SSBpR0
• [x] 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
• [x] 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

[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 https://github.com/nmfs-stock-synthesis/stock-synthesis/issues/191#issuecomment-1429004201, 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]

R0h and AB are looking identical except when R0 parameter is time-varying. Basically R0h and AB both define parameters that are relevant to virgin (year 1) biology, then implement the equation using current SSB as the driving factor. Note that the use of the regime offset parameter does not make R0 time-varying.

When R0 is time-varying, then SS3 uses the current year's biology to calculate an equilibrium SSB and uses it in place of the virgin SSB that is normally used in the R0h formula. So, time-varying biology only adjusts SSB_virgin_adj if R0 is time-varying.

Next step is to create some scenarios to demo what is happening above. constant R0 and constant biology constant R0 and timevary bio time-vary R0 and constant bio time-vary R0 and constant bio. For each condition, implement with R0h and AB.

[timjmiller]

Below applies for Beverton-Holt or Ricker as defined in the paper Liz and I did. Liz, please correct any errors.

The R0/steepness parameterization requires a third variable: unfished SSB/R, which varies when post-recruit biology varies. R0, S0 (equilibrium unfished R and SSB) and steepness are functions of equilibrium unfished SSB/R and therefore will also vary for the same reason. The a/b (pre-recruit biology: density independent and dependent mortality between egg and recruitment) parameterization requires no information about post-recruit biology. Therefore, the shape of the SR curve should only change when a or b changes. Post-recruit biology is only needed to define the equilibrium points on the curve (at F=0 or F=Fmsy, etc.) and changes in post-recruit biology only change the equilibrium points on the curve.

However, variation in a and b will also cause variation in equilibrium points because steepness (a only) and R0 are functions of these as well as equilibrium unfished SSB/R.

If R0 is time varying because a or post-recruit biology is varying then steepness should also be varying because steepness is also a function of these. If R0 is varying only because b is varying, then steepness would still be constant because steepness is not a function of b. If both pre- and post-recruit biology (a,b, and unfished SSB/R) are constant then R0 and steepness are constant. It isn't possible to have constant R0 and time-varying post-recruit biology unless variation in pre-recruit biology (a,b) occurred in a way that was completely correlated with the variation in post-recruit biology.

I don't think we typically model a,b as time varying in assessment models, but there are several papers that estimate state-space SR models with autoregressive parameters.

On Mon, May 20, 2024 at 5:52 PM Richard Methot @.***> wrote:

R0h and AB are looking identical except when R0 parameter is time-varying. Basically R0h and AB both define parameters that are relevant to virgin (year 1) biology, then implement the equation using current SSB as the driving factor. Note that the use of the regime offset parameter does not make R0 time-varying.

When R0 is time-varying, then SS3 uses the current year's biology to calculate an equilibrium SSB and uses it in place of the virgin SSB that is normally used in the R0h formula. So, time-varying biology only adjusts SSB_virgin_adj if R0 is time-varying.

Next step is to create some scenarios to demo what is happening above. constant R0 and constant biology constant R0 and timevary bio time-vary R0 and constant bio time-vary R0 and constant bio. For each condition, implement with R0h and AB.

— Reply to this email directly, view it on GitHub https://github.com/nmfs-ost/ss3-source-code/issues/191#issuecomment-2121270716, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEIGN7H6CN3UG6GHTCVPD3DZDJWBTAVCNFSM5ADSFN2KU5DIOJSWCZC7NNSXTN2JONZXKZKDN5WW2ZLOOQ5TEMJSGEZDOMBXGE3A . 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

[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 https://github.com/nmfs-ost/ss3-source-code/issues/191#issuecomment-2123485783, 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 https://github.com/nmfs-ost/ss3-source-code/issues/191#issuecomment-2123485783, 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 https://github.com/nmfs-ost/ss3-source-code/issues/191#issuecomment-2129958687, 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 = alphaspr0 / (4+alphaspr0) R0 = (alphaspr0 -1)/betaspr0

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 https://github.com/nmfs-ost/ss3-source-code/issues/191#issuecomment-2192685752, 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