SeptDoc.qmd

---
title: "Eastern Bering Sea pollock stock assessment model evaluations"
author:
  - name: James Ianelli
    id: ji
    orcid: 0000-0002-7170-8677
    email: jim.ianelli@noaa.gov
    affiliation: 
      - name: Alaska Fisheries Science Center
        city: Seattle
        state: WA
        url: https://www.fisheries.noaa.gov/alaska/population-assessments/north-pacific-groundfish-stock-assessments-and-fishery-evaluation
csl: cjfas.csl
date: now
date-format: "DD MMM YYYY"

format: 
  pdf:
    number-sections: true
    include-in-header: 
      text: |
        \usepackage{scrlayer-scrpage}
        \rohead{EBS Pollock assessment considerations}
        \lohead{Draft}
        \lofoot{September 2023 }
        \rofoot{BSAI Plan Team}
    includes:
      in_header: "doc/preamble.tex"
  
bibliography: doc/references.bib
---

```{r setup}
#| echo: false
#| warning: false

#source(here::here("R","prelims.R"))
library(ebswp)
load(file="doc/septmod.rdata")
load(file="doc/modtune.rdata")
thisyr=2022
nextyr=thisyr+1
#names(modlst) <-  c("Last year", "Diagonal cov BTS","GenGam","SSB0","SSB1","SSB2","AVO new","AVO full")
mod_names<-names(modlst)
mod_tune_names<-names(modtune)
knitr::opts_chunk$set(fig.pos = "H", out.extra = "")
#\label{tab:modcomp} 
#        #\usepackage{fancyhdr}
        #\pagestyle{fancy}
        #\fancyhead[C]{Draft}
        #\fancyfoot[C]{EBS pollock assessment}
```

# Background

Each year we attempt to show how the Eastern Bering Sea walleye pollock assessment model is affected by data and assumptions. In most recent years, the evaluations have arisen from alternative methods for data processing (e.g., use of the VAST index, including NBS region, etc.). This year we evaluate new treatment of the acoustic data collected from the bottom-trawl survey data (the acoustic vessels of opportunity (AVO) series). In addition, we updated and re-evaluated some of the other input data specifications. This included reconsidering the relative weights among different data sets. We present alternative process-error specifications (e.g., the amount of allowed variability among processes related to time-varying selectivity) and compare that with an alternative where additional observation error variance was estimated. We also began to take the steps needed to be able to compare the current model across different assessment software platforms. Table 1 below shows the model evaluations examined, the naming convention used, relevant comparisons, and notes.

These models are contrasted and categorized by the type of analysis involved. Specifically, we start with the 2022 assessment model and data configuration, then apply the revised new AVO series (AVON; @fig-schema). Here we assumed that the application of the more comprehensive (spatially) time series would be accepted by the SSC and Plan Team. The next step in the model process was to re-evaluate the relative weights specified as input variances and sample sizes. We attempted to apply different tuning approaches which achieved a balance between the observation errors and model process errors. We then have a section where we examined some new model options followed by sensitivities to spawning biomass weight-at-age assumptions.

The purpose of this document is to have a discussion on some of these options and direction for the November assessment and beyond; ideally, advice for bridging to a set of candidate base models for this year's assessment would be ideal.

For the 2023 assessment, we started responses to the SSC requests and tracking plans; these are available [here](https://afsc-assessments.github.io/ebswp/articles/Assessment-planning.html "Link to assessment planning").

<!-- \begin{tabular}{p{1.5in}p{.5in}p{.75in}p{4.5in}} -->

```{=tex}
\hypertarget{tbl-plan}{}
\begin{table}[ht]
\centering
\caption{\label{tbl-plan}Preliminary model configuration comparisons for the EBS pollock stock assessment, September 2023.}
\tabularnewline
\scalebox{0.80}{
\begin{tabular}{p{1.5in}p{5.7in}} 
\hline
Description                    & Notes                                                                                                                                                                                                                                       \\
\hline
Base                           & Original from folder 07, now in folder base22 with new code that accommodates arbitrary number of age-error matrices by gear and year, and the generalized Gamma distribution for index data                                                                                                                                            \\
\hline 
AVO new series                 & Newly integrated acoustic backscatter from "Acoustic vessels of opportunity"--namely the bottom-trawl survey vessels. Includes data from 2009 (data from early and in 2011 and 2013 omitted). \\
AVO full                       & As in "AVO new" series but with data from 2006-2008, 2011, and 2013 used from previous series (rescaled to have the same mean)\\
\hline 
Tuned TV selectivity for ATS   & As in "AVO full" but with time-varying (TV) acoustic-trawl survey (ATS) selectivity variability "tuned"  to achieve model consistency with input variance terms\\
Tuned observation errors for 
indices                        & As in "AVO full" but with input index variability estimated for consistencey with model fit\\
\hline 
Ageing error                   & Include age-determination error matrix \\
Diagonal                       & As in base but uses only the diagonal of the covariance matrix for the bottom-trawl survey (BTS) biomass estimates                                                                                                                                                         \\
Generalized Gamma distribution applied to BTS data & This is based on the MCMC posterior samples generated from the process of accounting for density-dependence within tows; the goal is to evaluate annually varying distribution assumptions compared to the standard lognormal distribution. \\
\hline 
Alternative wt-age             & Show how stock changes  assuming a constant weight-at-age for SSB calculations \\
for SSB                        & Use empirical mean values from A-season fishery as a proxy   \\
                               & Use predicted mean values from the year and cohort random-effects model        \\
\hline 
Base Proposed                  & Full AVO series and with the smoothed mean-wt-age (random effect model)      \\
\hline 
\end{tabular}
}
\end{table}
```
<!-- AMAK platform                  & amak    & diag2         & As with "diag2" data but uses "AMAK2" software configured to be similar to base22 \\ -->

![Schematic of approach for evaluating some planned changes in the EBS pollock assessment. See text for details.](doc/figs/Schema.png){#fig-schema}

{{< pagebreak >}}

# Methods

## Revised opportunistic acoustic data series

The new time series of data for the AVO data shows a peak estimate in 2015 and less of an increase in 2019-2022 (@fig-avodata). These new data represent a broader area and are considered an improvement over the previously calibrated grid cells (N. Lauthenberger, AFSC pers. comm.). The scales of each of these series is set to have the same mean over the updated years of available data. The results from the newly analyzed back-acoustic scatter records are shown as "X" in the figure; the in-between and earlier years were set to the previous year's values (normalized so the means are the same).

![Time series of EBS pollock data from the acoustic vessels of opportunity (AVO) showing the years of new data compared to previous series and the full series.](doc/figs/AVO_Series_change.png){#fig-avodata}

## Model tuning: observation versus process-error specifications

To re-visit data weighting among indices and age compositions we attempted to develop a simplistic approach. We tried to show model consistency by either inflating (or reducing) as needed the input specified variances (model correctly specified but data conflicts suggests observation errors too low, @Francis2011). Alternatively, we assumed process errors were mis-specified and that the input data variance terms were "correct." The goal was to contrast different statistical properties common in modern stock assessments (i.e., estimating additional variance versus adjusting model flexibility (@stewart2017fra, @Martell2013a). Here we defined a model as being "consistent" when the specified (or estimated) observation-error variances were about equal to the residual variances from the model fit.

To achieve consistency using the observation-error estimation approach, we iterated on the magnitude of the mean coefficient of variation (CV) for each index. For example, if an index had an SDNR of 1.2, indicating that the model fit was inconsistent with the specified variance terms, then the annual CVs for that index would be all multiplied by 1.2. This approach assumes the process/model specification is appropriate and that the observation errors require tuning. We note that in other settings such "added variance" terms can optionally be freely estimated (e.g., @methot2013) but here it took only four iterations to achieve standard deviations of the normalized residuals equal to 1.0 (@Francis2011).

For the process error term we focused on the index that appeared to be mis-specified (i.e., the residual variance was greater than the specified input variances). This meant re-evaluating the degree of process error terms for the acoustic-trawl survey selectivity (in the base model, the selectivity was assumed to be constant). Results from initial trials suggested that adding flexibility to this process error term could achieve estimation consistency (where the variances specified are consistent with the model's fit to the data).

## Weight-at-age used for spawning biomass

In the base model used for 2022, the spawning biomass body weight-at-age used was set to be equal to the estimated annual fishery values. The rationale for using fishery data was partly based on the fact that sample sizes from the fishery are much larger than those available from any survey data and the fact that the fishery covers the season when peak spawning occurs (in early spring / late winter). Additionally, the methods used to estimate near-term weight-at-age considers the recent mean values, the estimation uncertainty of those means (i.e., inverse-variance weighted), and accounts for cohort-effects (e.g., see @ianelli2022). New developments (e.g., @cheng2023) include refinements in the ability to estimate mean weights-at-age with appropriate measures of uncertainty within assessment models. Consequently, we re-evaluated a how fishery data could best be applied as a measure of reproductive output (fish body weight-at-age). Specifically, we evaluated available data drawn from the spawning season.

We apply the original method developed in Ianelli et al. (2007). This approach uses actual observed fishery mean and variances in weight-at-age over the period 1991-`r thisyr-1` from the A-season (January-March). To supplement these data, we also include current-year survey observations (see @indivero2023). The model for this is: \begin{align}
\hat{w}_{ta} &= \bar w_a e^{\upsilon_t} & a=1, \, t \ge 1964 \\
\hat{w}_{ta} &= \hat{w}_{t-1,a-1} + \Delta_a e^{\psi_t} & a > 1, \, t > 1964 \\
\Delta_a     &= \bar w_{a+1} - \bar w_a & a<A  \\
\bar w_a     &= \alpha \left\{L_1+ \left(L_2-L_1\right)\left(\frac{1-K^{a-1}}{1-K^{A-1}}\right)\right\}^3  \\
\end{align} where the fixed effects parameters are $L_1, L_2, K,$ and $\alpha$ (length at age 1 and 10, growth rate, and scalar, respectively) while the random effects parameters are $\upsilon_t$ and $\psi_t$ for year and cohort effects.

The random effects model used to estimate weight-at-age for the spawning biomass was based on the A-season fishery and the summer survey data (scaled to have the same mean as the A-season means) which shows distinct year and cohort effects (@fig-ssbwt). @fig-wt_ssb shows how the patterns differ.

![Weight-at-age data (left and middle panel) and predictions (right-most panel) for EBS pollock as an alternative series for spawning biomass weight-at-age. Note that the color shadings indicate the anomaly within age (columns).](doc/figs/ssb_wtage_data_pred.pdf){#fig-ssbwt}

![Alternative scenarios for body weight-at-age used for computing spawning biomass for EBS pollock biomass weight-at-age.](doc/figs/wt_ssb.pdf){#fig-wt_ssb}

For assessment model evaluations we compared model runs where the mean body weight-at-age was constant over time with the empirical data from the A-season and additionally runs with the smoothed random-effect model result (the $\hat{w}_{ta}$ ). We consider the random effects approach to provide a simple generalized approach from which to use the available data and account for uncertainty across years and ages.

## Other model evaluations

Some model alternatives were examined to either explore new approaches or features that exist but are turned off. For example, the current pollock base model has age-determination error as an option but it is turned off partly due to earlier considerations from the work of @kimura1992. With the advent of the potential application of Fourier Transform near-infrared spectroscopy (FT-NIRS) for age determination (@helser2019; @healy2021), the model was modified such that age compositions could be fit by any number of error matrices over time and among gear types. For now, we used the previously estimated age-determination error to re-evaluate the effect on model fit and results (the FT-NIRS data are still in development).

Another feature of the current model is the use of the estimated covariance matrix over years from a combination of density-dependent corrections (see @kotwicki2014) and VAST spatial distribution modeling (@thorson2018b) for the BTS series. While this has been a feature of the model since 2015 and part of the current base-case model, it is rare to find other assessment software platforms that have this capability. Consequently, in an effort to compare the current model with other software platforms we tested a model where the variances of the individual years in the bottom-trawl survey (BTS) are treated independently (i.e., using just the diagonal of the covariance matrix).

Finally, Monnahan (pers. comm.) has developed an approach for using a generalized gamma distribution which allows for different scale and shape parameters for each year of the bottom trawl survey. We provide that for exploration purposes since there is evidence that the statistical distribution of the CPUE from survey trawls can vary from the assumed log-normal error distribution.

Across model alternatives, we attempted to provide consistent tables for judging goodness of fit, impact on status measures (e.g., current biomass and reference points) and figures showing impact on trends in spawning biomass.

# Results

## Alternative AVO data treatment

The relative fits to the AVO data degraded with the new short and full series (@tbl-modfitsavo). However, there was relatively minor impact on biological reference points and stock status (@tbl-modrefavo).

```{r modfitsavo}
#| label: tbl-modfitsavo
#| tbl-cap: "Goodness-of-fit measures to primary data for different assessment model configurations. RMSE=root-mean square log errors, NLL=negative log-likelihood (may not be comparable across model configurations), SDNR=standard deviation of normalized residuals, Eff. N=effective sample size for composition data)"
#| echo: false
#names(modlst[5:6]) = 
df <- tab_fit(modlst[c(2:4)])
df |> gt::gt() |>  gt::fmt_number( rows = c(8:21),decimals=0)
```

```{r tbl-modrefavo}
#| label: tbl-modrefavo
#| tbl-cap:  "Summary of model results and the stock condition for EBS pollock. Biomass units are thousands of t."
#| echo: false
df <- tab_ref(modlst[c(2:4)])
df |> gt::gt() |>  gt::fmt_markdown()
```

## Model tuning

Assuming the AVO revised data will be accepted, we began "tuning" the model with this series adopted (instead of working from the 2022 base model and data). As noted above, this entailed examining trade-offs between observation and process error assumptions. We modified the ability for one index (ATS) to have time-varying selectivity-at-age and iterated on the process-error variance term to satisfy the specified observation errors via conforming to standard deviations of normalized residuals being 1.0 or close (@Francis2011).

The alternative tuning approaches affected the relative fits and uncertainty assumptions. The observation errors for the AVO data increased in both model sets of tuning approaches whereas the error bars for the process-error approach increased only slightly for the ATS index @fig-mod_index_fits_tuned). For the BTS data, the fits and impact on error-terms were very similar (presumably due to the covariance matrix used for that data component (@fig-mod_bts_fit_tuned). The process-error tuned model resulted in ATS selectivity estimates that varied slightly over time (@fig-sel_ats_tuned).

<!-- Updating the ATS data resulted in a different fit to the data trend of those data due to the lower overall uncertainty of the observations specified in revision (@fig-atsalt).  -->

<!-- This change had a minor impact on the recent estimates of spawning biomass (@fig-atsaltssb). -->

![Comparison of the model fit to acoustic trawl survey (top) and acoustic vessels of opportunity (AVO; bottom) for the 2022 assessment and tuning process error terms or observation errors.](doc/figs/mod_index_fits_tuned.pdf){#fig-mod_index_fits_tuned}

![Fit of EBS pollock model to bottom trawl survey data from the 2022 error-term specification.](doc/figs/mod_bts_fit_tuned.pdf){#fig-mod_bts_fit_tuned}

{{< pagebreak >}}

![Time varying acoustic trawl selectivity estimates for the process-error tuned model.](doc/figs/sel_ats_tuned.pdf){#fig-sel_ats_tuned}

{{< pagebreak >}}

```{r modfittune}
#| label: tbl-modfitstune
#| tbl-cap: "Goodness-of-fit measures to primary data for different assessment model configurations. RMSE=root-mean square log errors, NLL=negative log-likelihood (may not be comparable across model configurations), SDNR=standard deviation of normalized residuals, Eff. N=effective sample size for composition data)"
#| echo: false
#names(modlst[5:6]) = 
 tab_fit(modtune[c(1:2)]) |> gt::gt() |>  gt::fmt_number( rows = c(8:21),decimals=0) 
# tab_style( style = cell_borders( sides = c("bottom"), color = "black", weight = px(1.5), style = "solid" ),
#     locations = cells_body( columns = everything(),
#       rows = c(4,7,10) ) ) 
```

{{< pagebreak >}}

```{r tbl-modreftune}
#| label: tbl-modreftune
#| tbl-cap:  "Summary of model results and the stock condition for EBS pollock. Biomass units are thousands of t."
#| echo: false
tab_ref(modtune[c(1:2)]) |> gt::gt() |> gt::fmt_markdown()
```

## Spawning biomass weight-at-age

Compared to the values used in the base and "tuned" model, the trends are very similar but the scales change substantially in some years (@fig-ssbssb ; including the "SSB0" scenario). As expected, changing the spawning biomass weight-at-age used to calculate spawning biomass had no effect on the model fit (@tbl-modfitsssb). Also, the stock status and reference points varied among assumptions, but were relatively similar (@tbl-modrefssb).

![Results of the EBS pollock model for recent spawning biomass estimates comparing the base model under different assumptions about mean body weight-at-age during spawning.](doc/figs/mod_ssb_ssb.pdf){#fig-ssbssb}

{{< pagebreak >}}

```{r modfitsssb}
#| label: tbl-modfitsssb
#| tbl-cap: "Goodness-of-fit measures to primary data for different assessment model configurations. RMSE=root-mean square log errors, NLL=negative log-likelihood (may not be comparable across model configurations), SDNR=standard deviation of normalized residuals, Eff. N=effective sample size for composition data)"
#| echo: false
#names(modlst[5:6]) = 
tab_fit(modlst[c(5,9:11)]) |> gt::gt() |>  gt::fmt_number( rows = c(8:21),decimals=0)
```

{{< pagebreak >}}

```{r tbl-modrefssb}
#| label: tbl-modrefssb
#| tbl-cap:  "Summary of model results and the stock condition for EBS pollock. Biomass units are thousands of t."
#| echo: false
tab_ref(modlst[c(5,9:11)]) |> gt::gt() |>  gt::fmt_markdown()
```

{{< pagebreak >}}

## Alternative error-term specifications

In this section we present the updated error terms associated with the acoustic trawl data (the observations used previously were revised). 
We also show the impact of using the diagonal covariance matrix and an alternative for the bottom trawl survey data based on the generalized 
gamma distribution in place of the log-normal. Changing the bottom-trawl survey error to be based on independent annual survey index 
values (the base model used a covariance matrix estimated for the series) was nearly identical compared to the base tuned model (@fig-diagssb). 
However, the fit to the BTS index was worse than for the other alternatives based on the root-mean squared error (RMSE; @tbl-modfitserr; note 
that the for the "NLL" for the  for the BTS differ because the form of the likelihood equation changed). This degradation of fit to the BTS index
is also apparent in @fig-mod_bts_diag. 
The results on goodness of fits for the Generalized Gamma distribution was improved compared to the base tuned model and the others 
and the pattern of fits of the model predictions varies compared to the base tuned model (@fig-mod_bts_gengam). Presumably this is due to the
externally estimated shape and scale of the annually varying distribution (while for the base model the "shape" is constant
by assuming a multivariate log-normal distribution).

For the runs with the age-determination error option 
turned on, the fits to the data degraded for most of the age composition data (based 
on the NLL and the "effective-sample size") but improved fits 
slightly to the indices  (e.g., in @fig-mod_bts_ageerr, @tbl-modfitserr). 

For the models with modifications to the error-term specifications, the changes for biomass
reference points varied slightly (@tbl-modreferr).

![Results of the EBS pollock model for recent spawning biomass estimates comparing the base model using the covariance matrix with the one where only the diagonal is applied.](doc/figs/mod_diag_ssb.pdf){#fig-diagssb}

![Results of the EBS pollock model fit to the BTS data comparing a model using the covariance matrix (tuned22) with the one where only the diagonal is applied ("Diag cov BTS").](doc/figs/mod_bts_diag.pdf){#fig-mod_bts_diag}

![Results of the EBS pollock model fit to the BTS data comparing a model using the standard log-normal distribution is used (tuned22) with the one where the annually-varying generalize gamma distribution is applied ("GenGam").](doc/figs/mod_bts_gengam.pdf){#fig-mod_bts_gengam}

![Results of the EBS pollock model fit to the BTS data comparing a model using the standard log-normal distribution is used (tuned22) with the model where the age-determination error option is turned on    ("Ageing Error").](doc/figs/mod_bts_ageerr.pdf){#fig-mod_bts_ageerr}

{{< pagebreak >}}

```{r modfitserr}
#| label: tbl-modfitserr
#| tbl-cap: "Goodness-of-fit measures to primary data for different assessment model configurations. RMSE=root-mean square log errors, NLL=negative log-likelihood (may not be comparable across model configurations), SDNR=standard deviation of normalized residuals, Eff. N=effective sample size for composition data)"
#| echo: false
tab_fit(modlst[c(5:8)]) |> gt::gt() |>  gt::fmt_number( rows = c(8:21),decimals=0) 
# |>  tab_style( style = cell_borders( sides = c("bottom"), color = "black", weight = px(1.5), style = "solid" ),
    # locations = cells_body( columns = everything(),
      # rows = c(4,7,10) ) ) 
```

{{< pagebreak >}}

```{r tbl-modreferr}
#| label: tbl-modreferr
#| tbl-cap:  "Summary of model results and the stock condition for EBS pollock. Biomass units are thousands of t."
#| echo: false
tab_ref(modlst[c(5:8)]) |> gt::gt() |>  gt::fmt_markdown()
```

# Summary

For the models we explored there was generally improvements to the model specifications that could be made. These were relatively minor impacts on the stock status and trends. For discussions, we therefore recommend:

-   Adopting the use of the full revised AVO time series. *The data have been re-calibrated to the acoustic trawl survey and cover a larger area than previous AVO series. This may improve the ability to track expansion and contraction of the pollock stock in mid-water.*

-   Allowing modest process error terms consistent with the observation-error specifications. *We prefer to follow allow slight deviations in process errors (here selectivity/availability in the ATS) to achieve consistency with specified survey-based observation-errors*

-   Adopting the use of the RE model for A-season fishery mean body weight-at-age. *These smoothed values are based on uncertainty-weighted observations and reflect the pattern most available to the spawning season of pollock.*

-   Pursuing model configurations that can best be used to bridge between other software platforms (e.g., WHAM, stock synthesis, and AMAK2). *This practice should help confirm that the model used for this assessment concurs with other software (e.g., as found in* @li2021) *and also provide flexibility and transparency in transferring this assessment to future analysts.*

For future work we plan to adopt a new method for estimating cohort and year random effects for processes such as selectivity (e.g. @cheng2023).

<!-- References to add @oleary2022 @cheng2023 @maunder2020 @siskey2022 @indivero2023 @bailey2000  @defilippo2023 -->

We also set up application of @trijoulet2023 and this is provided [here](https://afsc-assessments.github.io/ebswp/articles/EBS_wp_OSA.html "OSA") (https://afsc-assessments.github.io/ebswp/articles/EBS_wp_OSA.html)

# References

::: {#refs}
:::

<!-- {{< pagebreak >}} -->

<!-- ![One-step ahead residual plots, see @trijoulet2023 for details.](doc/figs/EBSpollock_tuned_ats_age.pdf){#fig-OSAats} -->

<!-- ::: callout-warning -->

<!-- Here is a warning. -->

<!-- ::: -->

<!-- ::: callout-note -->

<!-- shit Note that there are five types of callouts, including: `note`, `tip`, `warning`, `caution`, and `important`. -->

<!-- ::: -->

<!-- ::: callout-important -->

<!-- Note that there are five types of callouts, including: `note`, `tip`, `warning`, `caution`, and `important`. -->

<!-- ::: -->

<!-- @francis2017fr provides -->

<!-- and this @afscpro2011 -->

<!-- @defilippo2023 -->