Report parameters estimated at or near bounds

[Cole-Monnahan-NOAA]

Sometimes the MLE of a parameter will be at or very close to a bound. This can make SE calculations unreliable. It also depends on the bounding function used (whether -hbf was invoked).

It would be useful to give the user some sort of feedback of when this occurs. For instance after optimization it could print to the console and then add a column to the .std file which indicates an issue or not.

Defining what qualifies as near a bound is not obvious, but is something we could probably work out with some trial and error.

@johnoel do you have any thoughts on how to implement something like this?

[Cole-Monnahan-NOAA]

Maybe the simplest solution is to take the untransformed parameters and check those values. This will depend on the bounding function used (hybrid_bounded_flag value) but should in theory be easy. Loop over each parameter and if bounded and large unbounded value then print it.

This may only be a few lines of core code and so pretty easy to do right before printing final run times and such: https://github.com/Cole-Monnahan-NOAA/admb/blob/master/src/nh99/modspmin.cpp#L97

Might be worth trying real quick instead of pursuing a more involved solution.

[Cole-Monnahan-NOAA]

Here's a first pass at this on this branch. @kellijohnson-noaa and @iantaylor-NOAA and @janesullivan-noaa is this what you had in mind?

I made a toy example and it looks like this

Starting standard error calculations...  done! (0 s)

Checking for parameters on bounds (experimental)...
Par 1 appears to be on upper bound: 10
Par 2 appears to be on lower bound: 0.100001

Finished running model 'catage' after 0.24 s.

Related to https://github.com/ss3sim/ss3sim/issues/230

[Cole-Monnahan-NOAA]

OK I did another pass on this this morning. I think it's pretty close. It will only show if you use the new output_flag option (proposed default but -output 1 on command line). I split into a "near" and "on" case. Here "close" is defined on the (0,1) bounded scale as 0.001 or 0.999, and "on" is 0.00001 or 0.9999. I created a fake model to test this (pasted below) which outputs:

bound_checks -iprint 0

Starting optimization of 'bound_checks' in phase 1 of 1 at Fri Dec 10 08:51:50 2021
Optimization completed after 0 seconds with final statistics:
 nll=12.322236 | mag=6.05871e-05

Calculating Hessian: 0%, 20%, 40%, 60%, 80%, 100% done! (0 s)
Inverting Hessian: 0%, 20%, 40%, 60%, 80%, 100% done! (0 s)
Starting standard error calculations...  done! (0 s)

Checking for parameters on bounds (experimental)...
 Par 1 appears to be near lower bound: 0.000100554 (-0.987232)
 Par 2 appears to be on lower bound: 2.00202e-06 (-0.998198)
 Par 3 appears to be near upper bound: 0.999899 (0.987232)
 Par 4 appears to be on upper bound: 0.999998 (0.998198)
 Par 5 appears to be on lower bound: 2.00854e-08 (-0.99982)
 Par 6 appears to be near lower bound: 9.63058e-05 (-0.987505)
 Par 7 appears to be near lower bound: 9.76947e-05 (-0.987415)
 Par 8 appears to be near lower bound: 0.000113818 (-0.986416)
 Par 11 appears to be on upper bound: 1 (0.999819)
 Par 12 appears to be near upper bound: 0.999898 (0.987164)
 Par 13 appears to be near upper bound: 0.999902 (0.987415)
 Par 14 appears to be near upper bound: 0.999886 (0.986416)

and

bound_checks -iprint 0 -hbf

Starting optimization of 'bound_checks' in phase 1 of 1 at Fri Dec 10 08:52:41 2021
Optimization completed after 0 seconds with final statistics:
 nll=12.322198 | mag=8.39204e-05

Calculating Hessian: 0%, 20%, 40%, 60%, 80%, 100% done! (0 s)
Inverting Hessian: 0%, 20%, 40%, 60%, 80%, 100%

 Error: Estimated variance of parameter 9 is -872437, failed to invert Hessian.
        No uncertainty estimates available. Fix model structure and reoptimize

Checking for parameters on bounds (experimental)...
 Par 1 appears to be near lower bound: 0.000162438 (-8.72505)
 Par 2 appears to be on lower bound: 6.39817e-07 (-14.2621)
 Par 3 appears to be near upper bound: 0.999838 (8.72505)
 Par 4 appears to be on upper bound: 0.999999 (14.2621)
 Par 5 appears to be on lower bound: 1.24045e-08 (-18.2052)
 Par 6 appears to be near lower bound: 1.0463e-05 (-11.4677)
 Par 7 appears to be near lower bound: 1.04631e-05 (-11.4676)
 Par 8 appears to be near lower bound: 1.04748e-05 (-11.4665)
 Par 9 appears to be near lower bound: 1.17381e-05 (-11.3527)
 Par 11 appears to be on upper bound: 1 (18.2052)
 Par 12 appears to be near upper bound: 0.99999 (11.4677)
 Par 13 appears to be near upper bound: 0.99999 (11.4676)
 Par 14 appears to be near upper bound: 0.99999 (11.4665)
 Par 15 appears to be near upper bound: 0.999988 (11.3527)

Notice how there are slight differences. I think this is due to either inherent differences in bounding functions or some internal penalty. If I try the same thing with a bounded_dev_vector it breaks pretty bad (not shown).

Test model:

// model to test bound hitting

DATA_SECTION

INITIALIZATION_SECTION

PARAMETER_SECTION
  init_bounded_number par1(0,1);
  init_bounded_number par2(0,1);
  init_bounded_number par3(0,1);
  init_bounded_number par4(0,1);
  init_bounded_vector parvec(1,12,0,1);
 // init_bounded_dev_vector parvec2(1,12,0,1); 
  objective_function_value f

PROCEDURE_SECTION
  f=0;
  f+= dnorm(par1, -.01,1);
  f+= dnorm(par2, -.5, 1);
  f+= dnorm(par3, 1.01 ,1);
  f+= dnorm(par4, 1.5,1);

  f+= dnorm(parvec(1), -.5,.1);
  f+= dnorm(parvec(2), 0,.1);
  f+= dnorm(parvec(3), 0.0000001,.1);
  f+= dnorm(parvec(4), 0.00001,.1);   
  f+= dnorm(parvec(5), 0.001,.1);     
  f+= dnorm(parvec(6), 0.1,.1);   

  f+= dnorm(parvec(7), 1.5,.1);
  f+= dnorm(parvec(8), 0.999999999,.1);
  f+= dnorm(parvec(9), 0.9999999,.1);
  f+= dnorm(parvec(10),0.99999,.1);   
  f+= dnorm(parvec(11),0.999,.1);     
  f+= dnorm(parvec(12),0.9 ,.1);      

  // f+= dnorm(parvec2(1), -.5,.1);
  // f+= dnorm(parvec2(2), 0,.1);
  // f+= dnorm(parvec2(3), 0.0000001,.1);
  // f+= dnorm(parvec2(4), 0.00001,.1);   
  // f+= dnorm(parvec2(5), 0.001,.1);     
  // f+= dnorm(parvec2(6), 0.1,.1);   

  // f+= dnorm(parvec2(7), 1.5,.1);
  // f+= dnorm(parvec2(8), 0.999999999,.1);
  // f+= dnorm(parvec2(9), 0.9999999,.1);
  // f+= dnorm(parvec2(10),0.99999,.1);   
  // f+= dnorm(parvec2(11),0.999,.1);     
  // f+= dnorm(parvec2(12),0.9 ,.1);

[kellijohnson-NOAA]

What about adding output to a file for those that are not paying attention to the console output.

[Cole-Monnahan-NOAA]

@kellijohnson-NOAA I'm always hesitant to add new output files. This calculation is not always done so you could have a stale file around which would be misleading. Really I'd like to tack it on as a new column in the .std file, but that's dangerous since presumably there's a lot of downstream code that parses that and this would break that. Alternatively you could just do these calculations outside of ADMB, all you need are the bounded and unbounded parameter vectors and the HBF value. I could just write the unbounded ones to the admodel.hes file (bounded ones already there). Trivial to do in R or other software that way.

