[Feature]: implement survey errtype that adds lognormal bias adjustment

[Rick-Methot-NOAA]

Describe the bug

Expected to see different survey Q when bias_adj flag is turned on

To Reproduce

Examine the survey INDEX01 output for a survey that does and a survey that doe not have different se among the obs. Do this with four conditions: bias_adj on vs. off; float on vs. off. use this example: simple_small_biasadjQ_varySE.zip

Expected behavior

different Q

ToDo

• [ ] create new options for svy_errtype: -2=gamma (Cole); -1=normal; 0=lognormal; >1XXX=T. Add 1 = lognormal with bias adj.. Note that this means that the T-dist cases must be changed to use the revised code.
• [ ] Explain to user the use of bias adjustment inherent to error_type as input in data file, vs. explicit bias_adjustment flag in control file.

[Rick-Methot-NOAA]

@iantaylor-NOAA @Cole-Monnahan-NOAA @e-gugliotti-NOAA Let's use this info to augment the Q section of the user manual. fleet 1 has varying se; fleet 2 has same se each year

err_type	bias_adj	float	fleet1_ln(q)	fleet2_ln(q)	fleet1_logL	fleet2_logL
0	0	0	0.670859	-6.31285	-1.17949	-3.33496
0	0	1	0.670859	-6.31285	-1.17949	-3.33496
0	1	0	0.670859	-6.31285	-1.17949	-3.33496
0	1	1	0.797298	-6.06785	-1.05305	-2.59996
1	0	0	0.791915	-6.07255	-1.38859	-3.33224

only one scenario shows effect of bias adjustment. Note that the bias_adjustment to Q is only applied when doing the "float" approach. In the code, I see the following note:

            if (Svy_use(f, i) > 0)
            {
              surv_like(f) += 0.5 * square((Svy_obs_log(f, i) - Svy_est(f, i)) / Svy_se_use(f, i)) + sd_offset * log(Svy_se_use(f, i));
              //            should add a term for 0.5*s^2 for bias adjustment so that parameter approach will be same as the  biasadjusted scaling approach
            }

[e-perl-NOAA]

@Rick-Methot-NOAA Would it be worth including a warning when people turn on bias adjust but don't also turn on float, to go back and do that?

[Rick-Methot-NOAA]

Yes! That makes sense. Will do.

[Rick-Methot-NOAA]

Larry's work in 2012 also added lognormal with bias adjustment with intention of using it in the prior on Q: https://vscode.dev/github/nmfs-ost/ss3-source-code/blob/main/SS_objfunc.tpl#L1587-L1597

But I later wrote to him (Oct 2012): "Any help you can provide on this would be much appreciated. I still have the spreadsheet in which we compared the two methods for internal Q scaling to the parameter method; which disclosed that the parameter method needed a bias correction term. Regarding below: since the Q parameter has units of log(Q), why not put a normal prior on it rather than a lognormal prior?"

It now seems the the bias adjustment approach we worked out in 2012, and now raised again by Cole, was never implemented. I will create a new survey err_type option so we can test this out. Attaching the spreadsheet Larry and I were using in 2012.

from Cole:

[Rick-Methot-NOAA]

Qmethod2.xlsx

[Rick-Methot-NOAA]

Here is my current state of some progress:

Two surveys, one with constant observation se and the other with time-varying se Four scenarios: Float; float with bias corr; parm; parm with bias corr Note that bias corr with float is invoked in the control file; but bias adj with params is invoked with err_type in the data file!!! Who designed this?

Results:

1. Parm and Float produce identical q, expected values, logL and rmse. As expected. Plots show that the devs are centered on zero.

2. Float,bias shows changed q, degraded rmse and degraded logL. Plot shows that the devs are shifted towards negative because the q and expval were increased by the bias corr.

3. parm, bias is interesting. It is only scenario in which the model changes other parameters such that the vuln_bio is different. With this adjusted population, this scenario produces a similar q to the float, bias option. It also produces a similar rmse to float, bias with shift in the devs towards negative values. Even though the rmse is worse than without bias corr, the logL is surprisingly better. The plots do not seem to give a clue as to why.

