Implementing new recruitment dev vector approach

[k-doering-NOAA]

Imported from redmine, Issue #57889 Opened by @RickMethot on 2018-11-23 Status when imported: In Progress

create new recruitment deviation vector approach do_recdev==3

With this approach: exp_rec(y,4) = R0_adj * mfexp(dev(y)) where R0_adj is R0 adjusted for blocks, etc. on R0; and dev(y) is a simple vector of unconstrained deviations

Then in the objective function: rather than getting the -logL from the dev vector, now create a new deviation which is dev=log(exp_rec(y,3)/exp_rec(y,4)); where exp_rec(y,4) is shown above and exp_rec(y,3) is the expected value of recruitment adjusted for spawn_recr curve, regimes and env, and bias-adjustment.

Initial test shows that the -logL is about same as using do_recrdev==2 which is simple devs in which exp_rec(y,4) = exp_rec(y,3) * mfexp(dev(y))

[k-doering-NOAA]

comment from @RickMethot on 2018-11-23: adding a constraint such that the (sum(devs))^2 contributes to the total logL give result identical to the result obtained from do_recdev==1 (which invokes the ADMB zero-centered dev_vector). But with this constraint the convergence takes many more iterations. Will try to tweak this.

[k-doering-NOAA]

comment from @RickMethot on 2018-11-27: completed. Now have approach 3 and 4 that do the above and with approach 4 creating a penalty based on sum(devs)

do_recdev: 0=none; 1=devvector (R=F(SSB)+dev); 2=deviations (R=F(SSB)+dev); 3=deviations (R=R0*dev; dev2=R-f(SSB)); 4=like 3 with sum(dev2) adding penalty"

Performance is as expected and documented in attached spreadsheet. Need test now with MCMC

add output of N iterations to ss_summary add output of sum(recdevs) to ss_summary

[k-doering-NOAA]

comment from @RickMethot on 2018-11-27: Implementation of this feature required separating the recruitment FUNCTION into one function that does Rhat=F(SSB) and a separate function that applies the dev to Rhat

[k-doering-NOAA]

comment from @iantaylor-NOAA on 2018-11-29: The new options 3 and 4 don't work with recdevs starting prior to startyr, as indicated by the following crash associated with the attached files (which went away when I shifted the start of the early devs from 1946 to 1966. Some mcmc comparisons are in process, but figured I would report this issue before I forget to do so.

Estimating...please wait...
0 0 -log(L): 0  between
 MG setup OK
 growth1 OK
 natmort OK
 migr OK
 selex OK, ready to call ALK and fishselex
 ready for virgin age struc
 ready for initial age struc
 OK with initial conditions
 OK with time series
 did survey obj_fun  0 4.48942
 did agecomp obj_fun  6946.59 419.773
 initequ_catch -log(L) 0
 catch -log(L)  2.53129e-17 0
matrix row bound exceeded -- row index 1946 is lower than 1964 in linad99\fvma_acc.cpp:27.

[k-doering-NOAA]

comment from @iantaylor-NOAA on 2018-12-18: I've slowly been able to do some digging into this issue.

First, it seemed to work to recompile ADMB after commenting out the test for mc_phase in the third line below from nh99/model5.cpp as suggested by Johnoel.

    dvariable s=mean(*this); 
    pen+=10000.0*s*s; 
    if (!initial_params::mc_phase) 
    { 
      (*this)-=s; 
    } 

I think this fixes things based on no longer seeing a difference between the values stored in the binary ss.psv during the -mcmc run from those in posteriors.sso created during -mceval. There seemed to be some indication from Dave Fournier that this penalty keeping a dev_vector zero-centered caused problems in the MCMC, but I don't see any issue with the hake model that I was testing with. The MCMC takes a similar amount of time to run and has a similar acceptance rate.

The difference in overall results is also thankfully very small. The sum of the main recdevs as reported in ss.psv shows a range of the -0.023 to 0.024, definitely bigger than the range -1.481e-06 to 9.638e-07 for the values in posteriors.sso, but too small to have much impact on the model. I'm not sure why these results seem more similar than the model that uses simple (non-zero-centered) devs, so should look deeper into that.

Even if ADMB works with the change (I need to report to Johnoel and/or the ADMB issues list), it still may be worth pursuing a way to control these features within SS. To that end, I think the calculation that would make the results match the current dev structure would be

dev=log(exp_rec(y,4)/exp_rec(y,2));

instead of

dev=log(exp_rec(y,3)/exp_rec(y,4));

And it also seems like the bias adjustment is perhaps not being accounted for. This is based on simply trying to get the "dev" column in the recruitment output to match a model that uses the old approach.

[k-doering-NOAA]

comment from @RickMethot on 2018-12-20: Ian,

I think the matching of the dev vector is just a reporting issue. What gets reported in the spawn_recr section is the dev vector itself, not the new dev that is result of the dev=log(exp_rec(y,3)/exp_rec(y,4)) calculation.

Also, I have an active conversation with Tim Miller regarding time-varying biology. I just realized that with this new approach we can base the SR curve based on the average biology over a range of years because the obj_fun is not called until after the whole time series is processed.

[k-doering-NOAA]

comment from @k-doering-NOAA on 2020-04-01: @richard.methot, can you also add an update on this issue? Is it left open because it needs to be documented or because you plan to write an essay on it (or both?)

[k-doering-NOAA]

comment from @RickMethot on 2020-04-01: Kelli is developing a plan to test performance of the new alternatives, then will document with an essay.

[kellijohnson-NOAA]

Question for @Rick-Methot-NOAA regarding do_recdev == 2 How come the estimate of the recruitment deviation in the terminal year (regardless if it is a main or late rec dev) is zero. This is not the case when do_recdev == 1. Chantel and I tested it with a few models and with everything set exactly the same, when you change from simple to the constrained vector the terminal year estimate becomes non-zero.

[Rick-Methot-NOAA]

I suspect this is because there are no data to inform that last dev, so it is 0.0 if part of an unconstrained vector (option 2), but with the constrained vector ADMB is free to move it to achieve the sum to zero constraint.

Rick

On Tue, Mar 2, 2021 at 7:57 AM Kelli Johnson notifications@github.com wrote:

Question for @Rick-Methot-NOAA https://github.com/Rick-Methot-NOAA regarding do_recdev == 2 How come the estimate of the recruitment deviation in the terminal year (regardless if it is a main or late rec dev) is zero. This is not the case when do_recdev == 1. Chantel and I tested it with a few models and with everything set exactly the same, when you change from simple to the constrained vector the terminal year estimate becomes non-zero.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/nmfs-stock-synthesis/stock-synthesis/issues/37#issuecomment-789012398, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABPV4IGRS5DWJA6HN4FYJ3LTBUDHNANCNFSM4TLVLRJA .

[Rick-Methot-NOAA]

closely related to #29

[k-doering-NOAA]

Simulation testing and documentation will take place in #29

[Rick-Methot-NOAA]

adding the comparison spreadsheet compare_recdev_approaches.xlsx