review implementation of parameter random effects

[k-doering-NOAA]

Imported from redmine, Issue #39714
Opened by @RickMethot on 2015-10-29
Status when imported: In Progress

In 3.30, the deviations are assumed to be unit normal and are multiplied by a stddev before being used. However, they may not actually be unit normal in any instance and this caused consternation for at least one user, Grant Thompson.

[k-doering-NOAA]

comment from @RickMethot on 2019-02-13:
Also, need to check to see if the SD_offset is being used for parameter devs

[k-doering-NOAA]

comment from @RickMethot on 2020-08-11:
Find two improvements to make in the usage of parameter deviations:

1. the dev vectors were created with min-max bounds of -10 and 10. But the check on the current value relative to those bounds was set at -5, 5. So if a run produced large deviations, then when that ss.par was read back into ss, those large values would get truncated. Rectify by putting check at -10 and 10.

2. In write_report, the routine for checking parameters close to bounds was not being called for parameter deviations. Rectify this with a small bit of new code.

prepare these changes for 3.30.16.

[k-doering-NOAA]

comment from @RickMethot on 2020-08-12:
I also did some extensive experimenting and testing which results in the following observations:

1. in 3.24, the devs were in parameter units so were divided by the stddev value and the -logL included a component for -log(stddev). With that approach, if the stddev parameter was allowed to be estimated, then it tended to go towards 0.0 and -log(stddev) overwhelms the overall logL.
2. In 3.30, the devs are designed as unit normal, so the parameter change uses dev(y)*stddev. The -logL does not use -log(stddev) because that term is 0.0. With this approach, if stddev is allowed to be estimated, the value of stddev gets very large because there is not penalty on it, and the rmse of the estimated devs gets small. These small devs produce very low contribution to the logL, and these small devs times a large stddev give SS maximum flexibility to fit the data.
3. I experimented with addition of a penalty: square(10*(1-rmse)) which tends to keep rmse near value of 1.0 and enables estimation of the stddev now conditioned on the devs having rmse near 1.0. Works, but convergence is horrible because there is extreme high correlation among the devs with this constraint.
4. Also did testing of the various dev methods. Find that both method 4 (mean-reverting random walk) and new method 5 (MRRW with logit transform) work quite well. Method 5 seems somewhat better, but still testing.

[k-doering-NOAA]

comment from @RickMethot on 2020-08-12:
extending the last message:
5. It would be good to get a test of the ability to estimate the stddev parameter using MCMC w/ NUTS. However, still see nothing about the design of the approach that inherently will give a vector of devs with rmse near 1.0 and a stddev that scales those devs to achieve the desired results. I think we should consult with Tim Miller to see how he implements random effects in WHAM.
6. For now, with MLE in SS, I think the best approach is to provide a fixed stddev value such that baseparm +- 2.0*stddev will provide a reasonable range for the timevarying parameter.

[k-doering-NOAA]

comment from @iantaylor-NOAA on 2020-08-12:
Good work Rick making progress on this.

Unfortunately I don't know random effects well enough to comment on this, so I support talking to Tim or others with expertise.
I did just look at the WHAM code and see calculation of time-varying and age-dependent M here: https://github.com/timjmiller/wham/blob/master/src/wham_v0.cpp#L434-L483
where Sigma_M is the random effect variance based on notes at https://github.com/timjmiller/wham/blob/master/src/wham_v0.cpp#L139-L142.
On first glance it does not appear that WHAM includes a generalized function for random effects that can be applied to different processes.

[k-doering-NOAA]

comment from @RickMethot on 2020-08-12:
Thanks Ian. I took a look and the logL of the dev vector seems wrapped up
in the Laplace transformation and is a bit opaque to me yet.

Here is the example I am using and the similar results using MRRW and the
logit version of MRRW.
Both have stddev set to a large value 14 and fixed.
logit version produces lower logL because its devs in log space have
smaller absolute value, so closer to 0.0, so have logL that is smaller. But
these same devs produce a slightly better fit to the agedata.
Can you give this a try using NUTS and with estimation of the stddev
parameters?

