inclusion of priors in likelihood?

[iantaylor-NOAA]

In the Big Skate STAR we were asked to produce catchability profiles both with and without the influential prior (shown below).

• I like the with-prior version (right) because the minimum of the total likelihood curve matches the MLE and helps see how influential the prior is vs the data.
• The review panel liked the without-prior version because you can see more detail and see where the data would have pointed that parameter had we not included a prior.
• Both versions can be useful although it would be convenient not to include two figures for every profile over a parameter with a prior.

Going back and forth can be done without re-running the profile if you initially include the prior and then optionally subtract the prior likelihood for that parameter from both the total and the prior likelihood (which may include contributions from other priors).

@chantelwetzel-noaa, sorry to keep dragging you into lingcod-related stuff, but can you help me understand the nwfscDiag settings related to priors. Looking at https://github.com/nwfsc-assess/nwfscDiag/search?q=prior_like, it appears that the default is to turn off the prior contribution to the likelihood for fixed parameters, which would get you the without-prior version of the profile, but then how does it add the extra model for the MLE to the curve without including the prior? Does it include the point estimate of the MLE for that parameter as a fixed value within the vector you provide? Will changing to model_settings$prior_like <- 1 include the prior if I want it?

[Edit: checklist added on 8 June 2023]

• remove use_prior_like option (and associated code related to prior_used variable in the profile wrappers)
• do some testing

[chantelwetzel-noaa]

If you would like to include the prior likelihood contribution, yes you would add model_settings$prior_like <- 1 to your model_settings object. I think there are different instances when we may want to include or exclude the prior likelihood. Much of my experience comes from thinking about steepness or natural mortality where we often what to understand what the data are saying about a parameter without the influence of the prior. In the 2017 SSC stock assessment review for POP the guidance was them was very clear that the prior likelihood should not be included in at least the steepness profiles.

However, there are instances when one may be estimating a parameter with a prior where not including the prior likelihood contribution in the profile can result in the MLE estimate not falling along the profile like across the total likelihood or the derived quantities.

After the assessment dust settles we should have a discussion as a group as to when we should include or exclude the prior likelihood. Additionally, we could add additional functionality into the package where we could run profiles both ways when a non-uniform prior is present in a single call to the run_diagnostics function.

[iantaylor-NOAA]

Thank you @chantelwetzel-noaa. Another idea would be to have separate lines for total with prior and total without prior. Regardless, the nwfscDiag package is working really well for lingcod and our reviewers will likely be similar to your POP reviewers and happy with what we've produced already.

[chantelwetzel-noaa]

@kellijohnson-NOAA @iantaylor-NOAA I have completed the addition of the ability to specify including or excluding the prior likelihood in a single call in commit bcf3d20. The changes that users need to be aware of is the added input to the get_settings_profile function (use_prior_like). Here is an example call which would result in two profiles over female M, one without the prior likelihood contribution and one with:

get = get_settings_profile( parameters = c("NatM_p_1_Fem_GP_1", "NatM_p_1_Fem_GP_1"),
low = c(0.08, 0.08),
high = c(0.14, 0.14),
step_size = c(0.01, 0.01),
param_space = c('real', 'real'),
use_prior_like = c(0, 1))

I am going to leave this issue open pending testing by the lingcod team.

[shcaba]

Do we have guidance on whether/when this should be on or off,? I assume consistency with the reference model is a base consideration, but I know there are thoughts pro and con to having it on.

[chantelwetzel-noaa]

@iantaylor-NOAA and I discussed this a bit above. I think whether you want the prior likelihood contribution included depends upon the situation. I typically run without including the prior likelihood because I want to see what the data support. However, excluding the prior likelihood can result in profile results that are not entirely consistent with the model estimation if the parameter is estimated with a prior.

[shcaba]

Thank you. That is where I thought we had left it; I wanted to make sure I hadn't missed any details.

[kellijohnson-NOAA]

@chantelwetzel-noaa the naming scheme of having the folders named differently with priors or without is going to cause some issues with our documents, i.e., we wouldn't be able to share code across assessments because they would be referencing different folders. Thoughts?

[chantelwetzel-noaa]

I was focused on allowing users the flexibility to run with and without the prior for the same parameter within a single call to the function. That required a distinct folder naming convention in order to avoid erasing the previous run and to also allow people to distinguish the runs from each other. You are right that this will impact existing documents that are referencing profile runs. I failed to factor that in and I apologize. I am unsure of an alternative approach at this moment. I am open to ideas.

[kellijohnson-NOAA]

I think your approach of being explicit with the name is best. There are many ways we can deal with searching for a name on the downstream end. I will brainstorm.

[chantelwetzel-noaa]

The new approach has been used and tested by the lingcod assessment team. Going to close this issue but can be opened if issues are encountered

[kellijohnson-NOAA]

I am reopening this issue to remind us that we need to document the "best practices" part of this in the argument documentation for use_prior_likelihood because I think that we can be more specific than "in certain instances".
https://github.com/nwfsc-assess/nwfscDiag/blob/2213b342dde0e680ea9c325e053d657dc9ff97e8/R/get_settings_profile.R#L24-L28

[iantaylor-NOAA]

Comments during team meeting discussion today:

@iantaylor-NOAA suggested that we create a new figure that has with the prior and without the prior. 🏆 for everyone!

If you fix the value you have a "point prior" because you are fixing it rather than the actual prior.

3 situations which we want to work:

• parameter estimated with a prior (often M)
• parameter which is fixed (do we want to present the point prior?)
• parameter which is estimated with no prior (but influenced by priors on other parameters), like log(R0)

[iantaylor-NOAA]

I'm working on the "🏆 for everyone" solution to this issue in r4ss/r4ss#750.
@chantelwetzel-noaa, if there are changes to the implementation that would make it easier to incorporate into {nwfscDiags} (or just make you happier), please let me know.

[iantaylor-NOAA]

Low priority future change:

This issue got closed automatically when I merged r4ss/r4ss#822. However, to fully resolve this, I think that {nwfscDiag} should be simplified to remove the option to exclude the prior likelihood now that the figure shows the total with and without any prior for the parameter being profiled. I added that change as a checklist item on the top comment and am assigning @chantelwetzel-noaa, but obviously this is not a top priority as the function will still work as it stands. If the user sets use_prior_like = 1 then they should see the total with and without prior in the plot. If they set use_prior_like = 0 they will only see the without prior total.

@chantelwetzel-noaa, if you're happy with that option remaining, please close this issue and don't think about it again.

[chantelwetzel-noaa]

I agree @iantaylor-NOAA. Thank you for re-opening this issue. I plan to work on this issue next week.

[iantaylor-NOAA]

I just discovered a logical flaw in my thinking around this issue while writing up the petrale report.
I thought that the dashed line in the plot below showed that the M prior had relatively little influence on the model results, which was surprising. However, female and male M are highly correlated (0.96) and the profile on female M leaves the prior in place on male M, so effectively there's still a prior.

Three solution that I'm not going to pursue today are:

• turn off the male M prior for the female M profile
• profile in 2 dimensions over male and female M
• set male M as an offset to female M instead of an independent parameter