Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Assorted cleanup and updating #19

Merged
merged 11 commits into from
Jan 31, 2025
46 changes: 26 additions & 20 deletions 01-data-sources.qmd

Large diffs are not rendered by default.

22 changes: 17 additions & 5 deletions 02-model-choices.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ The comments below are based on a team discussion on how to select a population

**Choice 3, Data length bins** should probably be big enough that the plus group doesn't have much more in it than the previous bin. We should check to make sure that the selection of length and age bins are consistent with each other. Typically, we often have max length bins where there are only a few fish, but a larger proportion of data in the data age plus group.

**Choice 4, Data age bins** should probably be big enough that the plus group doesn't have much more in it than the previous bin. In regards to population maximum age, there is no negative repercussions within SS3 for having a larger value, beyond a slower run time. Revising the age bins based on the data, rather than a priori rules about how to set this, may be considered "data mining". Could create a check within the survey and PacFIN comp codes that creates a flag when the data plus group has more than a certain percentage (i.e., 5\%). Also, add a warning to [{r4ss}](https://github.com/r4ss/r4ss) about the percentage in the population and data plus groups.
**Choice 4, Data age bins** should probably be big enough that the plus group doesn't have much more in it than the previous bin. In regards to population maximum age, there is no negative repercussions within SS3 for having a larger value, beyond a slower run time. Revising the age bins based on the data, rather than a priori rules about how to set this, may be considered "data mining". Could create a check within the survey and PacFIN comp codes that creates a flag when the data plus group has more than a certain percentage (i.e., 5%). Also, add a warning to [{r4ss}](https://github.com/r4ss/r4ss) about the percentage in the population and data plus groups.

## Determining prior for natural mortality

Expand All @@ -26,11 +26,11 @@ Owen's Advice on $M$; July 6, 2022

## Length-weight relationships

