[Feature]: Mean recruitment adjustment

[Rick-Methot-NOAA]

Describe the solution you would like.

see discussion in #593
When using recdev option 1, the dev_vector with tuned bias adjustment forces the arithmetic mean recruitment of the time series to be consistent with the mean from the spawner-recruitment relationship.
This is because the spawner-recruitment relationship is used to produce MSY; MSY is not based on the time series of estimated recruitments.
With the transition to recommending that recdevs be based on recdev option = 2, which is not zero-centered, there is a concern that the MSY calculations will not be consistent with the mean of the recruitments.
It seems that MSY should remain based on the SRR, as adjusted by the mean deviation of the recdevs.
Solution 1: add a penalty to mean recdev to improve degree to which it is 0.0
Solution 2: base the benchmark quantities on log(R0) - mean(dev)

Describe alternatives you have considered

Rely on the sigma penalty to keep mean recruitment close enough to zero. However, it is clear that his does not always happen. A change to log(R0) can be perfectly offset by a drift in mean recruitment.

Statistical validity, if applicable

No response

Describe if this is needed for a management application

definitely important to not allow a know bias in this calculation

Additional context

No response

[Cole-Monnahan-NOAA]

Isn't option 1 just recreating the dev_vector, but more explicitly? Perhaps a better alternative for this option is to get the dev_vector working better in ADMB in general? Seems unlikely given the ADMB dev status.

I've always thought option (2) made the most sense. It should mimic the dev_vector approach in both the estimates and derived quantities, with no apparent downsides that I can anticipate. @jimianelli and I were discussing this recently. Any thoughts Jim?

[Rick-Methot-NOAA]

In principle yes, but some situations cause the mean of the devs to drift away from zero.

[Cole-Monnahan-NOAA]

Maybe I don't understand (likely) but I don't see why it matters whether there is a drift or not. Can we not consider the two approaches (dev_vector vs. post-hoc centering) just a different way of implementing the same feature? Or as you put it "A change to log(R0) can be perfectly offset by a drift in mean recruitment."

How can I help with this issue? Do some analyses need to be done?

[allanhicks]

I'm very glad to see this issue being discussed. I understand that a dev vector is an issue in MCMC when using ADMB, but when the recruitment deviates do not sum to zero, the benchmark calculations do not make sense (as Rick explained earlier). The CASAL manual <https://niwa.co.nz/sites/default/files/sites/default/files/casalv230-2012-03-21.pdf> (page 32) describes this and some solutions that were proposed many years ago. The current Pacific hake model uses option 2 and the sum of the deviates is close to zero, but not exactly zero. My concern is that other assessments with potential changes in average recruitment and lack of data or short "main" periods of recruitment may result in a large difference in average recruitment compared to R0. One test would be to simulate a model many years into the future to determine the benchmark quantities (B0, MSY, etc.). I believe that without a sum to zero recruitment vector, you cannot be guaranteed to achieve those benchmarks. I think that when the AFSC assessments assume steepness of 1 and use external projection software, this is not a concern. Nevertheless, I think some testing would be useful. I have not experimented with other recruitment types in SS, but would options 3 or 4 be worth considering? I'll try to cook up some examples to show where this is a concern, but if I can be of any help, please let me know.

…

On Mon, Apr 29, 2024 at 4:57 PM Cole Monnahan ***@***.***> wrote: Maybe I don't understand (likely) but I don't see why it matters whether there is a drift or not. Can we not consider the two approaches (dev_vector vs. post-hoc centering) just a different way of implementing the same feature? Or as you put it "A change to log(R0) can be perfectly offset by a drift in mean recruitment." How can I help with this issue? Do some analyses need to be done? — Reply to this email directly, view it on GitHub <#594 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABPA4ZNMJMTE67RVEFOKBJDY73M6XAVCNFSM6AAAAABG3RSILSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDAOBTHA4TIMJXG4> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[Rick-Methot-NOAA]

