bin-specific comp logL wrong in compreport for D-M; make comp logL code modular with function

[k-doering-NOAA]

Thanks Grant Thompson for reporting the following:

...I was hoping that the log likelihoods listed in CompReport.sso for a given year/fleet/kind would sum to the value listed in the FIT_LEN_COMPS or FIT_AGE_COMPS section of Report.sso, but they are not even close.

I am using the Dirichlet-multinomial option, and I note that the SS change log includes an entry for 4/5/21 referencing "change logL for D-M composition," so I was wondering if perhaps the change did not get carried over into the code that generates the CompReport.sso file (just a wild guess on my part). ...

Investigate this and see if that is what is going on. (KD: we don't test CompReport.sso at all in our testing routines, so this is certainly possible)

[Rick-Methot-NOAA]

We should also check the fit_agecomp and fit_lencomp tables in report.sso. Those tables report logL for each observation. logL in compreport.sso is only useful if the user wants a logL contribution for each bin of each comp observation.

…

On Tue, Oct 26, 2021 at 6:41 AM Kathryn Doering ***@***.***> wrote: Assigned #219 <#219> to @Rick-Methot-NOAA <https://github.com/Rick-Methot-NOAA>. — You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <#219 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABPV4IBF4QMK2R3TRCDDSJLUI2VZXANCNFSM5GX3C7SA> . 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>.

[Rick-Methot-NOAA]

I think we need a code review here.

In ss_prelim, there is:

      int z1=tails_a(f,i,1);
      int z2=tails_a(f,i,2);
      offset_a(f,i) += gammln(nsamp_a(f,i) + 1.) -
        sum(gammln(1. + nsamp_a(f,i)*obs_a(f,i)(z1,z2)));

then in ss_objfun, there is:

          age_like(f,i) = -offset_a(f,i);  //  so a perfect fit will approach 0.0

and

              if(Comp_Err_A(f)==1) dirichlet_Parm=mfexp(selparm(Comp_Err_Parm_Start+Comp_Err_A2(f)))*nsamp_a(f,i);
              if(Comp_Err_A(f)==2) dirichlet_Parm=mfexp(selparm(Comp_Err_Parm_Start+Comp_Err_A2(f)));

and

                temp+=sum(gammln(nsamp_a(f,i)*  obs_a(f,i)(tails_w(1),tails_w(2))
                               + dirichlet_Parm*exp_a(f,i)(tails_w(1),tails_w(2))));
                temp-=sum(gammln(dirichlet_Parm*exp_a(f,i)(tails_w(1),tails_w(2))));

and

              age_like(f,i)-=temp;

then in write_report, the age_like are output for each obs in fit_age table, and these sum to the total reported logL

          age_like(f,i) = -offset_a(f,i);  //  so a perfect fit will approach 0.0

then later in write_output is the compreport where bin specific logL is reported. This seems to not be doing D-M correctly.
The comp output section includes code to do DM for pearson residuals, but not for the logL. This is the problem noted by Grant.

[iantaylor-NOAA]

I'm worried that I made an error in the changes to the DM likelihood associated with change in the code to add the constant of integration as discussed under #31.
Grant confirmed by email back in August that a model run in versions before and after that change had the same parameter estimates, so the issue seems to just be with the added constant.

[Rick-Methot-NOAA]

I created a test case from "simple" in which I verified that for multinomial with both length and age comp:
a) logL reported in fit_age and fit_length matches the sum of the bin specific logL in compreport.sso
b) by creating an input observation from the expected values output, verify that the model fits this perfectly and that it's logL is near 0.0, as it should be because it includes subtraction of an offset which is the logL for a perfect fit.

Then switched to test case from hake_2018 that has only age comp and uses D-M.
It passes neither of the comparisons that are passed for multinomial.
Problem 1: calculation of offset seems not to be using the D-M parameter. Bigger issue is that the offset cannot use the D-M parameter because the D-M parameter can be estimated, hence change during the run. @iantaylor-NOAA What logic did you use when you created the current offset????
Problem 2: reporting of bin-specific logL is not using D-M calcs. But Pearson residuals do appear to be using D-M calcs
Problem 3: logL of D-M has a component that is not bin-specific. If we want the compreport's logL output to sum to match the overall logL of that observation, then this component will need to be spread across the bins of that observation.
@kellijohnson-NOAA bringing this to your awareness.

[iantaylor-NOAA]

Sorry for all these issues. I just looked through the changes in https://github.com/nmfs-stock-synthesis/stock-synthesis/pull/132/files. I thought that the whole point of the offset was that it represents the part of the likelihood which is independent of the parameters, and thus can be calculated in SS_prelim.tpl.

I tried to find notes related to checking the likelihood against the Thorson et al. paper but unfortunately haven't found them yet. Also, I'm on leave the next two weeks but can prioritize working on this for the week of Jan 3.

[Rick-Methot-NOAA]

The offset for D-M that is in place now may well be as good as it can be. Don't sweat it.

[iantaylor-NOAA]

@Rick-Methot-NOAA, I finally took the time to double-check the implementation of the Dirichlet-multinomial likelihood and reread your note above #219 (comment) from December.

I verified outside of SS3 that the DM likelihood (when choosing the linear parameterization) matches equation (10) from Thorson et al. (2017) shown below.

When revising the DM offset in Feb 2021 (#31), I failed to understand the role of the offset for multinomial likelihood as you noted in the comment referenced above. That is, for the multinomial, I now see that the goal of the offset is to make the logL approach 0.0 for a perfect fit, not to match the full likelihood equation including normalizing constant. So given the important log-likelihood term of n*obs*log(exp), the offset is set to -n*obs*log(obs).

Having the logL approach 0 for perfect fit probably has lots of benefits, including keeping the scale of all the likelihood contributions on a similar order of magnitude (in contrast to the current DM likelihood which can have really large values, and can be negative as well). However, this approach gets in the way of the one perceived benefit of having the DM likelihood include the constant term, which is facilitating comparing the DM vs multinomial likelihoods, as noted in #31 (comment).

Now that I understand the difference between the two implementations I have a few questions on possible few ways forward.

1. Would there be any purpose in giving the user the option of the multinomial likelihood with the full constant term to facilitate comparison of the two likelihood options? I presume this comparison would only work if the input sample sizes were not adjusted, so maybe this wouldn't be valuable after all since we typically always weight the comp sample sizes when using the multinomial likelihood.
2. Is there a better offset for the DM likelihood that would bring the scale of the values down, or should we just stop messing with it and leave it as is, matching the published equations?

@James-Thorson-NOAA, do you have any thoughts on this?

[iantaylor-NOAA]

Looking back over #219 (comment) above yet again, I see that the DM likelihood values should change in future models no matter what we decide on the two questions I just posed in order to resolve the problems related to bin-specific logL so that the reported values are correct and the compreport's logL output sums to match the overall logL of that observation.