Use the [getWLpars()](https://pfmc-assessments.github.io/PacFIN.Utilities/reference/getWLpars.html) function in {PacFIN.Utilities} which provides estimates for the parameters required by Stock Synthesis. It is common to rely only on data from the WCGBT Survey to get the parameters used in the model under the assumption that this survey is most representative of the population as a whole. These parameters should always be fixed in the model as no data type is available in SS3 to accurately estimate them internally.
Use the [nwfscSurvey::estimate_weight_length()](https://pfmc-assessments.github.io/nwfscSurvey/reference/estimate_weight_length.html). It is common to rely only on data from the WCGBT Survey to get the parameters used in the model under the assumption that this survey is most representative of the population as a whole. These parameters should always be fixed in the model as no data type is available in SS3 to accurately estimate them internally. Other data sources (e.g., PacFIN if weights are available, etc.) can be passed to this function to estimate weight-length relationships. The `col_length` and `col_weight` arguments should be used to specify the specific string used to identify these quantities in your data.

## Estimating ageing error

See the [{nwfscAgeingError}](https://github.com/pfmc-assessments/nwfscAgeingError) package.
See the [{AgeingError}](https://github.com/pfmc-assessments/AgeingError) package.

## Maturity

Expand All @@ -46,11 +46,11 @@ Note: John Field has expressed concern that we are too focused on recent samples

[Dick et al. (2017)](http://www.sciencedirect.com/science/article/pii/S0165783616303745) has estimates for the $a$ and $b$ parameters in the functional form $F = aL^b$ for many rockfish. The estimates are based on length in mm and predicted number of eggs, fitted in log-space. If you use length in cm (like in SS3), and don't want huge spawning biomass values (eggs), you can convert the values in the paper to units of cm and millions, billions, or trillions of eggs. First, find the values "$\exp(a)$" (referred to as $a'$ below) and "$b$" from Table 6 in Dick et al. for your species (or subgenus, if no estimate for your species is reported). If your subgenus is not reported, you can use the "*Sebastes*" values.

The correct value of the "$a$" parameter in Stock Synthesis (using fecundity at length option 2: eggs=a*L^b) with length in cm and spawning output in **millions of eggs** is:
The correct value of the "$a$" parameter in Stock Synthesis (using fecundity at length option 2: eggs=a\*L\^b) with length in cm and spawning output in **millions of eggs** is:

$$a = \frac{a'\cdot10^b}{1000}$$

The division is by 1 thousand instead of 1 million because recruitment in SS3 is in thousands. The value of $b$ is unchanged, and can be used directly in the assessment.
The division is by 1 thousand instead of 1 million because recruitment in SS3 is in thousands. The value of $b$ is unchanged, and can be used directly in the assessment.

Other options include spawning output in **billions of eggs** via:

Expand Down Expand Up @@ -185,3 +185,15 @@ Below is detailed information regarding the meta-analysis approach applied betwe
- 2015: mean = 0.77, sd = 0.147

- 2017: mean = 0.72, sd = 0.158

## Projections (a.k.a. Forecasts)

The time-varying Pstar for your forecast.ss file can be calculated using `PEPtools::get_buffer()`,see [PEPtools github repo](https://github.com/pfmc-assessments/PEPtools/).

The default `sigma` input to that function is 0.5 for Category 1 stocks and 1.0 for Category 2 but if uncertainty in the model (one of `model$Pstar_sigma` or `model$OFL_sigma` as output from `r4ss::SS_output()`--check with SSC folks on which) is larger, use that instead. The `pstar` is set by the Council and the adopted value can be found online table [GMT015](https://reports.psmfc.org/pacfin/f?p=501:5301:6060050528090:INITIAL::::). It is recommended to get confirmation from the GMT representative as well.

Make sure you are using the right control rule for the PFMC in your forecast file: `3 # Control rule method (0: none; 1: ramp does catch=f(SSB), buffer on F; 2: ramp does F=f(SSB), buffer on F; 3: ramp does catch=f(SSB), buffer on catch; 4: ramp does F=f(SSB), buffer on catch)`

Fixed catches for first 2 years of forecast (the current year and next) should be provided by GMT rep to STAR Panel. The values are input at the bottome of the SS3 forecast file and should be split by fleet in the model. Adopted OFL/ACL values are in [GMT016 table](https://reports.psmfc.org/pacfin/f?p=501:1000:4357688221402:::::).

Creating the forecast table of OFLs and ABCs (ACLs if the relative stock size is below the management target) for your assessment document you can use [`r4ss::SSexecutivesummary()`](https://r4ss.github.io/r4ss/reference/SSexecutivesummary.html) to help with this task. You may want to specify the OFLs (and ACLs) for the next two years: `r4ss::SSexecutivesummary(model, forecast_ofl = c(3763, 3563))`
15 changes: 9 additions & 6 deletions 03-stock-synthesis.qmd
Original file line number Diff line number Diff line change
@@ -1,25 +1,28 @@
# Running Stock Synthesis

Note that there's lots more information on the SS3 website: <https://nmfs-ost.github.io/ss3-website> including the User Manual.

## Model setup

1. **Run the model without estimating anything at first.** The best way to do this is to run as `ss -stopph 0 -nohess` where `stopph` (short for "stop phase" is equivalent to setting the maximum phase in the starter file).

2. **Make sure you pay attention to any notes or warnings in "warning.sso".** If you don't understand the warnings, find out, but don't just ignore it.

3. **Debugging models that don't run**. The first place to look is "echoinput.sso". Start at the bottom and scan upwards until things start to look right or start at the top and scan downwards until things start to look wrong. It's often obvious when you have an extra input and things model starts to go awry. Use this information to fix your input files and try again. There are some additional debugging tips in the [Stock Synthesis User Manual](https://nmfs-stock-synthesis.github.io/doc/SS330_User_Manual.html#debugging-tips). Consider also reading the input files into R using `r4ss::SS_read()` which may help you find mismatched columns or bad values in the inputs.
3. **Debugging models that don't run**. The first place to look is "echoinput.sso". Start at the bottom and scan upwards until things start to look right or start at the top and scan downwards until things start to look wrong. It's often obvious when you have an extra input and things model starts to go awry. Use this information to fix your input files and try again. There are some additional debugging tips in the [Stock Synthesis User Manual](https://nmfs-ost.github.io/ss3-doc/SS330_User_Manual_release.html#debugging-tips). Consider also reading the input files into R using `r4ss::SS_read()` which may help you find mismatched columns or bad values in the inputs.

4. **Once the model runs, look at the ".ss_new" files.** These files contain rich comments and often better formatting than your own input files. They are also good for debugging, because sometimes a model will run, but the parameter lines are associated with different fleets, or have different roles than you expected. Check the parameter names on the right hand side of "control.ss_new" to make sure everything looks right. You can then either replace your input files with the ".ss_new" files or just copy and paste elements that you want to keep. Note that if you've estimated any parameters, then the initial values in "control.ss_new" have been updated to these estimates.

5. **Pull in parameter bounds that are way too wide** -- if you aren't anywhere near them during minimization, extremely wide bounds (like -15 to 15 on recruit deviations, or 3 to 31 for log-R0) just slow minimization and may result in poorer convergence properties.

## Model tuning

6. **Indices are typically tuned via the extra standard deviation parameter.** There are many reasons to expect that the input uncertainty values on indices of abundance are underestimates of the true uncertainty. Estimating an extra uncertainty parameter has worked well in a number of west coast groundfish assessments. However, the 2021 best practices document says, "STATs should be cautious to avoid adding variability to an index as a means of resolving model structure issues such as conflicts among data sources. Rather, variability should be added to account for sampling variance underestimating index uncertainty. *STATs should provide a priori reasons for why the index variability input to the model has been underestimated (or underspecified).*" Note that the extra SD parameter should reflect the observed variability in survey indices rather than poor fit to the observed trend in survey indices. Resist adding SD to surveys where there are trends in residuals without evidence of hyperdepletion or hyperstability, in which case a non-linear relationship between indices and stock size is more appropriate.
6. **Indices are typically tuned via the extra standard deviation parameter.** There are many reasons to expect that the input uncertainty values on indices of abundance are underestimates of the true uncertainty. Estimating an extra uncertainty parameter has worked well in a number of west coast groundfish assessments. However, the 2021 best practices document \[BUT REVISIONS TO THIS IN PREP FOR 2025\] says, "STATs should be cautious to avoid adding variability to an index as a means of resolving model structure issues such as conflicts among data sources. Rather, variability should be added to account for sampling variance underestimating index uncertainty. *STATs should provide a priori reasons for why the index variability input to the model has been underestimated (or underspecified).*" Note that the extra SD parameter should reflect the observed variability in survey indices rather than poor fit to the observed trend in survey indices. Resist adding SD to surveys where there are trends in residuals without evidence of hyperdepletion or hyperstability, in which case a non-linear relationship between indices and stock size is more appropriate.

7. **Composition data is typically tuned by either iterative Francis weighting or estimating Dirichlet-multinomial parameters.** The McAllister-Ianelli method has not performed as well in simulation testing. See the [data weighting section of the Stock Synthesis Manual](https://nmfs-stock-synthesis.github.io/doc/SS330_User_Manual.html#DataWeight) for more info (Dirichlet-multinomial guidance was updated by Jim Thorson in September 2022).
7. **Composition data is typically tuned by either iterative Francis weighting or estimating Dirichlet-multinomial parameters.** The McAllister-Ianelli method has not performed as well in simulation testing. See the [data weighting section of the Stock Synthesis Manual](https://nmfs-stock-synthesis.github.io/doc/SS330_User_Manual.html#DataWeight) for more info (Dirichlet-multinomial guidance was updated by Jim Thorson in September 2022). The choice between these two methods often depends on qualitative evaluation of fits to the comp data and index data and how it differs among tuning approaches.

8. **Discard ratios and mean body weight** These data SHOULD be tuned, but we don't typically do so. Kelli Johnson tuned these for Sablefish in 2019 with the following description:
8. **Discard ratios and mean body weight** These data SHOULD be tuned if they are used in the model, but haven't often done so. In August 2023 the `r4ss::SS_output()` function was augmented to return `$discard_tuning_info` and `$mnwgt_tuning_info` based on code written by Kelli Johnson for the 2019 sablefish assessment. The following description for tuning was provided in the 2019 assessment of sablefish:

> Added variances for discard rates and mean body weights were set using values calculated iteratively using the RMSE of differences between input and estimated values derived from SS3. Variances were parameterized in terms of standard deviation and coefficient of variation, respectively.
> Added variances for discard rates and mean body weights were set using values calculated iteratively using the RMSE of differences between input and estimated values derived from SS3. Variances were parameterized in terms of standard deviation and coefficient of variation, respectively.


9. **Think about sigmaR.** This could be an arbitrarily chosen value, freely estimated, or iteratively tuned. Methot and Taylor (2011) suggest a way that the tuning could be done. Whatever you choose, put a little thought into it. SigmaR should be greater than the SD of the estimated recruitment deviations. The table `$sigma_R_info` output by `r4ss::SS_output()` provides information on tuning sigmaR.
9. **Think about sigmaR.** This could be an arbitrarily chosen value, freely estimated, or iteratively tuned. [Methot and Taylor (2011)](https://cdnsciencepub.com/doi/10.1139/f2011-092) suggest a way that the tuning could be done. Whatever you choose, put a little thought into it. SigmaR should be greater than the SD of the estimated recruitment deviations. The table `$sigma_R_info` output by `r4ss::SS_output()` provides information on tuning sigmaR.
Loading