@johnoel thoughts about adding this to a new file? Really it'd be ideal to have a flag so the user could have the console output piped to a file in addition to the console. That may be a separate issue though.

[kellijohnson-NOAA]

Thanks Cole for the reminder about downstream code and not wanting to goof that up. I like the idea of perhaps providing instructions to users on how to write their own summary. For example, we could add it to the warning.sso file that SS3 produces.

[Rick-Methot-NOAA]

SS3 already produces a report showing parameters that are within 1% of bound and displays to the screen at the end of a run the number of such parameters.

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

Mobile: 301-787-0241

On Fri, Dec 10, 2021 at 10:40 AM Kelli Johnson @.***> wrote:

Thanks Cole for the reminder about downstream code and not wanting to goof that up. I like the idea of perhaps providing instructions to users on how to write their own summary. For example, we could add it to the warning.sso file that SS3 produces.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/admb-project/admb/issues/221#issuecomment-991206682, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABPV4IANQ5NGA5AWZ25FAF3UQJCR5ANCNFSM5JCQJQDQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

[k-doering-NOAA]

What about adding output to a file for those that are not paying attention to the console output.

I think from the command line you can always pipe output from a call to a file.

For example, I ran ss_3.30.18 > test.out in Powershell (shere ss_3.30.18 was the compiled admb model) and it created a file called test.out where the console output was piped to.

I think there are similar options on other command lines.

[Cole-Monnahan-NOAA]

@Rick-Methot-NOAA This is a general upgrade for ADMB so I want it to work for other models. In fact it probably makes sense for SS to use the new -output 0 option since you print a lot of stuff to screen/file anyway. You should be able to modify that in SS to have 0 as default unless the user specifies it on the command line.

@kellijohnson-NOAA The whole point of me making all these source changes is so that people do pay attention to the console. I can see in a simulation setting that you'd want to save these things to file for sure. Piping should work in those cases I suppose.

[johnoel]

An option would be to write to the log file instead of creating yet another output file.

[Cole-Monnahan-NOAA]

@johnoel Good idea. I modified the code so now the log file looks like below. I also modified so it is only run after doing the SD calculations. This prevents it from running during a mode like MCMC or likelihood profile. Or if the Hess calcs fail.

I submitted a PR for review: https://github.com/admb-project/admb/pull/222. It would be good for others to test this with different options and models.

Edit: Updated PR link. I should also say it prints to both console and log file. A few minor cleanups to output otherwise it is as above.

Checking for parameters on bounds (experimental; hbf=0)...
 Par 1 appears to be near lower bound: 0.000100554 (unbounded=-0.987232)
 Par 2 appears to be on lower bound: 2.00202e-06 (unbounded=-0.998198)
 Par 3 appears to be near upper bound: 0.999899 (unbounded=0.987232)
 Par 4 appears to be on upper bound: 0.999998 (unbounded=0.998198)
 Par 5 appears to be on lower bound: 2.00854e-08 (unbounded=-0.99982)
 Par 6 appears to be near lower bound: 9.63058e-05 (unbounded=-0.987505)
 Par 7 appears to be near lower bound: 9.76947e-05 (unbounded=-0.987415)
 Par 8 appears to be near lower bound: 0.000113818 (unbounded=-0.986416)
 Par 11 appears to be on upper bound: 1 (unbounded=0.999819)
 Par 12 appears to be near upper bound: 0.999898 (unbounded=0.987164)
 Par 13 appears to be near upper bound: 0.999902 (unbounded=0.987415)
 Par 14 appears to be near upper bound: 0.999886 (unbounded=0.986416)

size of file gradfil1.tmp = 0
size of file gradfil2.tmp = 0
size of file varssave.tmp = 0
size of file cmpdiff.tmp = 0

[johnoel]

Merge in pull #235