Thanks Allan,
My proposal is to modify the -logL for recdevs to include this term:
square((mean_dev - 0.0) / se))
where se is sigmaR / N.

current code (w/o rho) is:
recr_like = sd_offset_rec * log(sigmaR); // where sd_offset_rec is N adjusted for bias adjustment = sum(biasadj)
recr_like += norm2(recdev(recdev_first, recdev_end)) / two_sigmaRsq;

Note that all biasadj set to 1.0 for MCMC, so sd_offset_rec equals N years when going to MCMC!

[allanhicks]

That sounds reasonable. I imagine that N is the number of years in the main period. Should the se = sigmaR/sqrt(N)? It will be interesting to see what the behavior is and if there are still cases where this penalty will not be strong enough to keep the mean near zero. Dividing sigma by N means that the se is larger at smaller sample sizes (as expected), but does that mean that with a shorter main period this penalty is less influential and could more easily depart from zero? I may be heading towards some recommendation of the minimum length of the main period of recruits (which may already be specified). This seems like a good compromise between simple random deviates and a deviate vector forced to sum to zero.

…

On Tue, Apr 30, 2024 at 8:23 AM Richard Methot ***@***.***> wrote: Thanks Allan, My proposal is to modify the -logL for recdevs to include this term: square((mean_dev - 0.0) / se)) where se is sigmaR / N. current code (w/o rho) is: recr_like = sd_offset_rec * log(sigmaR); // where sd_offset_rec is N adjusted for bias adjustment = sum(biasadj) recr_like += norm2(recdev(recdev_first, recdev_end)) / two_sigmaRsq; Ah-Ha: sd_offset_rec should change to equal N years when going to MCMC! — Reply to this email directly, view it on GitHub <#594 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABPA4ZJ6XDGL3WFO5WEHXX3Y76ZN5AVCNFSM6AAAAABG3RSILSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDAOBVGY2TINBSHE> . You are receiving this because you commented.Message ID: ***@***.***>

[Cole-Monnahan-NOAA]

I don't know CASAL at all but from the manual it seems like they already implement Rick's option 2.

Here YCS is year class strength which when put into log space would be subtracting the mean.

Do I have that right?

@Rick-Methot-NOAA can you clarify option 2? I thought this would not change anything in the likelihoods or other calculations, and instead only impact benchmark calculations (presumably all in the forecast module?) by defining a new "mean recruitment" as lnR0-mean(devs) where the parameters would be estimated using the simple deviations. In other words, it wouldn't affect estimation, only post-hoc calculations. Am I missing something here?

[Rick-Methot-NOAA]

We're getting close...
Either:
(a) implement a penalty on mean(recdev) such that it is close enough to zero and we treat the remaining delta as noise;
or
(b) accept that with recdev option 2 there could be drift in mean(recdev), so create an adjusted log(R0)' = log(R0) + mean(recdev) to use in benchmark and forecast

We'll need to worry about how this plays with bias_adjustment when transitioning from MLE to MCMC

[Cole-Monnahan-NOAA]

With (b) can you confirm (intuitively at least -- we should check) that you should get the same benchmarks using dev_vectors and simple vectors when implementing the new option 2? In other words, will this break backwards compatibility?

[Rick-Methot-NOAA]

I would implement this as a new, not revised, option so we could test.
Logically, it should be the same as dev_vector.
Where it gets squirrely is in calculation of initial population size and then annual depletion, e.g. SSB(y)/SSB0. So I think that the offset needs to be applied there also to achieve the current total internal consistency of the SS3 time series.
Of course, this is just cementing further the SS3 dependence on the R0,h formulation of spawner-recruitment. There is another branch in which I explore the a,b formulation per Miller and Brooks. Which is totally feasible, except for reworking all of the depletion dependencies. I think I can make that work. Basically, the a,b approach lets you calculate the R0,h conditioned on start year biology then I need to work out what happens when biology changes during time series. If you have comments, please put them at #191