Inside the code: expval code calculates Q if Q is based on parameters, calculates the survey abundance without Q, and if Q is based on parameters it applies that Q to create expected values. This section of code does not differentiate between bias corrected or not.

objfun code starts by calculating float q if that option is taken. If bias corr is selected, then the q calc is weighted by the se values. The svy_est values are based on the bias corr, or not, q values.

Then the logL is calculated. Here there is a lognormal and a bias corrected lognormal option. So, it float w/bias is used, then it would be double correcting if lognormal w/bias was selected. This should warrant a warning.

Qbias_test.xlsx

[Cole-Monnahan-NOAA]

Thanks Rick I'll take a look. While the expected indices would match I'd expect the float option to produce different uncertainties. Is that true?

I'm also not convinced you can compare those loglik values without including all the constants. The LN case will have a different constant. According to Burnham and Anderson we can compare AIC (thus likelihoods) under the conditions that they have the same units (discrete vs continuous) and the constants are included. I'm still struggling to work through a math/stats reason for not including the bias correction and what that implies about the distribution of the observed index. I'll keep thinking on that but open to other perspectives or help working through it. I couldn't find a good justification in the literature.

[Rick-Methot-NOAA]

No difference in the se on any derived quantities. Expect for the expected value of the survey itself, which I can rationalize but may not have predicted.

[Cole-Monnahan-NOAA]

I'm surprised by that. I would think that there'd be less uncertainty in expected index using the float option, and that would propagate through to the rest of the model. How can you be less certain about the fit to one data set and have it not impact the rest? Tough to think about.

[Rick-Methot-NOAA]

Same is true for F(yr). These can be either internal tuning coefficients that SS3 constantly adjusts to come close to matching the input catch values, or can be model parameters. Even have option to start as coefficients and change to parameters during a specified between_phases. Either way, the derived quantities (like end year SSB) have the same variance.

[Cole-Monnahan-NOAA]

If you output a sdreport for q_float does that have the same SE as if it's an estimated parameter? It's a bit off topic but I'm still wrapping my head around that.

I wrote a short script to demonstrate the constant difference I mentioned above. You can assume that some r.v. X is lognormally distributed, X~LN(\mu,\sigma). Or you can assume that log(X)~N(\mu,\sigma). These are different. The latter transforms the state variable and thus the likelihoods are not strictly comparable. In this case they differ by a constant from a likelihood perspective, but I think it's important to think of them as different. You can see the impact of the bias correction so it's no surprise that you get a different fit.

## Play around with likelihood functions
mhat <- 3                               # DB estimate of biomass
shat <- 1                               # SE(mhat)
CV <- shat/mhat
shatlog <- sqrt(log(CV^2+1))     # SE in log space
Iseq <- seq(1,5, len=1000)              # expected index from assessment
ll1 <- data.frame(I=Iseq, nll=-dnorm(mhat, Iseq, shat, log=TRUE), v='normal')
ll2 <- data.frame(I=Iseq, nll=-dnorm(log(mhat), log(Iseq), shatlog, log=TRUE), v='normal log')
ll3 <- data.frame(I=Iseq, nll=-dlnorm(mhat, log(Iseq), shatlog, log=TRUE), v='lognormal')
ll4 <- data.frame(I=Iseq, nll=-dlnorm(mhat, log(Iseq)-shatlog^2/2, shatlog, log=TRUE), v='lognormal bias')
out <- rbind(ll1,ll2,ll3,ll4) %>% group_by(v) %>%
  mutate(nlldiff=nll-min(nll), mode=I[which.min(nll)]) %>% ungroup %>%
  pivot_longer(cols=c('nll', 'nlldiff'))
ggplot(out, aes(I, y=value, color=v)) + geom_line() + facet_wrap('name')
filter(out, I==Iseq[1])

Results in this plot: