Merge pull request #102 from PIFSCstockassessments/FEAT-cookbook-reci…

…pe-articles

incorporating marcs edits

README.Rmd

@@ -26,7 +26,7 @@ evaluate a Stock Synthesis model. Diagnostics include residual analyses,
 hindcast cross-validation techniques, and retrospective analyses.
 Functions also allow users to reproduce the key model diagnostics plots
 that are presented in the paper ‘A Cookbook for Using Model Diagnostics in
-Integrated Stock Assessments’.
+Integrated Stock Assessments’ [(Carvalho et al. 2021)](https://www.sciencedirect.com/science/article/pii/S0165783621000874).
 
 The `ss3diags` Github repository provides step-by-step R recipes on how
 to:

vignettes/articles/Jitter.Rmd

@@ -18,7 +18,7 @@ library(tidyverse)
 Jitter analyses are commonly implemented in Stock Synthesis to ensure a model has reached global convergence. Jitter involves changing the parameter start values by a small increment and rerunning the model to see if that adjustment causes the model to converge at a lower likelihood. This can be useful for distinguishing if a model reached a local minimum or a global minimum. The number of jitter iterations should be anywhere between 50-100 to ensure a good spread of start values. If any of those runs has a lower likelihood than your current model, parameter start values should be adjusted to use those from the run with a lower likelihood. You can do this by adjusting the values in the control.ss file to match those in the ss.par_#_of_the_lower_likelihood_run. We provide the steps for running jitter analysis using `r4ss::jitter()` below. 
 
 ## Model inputs  
-To run a stock synthesis model, 4 input files are required: starter, forecast, control, and data. The input files for the example model can be found within the `ss3diags` package and accessed as shown below. Also, if you do not have `r4ss` installed, you will need to install for this tutorial. 
+To run a Stock Synthesis model, four input files are required: starter.ss, forecast.ss, control.ss, and data.ss. The input files for the example model can be found within the `ss3diags` package and accessed as shown below. Also, if you do not have the `r4ss` package installed, you will need to install it for this tutorial. 
 
 ```{r eval=FALSE}
 install_github("r4ss/r4ss")
@@ -124,7 +124,7 @@ converged_sum$likelihoods %>%
   )
 ```
 
-The figure plots the Total likelihood from each jitter run and the red line indicates the total likelihood value from the reference run. If there are any runs with points below the red line, use the parameter values from the run with the lowest likelihood value. 
+The figure plots the total likelihood from each jitter run and the red line indicates the total likelihood value from the reference run. If there are any runs with points below the red line, use the converged parameter values from the run with the lowest likelihood value as starting values in your base model. 
 We can also use `r4ss::SSplotComparisons()` to compare the spawning biomass trajectories between jitter runs to see the impact of different parameter values on the estimated quantities. 
 ```{r }
 SSplotComparisons(converged_sum,

vignettes/articles/Retrospective-Analysis.Rmd

@@ -13,10 +13,10 @@ knitr::opts_chunk$set(
 library(ss3diags)
 ```
 
-Retrospective analysis is commonly used to check the consistency of model estimates such as spawning stock biomass (SSB) and fishing mortality (F) as the model is updated with new data in retrospect. The retrospective analysis involves sequentially removing observations from the terminal year (i.e., peels), fitting the model to the truncated series, and then comparing the relative difference between model estimates from the full-time series with the truncated time-series. To implement a retrospective analysis with stock synthesis the `r4ss` package provides the `retro()` function. Here we provide a step-by-step example of how to run and analyze a retrospective analysis using a simple, example SS model.   
+Retrospective analysis is commonly used to check the consistency of model estimates such as spawning stock biomass (SSB) and fishing mortality (F) as the model is updated with new data. The retrospective analysis involves sequentially removing observations from the terminal year (i.e., peels), fitting the model to the truncated series, and then comparing the relative difference between model estimates from the full-time series with the truncated time-series. To implement a retrospective analysis with Stock Synthesis the `r4ss` package provides the `retro()` function. Here we provide a step-by-step example of how to run and interpret a retrospective analysis using a simple SS model.   
 
 ## Model inputs  
-To run a stock synthesis model, 4 input files are required: starter, forecast, control, and data. The input files for the example model can be found within the `ss3diags` package and accessed as shown below. Also, if you do not have `r4ss` installed, you will need to install for this tutorial. 
+To run a stock synthesis model, 4 input files are required: starter.ss, forecast.ss, control.ss, and data.ss. The input files for the example model can be found within the `ss3diags` package and accessed as shown below. Also, if you do not have the `r4ss` package installed, you will need to install it for this tutorial. 
 
 ```{r eval=FALSE}
 install.packages("pak")
@@ -41,21 +41,21 @@ r4ss::get_ss3_exe(dir = dir_retro, version = "v3.30.21")
 ```
 
 ## Retrospective Analysis
-Once you have the 4 input files and SS executable, you can run retrospective analysis as shown below. We are running it for 5 1-year peels, so with each run, *n* years of data are removed from the reference model and the model is re-run for a total of 5 times (i.e. peel 1 removes the last year of data, peel 2 removes the last 2 years of data, etc.) . The number of year peels can be adjusted with the `years` argument. If the SS executable file you are using is named something other than "ss3" (e.g. ss_opt_win.exe), you will need to specify this with the argument `exe = "ss_opt_win"`. Full documentation of the `retro()` function can be found on the [r4ss website](https://r4ss.github.io/r4ss/reference/retro.html).
+Once you have the 4 input files and the SS executable file, you can run a retrospective analysis as shown below. We are running it for five one-year peels, so with each run, one to five years of data are removed from the reference model and the model is re-run for a total of five times (i.e. peel 1 removes the last year of data, peel 2 removes the last 2 years of data, etc.) . The number of year peels can be adjusted with the `years` argument. If the SS executable file you are using is named something other than "ss3" (e.g. ss_opt_win.exe), you will need to specify this with the argument `exe = "ss_opt_win"`. Full documentation of the `retro()` function can be found on the [r4ss website](https://r4ss.github.io/r4ss/reference/retro.html).
 
 ```{r message=FALSE}
 r4ss::retro(dir = dir_retro, exe = "ss3", years = 0:-5, verbose = FALSE)
 ```
 ## Visualizing Output 
 
-To visualize the output and inspect for any patterns or biases, you need to load the report files into R and can use the `SSplotRetro()` function from `ss3diags`. The easiest way to load multiple report files is using `r4ss::SSgetoutput()` and `r4ss::SSsummarize()` functions. The default sub-directories for each peel, 0 to 5, are labeled `retro0` to `retro-5`.
+To visualize the outputted results and inspect them for any retrospective patterns, you will need to load the report files into R and use the `SSplotRetro()` function from `ss3diags`. The easiest way to load multiple report files is using `r4ss::SSgetoutput()` and `r4ss::SSsummarize()` functions. The default sub-directories for each peel, 0 to 5, are labeled `retro0` to `retro-5`.
 
 ```{r get_retro_files, message=FALSE, warning=FALSE}
 retro_mods <- r4ss::SSgetoutput(dirvec = file.path(dir_retro, "retrospectives", paste0("retro", seq(0, -5, by = -1))), verbose = F)
 retroSummary <- r4ss::SSsummarize(retro_mods, verbose = F)
 SSplotRetro(retroSummary, subplots = "SSB", add = TRUE)
 ```
-The default settings plot the spawning stock biomass time series for each peel, with the reference run (e.g. model with no years removed) as the "Ref" line and each successive peel as colored lines labeled by their end year. The solid line ends at the end year and the dashed line to the point shows the 1 year forecast. Displaying the projected SSB can help assess forecast bias. Note, forecasts are done automatically when using `r4ss::retro()` and are based on the settings in forecast.ss. The grey shaded area represents the 95% confidence intervals of uncertainty around the spawning biomass time series. Displayed in the center of the plot is the combined Mohn's $\rho$ for all retrospective runs, and in parentheses is the forecast Mohn's $\rho$. 
+The default settings plot the spawning stock biomass time series for each peel, with the reference run (i.e. original, full time series model) as the "Ref" line and each successive peel as colored lines labeled by their terminal year value. The solid line ends at the terminal year followed by the dashed line extending to a circle representing the one year forecast estimate. Displaying the projected SSB value helps in assessing forecast bias. Of note, forecasts are done automatically when using `r4ss::retro()` and are based on the settings in forecast.ss. The grey shaded area represents the 95% confidence intervals of uncertainty around the spawning biomass time series. Displayed in the center of the plot is the combined Mohn's $\rho$ for all retrospective runs, and in parentheses is the forecast Mohn's $\rho$. 
 
 ### Customizing the Plot  
 
@@ -74,7 +74,7 @@ retro2 <- SSplotRetro(retroSummary, subplots = "SSB", add = TRUE, uncertainty =
 retro3 <- SSplotRetro(retroSummary, subplots = "SSB", add = TRUE, uncertainty = FALSE, xlim = c(2015, 2022), forecast = FALSE)
 retro4 <- SSplotRetro(retroSummary, subplots = "SSB", add = TRUE, uncertainty = FALSE, xlim = c(2015, 2022), forecast = FALSE, showrho = FALSE, forecastrho = FALSE)
 ```
-Additionally, the fishing mortality can be plotted instead of spawning biomass by replacing `subplots = "SSB"` with `subplots = "F"`
+Additionally, fishing mortality can be plotted instead of spawning biomass by replacing `subplots = "SSB"` with `subplots = "F"`
 
 ### Summary Table  
 

vignettes/articles/aspm.Rmd

@@ -11,7 +11,7 @@ knitr::opts_chunk$set(echo = TRUE)
 The application of the Age-Structured Production Model (ASPM) approach as a diagnostic can help identify misspecification of the production function. If, in the absence of composition data (likelihood weighting set to 0), the ASPM fits well to the indices of abundance that have good contrast, then the production function is likely to drive the stock dynamics and indices will provide information about the absolute abundance ([Carvalho et al. 2017](https://www.sciencedirect.com/science/article/pii/S0165783616303113)). If there is not a good fit to the indices, then the catch data and production function alone cannot explain the trajectories of the indices of relative abundance. 
 
 ## Model inputs  
-To run a stock synthesis model, 4 input files are required: starter, forecast, control, and data. The input files for the example model can be found within the `ss3diags` package and accessed as shown below. Also, if you do not have `r4ss` installed, you will need to install for this tutorial. 
+To run a Stock Synthesis model, four input files are required: starter.ss, forecast.ss, control.ss, and data.ss. The input files for the example model can be found within the `ss3diags` package and accessed as shown below. Also, if you do not have the `r4ss` package installed, you will need to install it for this tutorial. 
 
 ```{r eval=FALSE}
 install.packages("pak")
@@ -39,7 +39,7 @@ r4ss::get_ss3_exe(dir = dir_tmp, version = "v3.30.21")
 ```
 
 ## ASPM 
-Once you have the 4 input files, you need to determine what components need to be turned off to run the ASPM. ASPM only depend on index of abundance and catch data, so any composition data, recruitment deviations, etc. need to be turned off. We provide an example that includes multiple types of data and recruitment deviations, however, the exact steps necessary for an individual model may vary depending on the complexity and components included. Therefore these steps may not be fully comprehensive for your model so be sure check what other components you may need to change.  
+Once you have the 4 input files, you need to determine what components need to be turned off to run the ASPM. ASPM only depend on index of abundance and catch data, so any composition data, recruitment deviations, etc. need to be turned off. We provide an example that includes multiple types of data and recruitment deviations, however, the exact steps necessary for an individual model may vary depending on the complexity and components included. Therefore these steps may not be fully comprehensive for your model so be sure to check what other components you may need to change.  
 Below, we show how to use the `r4ss` functions to make all the necessary changes to the `control.ss` and `ss.par` files.  
 
 ### Generate files  
@@ -93,7 +93,7 @@ SS_writestarter(starter,
 ```
 
 ### Control File 
-Some things in the control file that may need to be changed, are phases for parameters such as: 
+The estimation phases for the following parameters may need to be changed: 
 
 * length selectivity 
 * age selectivity

vignettes/articles/hcxval.Rmd

@@ -13,7 +13,7 @@ knitr::opts_chunk$set(
 library(ss3diags)
 ```
 
-Hindcast Cross-Validation is used to evaluate the prediction skill of a model by removing a portion of data, fitting the model to the subsequent time series, then projecting over the period that was removed and comparing the projected values to the observations that were omitted. Prediction skill is determined using the mean absolute scaled error (MASE) for input data (index of abundance, length composition, etc). In brief, the MASE score scales the mean absolute error (MAE) of forecasts (i.e., prediction residuals) to mean absolute error (MAE) of a naïve in-sample prediction, which is realized in the form of a simple 'persistence algorithm', i.e. index of abunduance will be the same next year as it is this year (see Eq. 3, p.5 in [Carvalho and Winker et al. 2021](https://www.sciencedirect.com/science/article/pii/S0165783621000874)). A MASE score > 1 indicates that the average model forecasts are worse than a random walk. Conversely, a MASE score < 1 indicates a level of prediction skill, for example, a MASE of 0.5 indicates that the model forecasts twice as accurately as a naïve baseline prediction; thus, the model has prediction skill.
+Hindcast cross-validation is used to evaluate the prediction skill of a model by removing a portion of data, fitting the model to the subsequent time series, then projecting over the period that was removed and comparing the projected values to the observations that were omitted. Prediction skill is determined using the mean absolute scaled error (MASE) for input data (index of abundance, length composition, etc). In brief, the MASE score scales the mean absolute error (MAE) of forecasts (i.e., prediction residuals) to the MAE of a naïve in-sample prediction, which is realized in the form of a simple 'persistence algorithm' where, for example, an index of abundance will be constant from one year to the next (see Eq. 3, p.5 in [Carvalho and Winker et al. 2021](https://www.sciencedirect.com/science/article/pii/S0165783621000874)). A MASE score > 1 indicates that the average model forecasts are worse than a random walk. Conversely, a MASE score < 1 indicates some prediction skills. For example, a MASE of 0.5 indicates that the model forecasts twice as accurately as a naïve baseline prediction; thus, the model has some prediction power.
 
 Implementing the Hindcast Cross-Validation (HCxval) diagnostic in Stock Synthesis requires the same model outputs generated when running a retrospective analysis. Therefore, no additional step is needed for HCxval if conducted in conjunction with retrospective analysis (see how to conduct a [restrospective analysis](https://pifscstockassessments.github.io/ss3diags/articles/Retrospective-Analysis.html)). For this example we will use the same output created in the retrospective example.
 
@@ -44,7 +44,7 @@ r4ss::sspar(mfrow = c(1, 2))
 SSplotHCxval(retroSummary, subplots = "cpue", add = TRUE)
 ```
 
-In the plots above, we see that for both fleets, the model has fairly good prediction skill (compared to a random walk) of index of abundance data. In the plots, the white points and white dashed line are the observed data that were included in the model with a truncated time series. The larger colored points are the observed data from each retrospective peel (i.e. data that was removed in that peel). The smaller colored point and dashed line show the model predicted value. The "Ref" line is the model run with the complete time series of data. The grey shaded areas represent the uncertainty of the data, with the darker portion indicating the portion that were included in the model and the lighter portion indicating which ones were removed and projected. The MASE scores displayed are the MASE and adjusted MASE in parentheses.  
+In the plots above, we see that, for both fleets, the model has fairly good prediction skills for index of abundance data (compared to a random walk). In the plots, the white points and white dashed line are the observed data that were included in the model with a truncated time series. The larger colored points are the observed data from each retrospective peel (i.e. data that was removed in that peel). The smaller colored points and dashed lines show the model predicted values. The "Ref" line is the model run with the complete time series of data. The grey shaded areas represent the uncertainty of the data, with the darker portion indicating the portion that were included in the model and the lighter portion indicating which ones were removed and projected. The MASE scores displayed are the MASE and adjusted MASE in parentheses.  
 
 ### Composition Data