Rick

Richard D. Methot Jr. Ph.D.NOAA Fisheries Senior Scientist for Stock
Assessments

Mobile: 301-787-0241

On Wed, Aug 12, 2020 at 4:28 PM vlab.redmine@noaa.gov wrote:

[k-doering-NOAA]

comment from @iantaylor-NOAA on 2020-08-13:
I need to get back up to speed on NUTS first, but I will try to work on this sometime tomorrow.

[k-doering-NOAA]

comment from @iantaylor-NOAA on 2020-08-14:
The model has been running with adnuts for about 4 hours using ss.exe compiled from the current repository (and ADMB 12.2).
I didn't think about how long it might take. There's an option for specifying a duration which in hindsight I could have used, but instead will just let it continue to run over the weekend. I'm on out Monday but will report on results on Tuesday.

Note: the model expected a wtatage.ss file so I used the attached wtatage_used_by_ian.ss file from a recent Simple model assuming that it will be similar enough to have little impact on this test.

[k-doering-NOAA]

comment from @iantaylor-NOAA on 2020-08-18:
When I kicked off the NUTS algorithm on Friday I failed to re-read the message and note the suggestion to include "estimation of the stddev". Of course that's the whole point of exploring Bayesian sampling in this case. Also, it took 33 hours to complete the first of three chains, so I stopped it part-way through chain 2 anyway.

I've now just started a new NUTS sampling with 1 chain and using the "duration" argument to limit the length to 24 hours so I can report on something at this time tomorrow.
I changed the AgeSel_P2_FISHERY1(1)_dev_se and AgeSel_P3_FISHERY1(1)_dev_se parameters to have MAX = 2.0 and INIT = 0.5 (as shown below) out of concern that sampling really high values will lose signal about the dev_se and make it take longer. I can try again with wider bounds if I'm misguided on what a reasonable range should be.

LO      HI      INIT    PRIOR   PR_SD   PR_type PHASE   #       parm_name
0.0001  2       0.5     0.5     0.5     0       6       #       AgeSel_P2_FISHERY1(1)_dev_se

[k-doering-NOAA]

comment from @RickMethot on 2020-08-18:
I suspect this upper bound will be hit. I suspect that without a better
likelihood formulation employing Laplace it will not work well. But good
to confirm. I attached the relevant wtatage.ss file

Richard D. Methot Jr. Ph.D.NOAA Fisheries Senior Scientist for Stock
Assessments

Mobile: 301-787-0241

On Tue, Aug 18, 2020 at 12:28 PM vlab.redmine@noaa.gov wrote:

[k-doering-NOAA]

comment from @iantaylor-NOAA on 2020-08-18:
I just got IT to turn on my desktop at NWFSC which had somehow shut down some weeks ago but it's now available for me to remote desktop into to avoid tying up my laptop with too much processing. I just started a new set there with the upper bound restored to 14. Hopefully I'll have some results tomorrow.

[k-doering-NOAA]

comment from @iantaylor-NOAA on 2020-08-20:
I finally looked at the results of the NUTS samples for the 2 dev_se parameters. The median values were 0.066 and 0.165 and no samples above 0.4 as shown in attached figure. The posteriors.sso containing these samples is also attached along with the control file used in the sampling.

I don't know if these results match what was expected, but it was useful for me to learn more about Cole's adnuts R package and also generated some ideas for future additions to the output for deviation parameters (to be added in a separate issue).

[k-doering-NOAA]

comment from @RickMethot on 2020-08-20:
Thanks Ian. I will be curious to investigate the rmse of the parmdev
vectors

Richard D. Methot Jr. Ph.D.NOAA Fisheries Senior Scientist for Stock
Assessments

Mobile: 301-787-0241

On Thu, Aug 20, 2020 at 4:11 PM vlab.redmine@noaa.gov wrote: