diff --git a/README.Rmd b/README.Rmd index fd6c16a..5bf1e0e 100644 --- a/README.Rmd +++ b/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: diff --git a/vignettes/articles/Jitter.Rmd b/vignettes/articles/Jitter.Rmd index db9974f..f322535 100644 --- a/vignettes/articles/Jitter.Rmd +++ b/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, diff --git a/vignettes/articles/Retrospective-Analysis.Rmd b/vignettes/articles/Retrospective-Analysis.Rmd index 7e29be5..65bb7bb 100644 --- a/vignettes/articles/Retrospective-Analysis.Rmd +++ b/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 diff --git a/vignettes/articles/aspm.Rmd b/vignettes/articles/aspm.Rmd index 019306f..86be0e4 100644 --- a/vignettes/articles/aspm.Rmd +++ b/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 diff --git a/vignettes/articles/hcxval.Rmd b/vignettes/articles/hcxval.Rmd index cb4d62d..aa748d4 100644 --- a/vignettes/articles/hcxval.Rmd +++ b/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 diff --git a/vignettes/articles/likelihood.Rmd b/vignettes/articles/likelihood.Rmd index c3facd3..be5bd7d 100644 --- a/vignettes/articles/likelihood.Rmd +++ b/vignettes/articles/likelihood.Rmd @@ -12,10 +12,10 @@ knitr::opts_chunk$set( ```{r setup} library(ss3diags) ``` -Likelihood profiling is a key model diagnostic that helps identify the influence of information sources on model estimates. R~0~ is a commonly profiled parameter because it represents a global scaling parameter. To conduct a profile, values of the parameter over a given range are fixed and the model re-run and then changes in total likelihood and data-component likelihoods are examined. It is recommended to run this throughout the development process, particularly after incorporating new data sets to understand how informative each component is on the estimation of R0 and if there is conflict between data sources. +Likelihood profiling is a key model diagnostic that helps identify the influence of information sources on model estimates. The number of recruits in an unfished population (R~0~) is a commonly profiled parameter since it dictates the absolute size of a fish population (serves as global scaling parameter). To conduct a profile, values of the parameter over a given range are fixed and the model re-run and then changes in total likelihood and data-component likelihoods are examined. It is recommended to run this throughout the development process, particularly after incorporating new data sets to understand how informative each component is on the estimation of R0 and if there is conflict between data sources. ## 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,8 +39,8 @@ You will need to make sure you have the [SS executable](https://github.com/nmfs- r4ss::get_ss3_exe(dir = dir_tmp, version = "v3.30.21") ``` -## R~0~ Profile -Once you have the 4 input files and SS executable, you can run a likelihood profile as shown below. The first step is to run the model you would like to do a profile for. We will do this in `dir_tmp` and then copy the necessary files to `dir_profile`. It's best to run the profile in a subdirectory of your main model run to keep the output files separate from other diagnostic tests you may run. +## Running a Profile on R~0~ +Once you have the four input files and SS executable file, you can run a likelihood profile as shown below. The first step is to run the model you would like to do a profile for. We will do this in `dir_tmp` and then copy the necessary files to `dir_profile`. It's best to run the profile in a subdirectory of your main model run to keep the output files separate from other diagnostic tests you may run. ```{r} r4ss::run(dir = dir_tmp, exe = "ss3", verbose = FALSE) @@ -48,7 +48,7 @@ files <- c("data.ss", "control.ss_new", "starter.ss", "forecast.ss", "ss.par", " file.copy(from = file.path(dir_tmp, files), to = dir_profile) ``` -Once you have the input files you need in `dir_profile` you will need to create a vector of values to profile across. The range and increments to choose depend on your model and the resolution at which you want to analyze the likelihood profile. For this example we will use a fairly course resolution to speed up total run time. +Once you have the input files in the `dir_profile` directory, you will need to create a vector of values to profile across. The range and increments to choose depend on your model and the resolution at which you want to analyze the likelihood profile. For this example we will use a fairly coarse resolution to speed up total run time. ```{r} CTL <- SS_readctl_3.30(file = file.path(dir_profile, "control.ss_new"), datlist = file.path(dir_profile, "data.ss")) @@ -68,7 +68,7 @@ START$ctlfile <- "control_modified.ss" SS_writestarter(START, dir = dir_profile, overwrite = TRUE, verbose = F) ``` -To run the profile, use `r4ss::profile()` and you will need to specify a partial string of the name of the parameter you are profiling over (in this case "SR_LN" will match with "SR_LN(R0)"), and the vector of values to profile across. The `newctlfile` is the control file that will be adjusted with values from the `profilevec` and can be named anything you prefer, it just needs to match what you put in the starter file for "cltfile". Full documentation of the `profile()` function can be found on the [r4ss website](https://r4ss.github.io/r4ss/reference/profile.html). +To run the profile, use `r4ss::profile()` and you will need to specify a partial string of the name of the parameter you are profiling over (in this case "SR_LN" will match with "SR_LN(R0)"), and the vector of values to profile across. The `newctlfile` is the control file that will be adjusted with values from the `profilevec` and can be named anything you prefer, it just needs to match what you put in the starter file for "ctlfile". Full documentation of the `profile()` function can be found on the [r4ss website](https://r4ss.github.io/r4ss/reference/profile.html). ```{r} profile( dir = dir_profile, @@ -81,7 +81,7 @@ profile( ``` ## Visualizing Output -Once the profile is finished running, we can visualize the results to determine if there is conflict between the data sources. If all data sources reach a minimum at the same R~0~ value, this indicates good agreement between them. However, more likely, one or more will be minimized at different R0 values from the global R~0~ value. This is a sign of conflict between your data sources and may require you to consider data weighting. +Once the profile is finished running, we can visualize the results to determine if there is conflict between the data sources. If all data sources reach a negative log-likelihood (NLL) minimum at the same R~0~ value, this indicates good agreement between them. However, more likely, one or more will reach this NLL minimum at different R0 values from the global R~0~ value. This is a sign of conflict between your data sources and may require you to consider data weighting. ```{r} profile_mods <- SSgetoutput(dirvec = dir_profile, keyvec = 1:length(r0_vec), verbose = FALSE) @@ -92,7 +92,7 @@ SSplotProfile(profile_mods_sum, ) ``` -The profile plot shows the changes in log-likelihood across the vector of values profiled over for the total likelihood and each of the contributing components. There is a minimum threshold that a component must contribute so if there is a data source that is in your model but does not show up in the plot, the contribution may not be large enough. The steepness of each trajectory indicates how informative (or not) that data source was. For example, the age data in the plot above is much steeper on the left side of the minimum R0 value than the index data, which suggests that age composition data is more informative in the model. +The profile plot above shows the changes in log-likelihood across the selected range of values for each of the contributing data components. There is an adjustable minimum contribution threshold that a component must meet to appear in this figure (if there is a data source that is in your model but does not show up in the plot, its contribution to the total likelihood may not be large enough). The steepness of each trajectory indicates how informative (or not) a data source is. For example, the age data in the plot above is much steeper on the left side of the minimum R0 value than the index data, which suggests that age composition data is more informative in the model. You can also plot data-type and fleet-specific profiles using `r4ss::PinerPlot()`. Below we are plotting the profile for the length composition data by fleet and the likelihood of survey data by fleet. This will allow us to see if there are conflicts and what sources are the main drivers.