Skip to content

Commit

Permalink
Merge branch 'FEAT-cookbook-recipe-articles' of https://github.com/PI…
Browse files Browse the repository at this point in the history 
…FSCstockassessments/ss3diags into FEAT-cookbook-recipe-articles
  • Loading branch information
@MOshima-PIFSC
MOshima-PIFSC committed Feb 9, 2024
2 parents ebc8fb1 + 6785f5e commit 986ae3f
Show file tree
Hide file tree
Showing 6 changed files with 150 additions and 131 deletions.
84 changes: 45 additions & 39 deletions vignettes/articles/Jitter.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,14 +25,13 @@ install_github("r4ss/r4ss")
```

```{r message=FALSE, warning=FALSE}
library(r4ss)
library(r4ss)
files_path <- system.file("extdata", package = "ss3diags")
dir_jitter <- file.path(tempdir(check = TRUE), "jitter")
dir.create(dir_jitter, showWarnings = FALSE, recursive = TRUE)
list.files(files_path)
file.copy(from = list.files(files_path, full.names = TRUE), to = dir_jitter)
```

You will need to make sure you have the [SS executable](https://github.com/nmfs-stock-synthesis/stock-synthesis) file either in your path or in the directory you are running the retrospective from (in this case `dir_jitter`). An easy way to get the latest release of stock synthesis is to use the `r4ss` function `get_ss3_exe()`.
Expand All @@ -44,31 +43,34 @@ r4ss::get_ss3_exe(dir = dir_jitter, version = "v3.30.21")
We will run the model in `dir_jitter` first to produce the necessary output files. It is recommended to do jitter runs in a subdirectory of your model run. This will keep all of the output files separate from other diagnostic tests you may run.
```{r}
r4ss::run(dir = dir_jitter, exe = "ss3", verbose = FALSE)
```

## Jitter
For this example, we will run 50 jitters. The `jitter()` function automates the entire process so you only need to give it a few arguments and it will run and produce the total likelihoods for each run. Full documentation of the `jitter()` function can be found at the [r4ss website](https://r4ss.github.io/r4ss/reference/jitter.html).

```{r run-jitter}
Njitter <- 50
jit.likes <- r4ss::jitter(dir = dir_jitter,
Njitter = Njitter,
init_values_src = 1,
jitter_fraction = 0.1,
exe = "ss3",
printlikes = FALSE,
verbose = FALSE)
Njitter <- 50
jit.likes <- r4ss::jitter(
dir = dir_jitter,
Njitter = Njitter,
init_values_src = 1,
jitter_fraction = 0.1,
exe = "ss3",
printlikes = FALSE,
verbose = FALSE
)
```

To analyze the output of all 50 runs, use `r4ss::SSgetoutput()` and `r4ss::SSsummarize()` as shown below.

```{r }
jit_mods <- SSgetoutput(keyvec = 0:Njitter, #0 to include reference run (Report0.sso)
getcomp = FALSE,
dirvec = dir_jitter,
getcovar = FALSE,
verbose = FALSE)
jit_mods <- SSgetoutput(
keyvec = 0:Njitter, # 0 to include reference run (Report0.sso)
getcomp = FALSE,
dirvec = dir_jitter,
getcovar = FALSE,
verbose = FALSE
)
jit_summary <- SSsummarize(jit_mods, verbose = FALSE)
```

Expand All @@ -85,18 +87,19 @@ You may also want to check that the models converged. To do this you can check t
# Maximum gradient
converged_grad <- which(jit_summary$maxgrad < 0.0001)
converged_ssb <- jit_summary$SpawnBio %>%
mutate(across(c(1:(Njitter+1)),
.fns = ~./replist0)) %>% # for each column, divide SSB by SSB from the reference run (replist)
select(-Label) %>%
pivot_longer(col = c(1:(Njitter+1)), names_to = "jitter", values_to = "SSB") %>%
pivot_wider(names_from = Yr, values_from = SSB) %>%
mutate(rownumber = seq(1, nrow(.))) %>%
tibble::column_to_rownames("jitter") %>%
filter_at(vars(1:(ncol(.)-1)), all_vars((.) < 2 & (.) > 0.5)) %>% #keep only rows where SSB is a reasonable value
select(rownumber) %>%
converged_ssb <- jit_summary$SpawnBio %>%
mutate(across(c(1:(Njitter + 1)),
.fns = ~ . / replist0
)) %>% # for each column, divide SSB by SSB from the reference run (replist)
select(-Label) %>%
pivot_longer(col = c(1:(Njitter + 1)), names_to = "jitter", values_to = "SSB") %>%
pivot_wider(names_from = Yr, values_from = SSB) %>%
mutate(rownumber = seq(1, nrow(.))) %>%
tibble::column_to_rownames("jitter") %>%
filter_at(vars(1:(ncol(.) - 1)), all_vars((.) < 2 & (.) > 0.5)) %>% # keep only rows where SSB is a reasonable value
select(rownumber) %>%
pull(rownumber)
converged_mods <- intersect(converged_grad, converged_ssb) #getting which models are in both groups
converged_mods <- intersect(converged_grad, converged_ssb) # getting which models are in both groups
converged_jitters <- jit_mods[converged_grad]
converged_sum <- SSsummarize(converged_jitters, verbose = FALSE)
```
Expand All @@ -105,26 +108,29 @@ converged_sum <- SSsummarize(converged_jitters, verbose = FALSE)
## Visualizing Output
To compare the likelihoods of all runs, we plot them as shown below. There are no built in functions (as of writing this vignette) in `r4ss` or `ss3diags` to generate a likelihood plot, therefore we provide code in the tidyverse syntax (using `ggplot2`) to visualize the results.
```{r }
converged_sum$likelihoods %>%
filter(str_detect(Label, "TOTAL")) %>%
select(-Label) %>%
pivot_longer(cols = everything(), names_to = "jitter", values_to = "likelihood") %>%
separate(jitter, into = c("replist", "jitter"), sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(jitter = as.numeric(jitter)) %>%
converged_sum$likelihoods %>%
filter(str_detect(Label, "TOTAL")) %>%
select(-Label) %>%
pivot_longer(cols = everything(), names_to = "jitter", values_to = "likelihood") %>%
separate(jitter, into = c("replist", "jitter"), sep = "(?<=[A-Za-z])(?=[0-9])") %>%
mutate(jitter = as.numeric(jitter)) %>%
ggplot(aes(x = jitter, y = likelihood)) +
geom_point(size = 3) +
geom_hline(aes(yintercept = likelihood[1]), color = "red") +
theme_classic() +
labs(y = "Total Likelihood",
x = "Jitter runs at a converged solution")
labs(
y = "Total Likelihood",
x = "Jitter runs at a converged solution"
)
```

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.
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,
subplots = 2,
ylimAdj=1,
new = FALSE)
SSplotComparisons(converged_sum,
subplots = 2,
ylimAdj = 1,
new = FALSE
)
```

14 changes: 5 additions & 9 deletions vignettes/articles/Retrospective-Analysis.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,14 +24,13 @@ pak::pkg_install("r4ss/r4ss")
```

```{r message=FALSE, warning=FALSE}
library(r4ss)
library(r4ss)
files_path <- system.file("extdata", package = "ss3diags")
dir_retro <- file.path(tempdir(check = TRUE), "retrospectives")
dir.create(dir_retro, showWarnings = FALSE)
list.files(files_path)
file.copy(from = list.files(files_path, full.names = TRUE), to = dir_retro)
```


Expand All @@ -46,17 +45,15 @@ Once you have the 4 input files and SS executable, you can run retrospective ana

```{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`.

```{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)
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$.

Expand All @@ -73,10 +70,9 @@ Examples of each of these changes are shown below, incrementally making each adj
```{r message=FALSE, warning=FALSE}
r4ss::sspar(mfrow = c(2, 2), plot.cex = 0.8)
retro1 <- SSplotRetro(retroSummary, subplots = "SSB", add = TRUE, uncertainty = FALSE)
retro2 <- SSplotRetro(retroSummary, subplots = "SSB", add = TRUE, uncertainty = FALSE, xlim = c(2015, 2022))
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)
retro2 <- SSplotRetro(retroSummary, subplots = "SSB", add = TRUE, uncertainty = FALSE, xlim = c(2015, 2022))
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"`

Expand Down
81 changes: 47 additions & 34 deletions vignettes/articles/aspm.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,10 +19,10 @@ pak::pkg_install("r4ss/r4ss")
```

```{r message=FALSE, warning=FALSE}
library(r4ss)
library(r4ss)
files_path <- system.file("extdata", package = "ss3diags")
dir_tmp <- tempdir(check = TRUE)
dir_tmp <- tempdir(check = TRUE)
dir_aspm <- file.path(dir_tmp, "aspm")
dir.create(dir_aspm, showWarnings = FALSE, recursive = TRUE)
list.files(files_path)
Expand Down Expand Up @@ -58,20 +58,24 @@ file.copy(from = file.path(dir_tmp, files), to = dir_aspm)
Set the recruitment devations in `ss.par` to 0.

```{r}
par <- SS_readpar_3.30(parfile = file.path(dir_aspm, "ss.par"),
datsource = file.path(dir_aspm, "data_echo.ss_new"),
ctlsource = file.path(dir_aspm, "control.ss_new"),
verbose = FALSE)
par <- SS_readpar_3.30(
parfile = file.path(dir_aspm, "ss.par"),
datsource = file.path(dir_aspm, "data_echo.ss_new"),
ctlsource = file.path(dir_aspm, "control.ss_new"),
verbose = FALSE
)
par$recdev1
par$recdev_forecast
par$recdev1[,"recdev"] <- 0
par$recdev1[, "recdev"] <- 0
#Would run if forecasts recdevs were not already 0
#par$recdev_forecast[,"recdev"] <- 0
SS_writepar_3.30(parlist = par,
outfile = file.path(dir_aspm, "ss.par"),
overwrite = T, verbose = FALSE)
# Would run if forecasts recdevs were not already 0
# par$recdev_forecast[,"recdev"] <- 0
SS_writepar_3.30(
parlist = par,
outfile = file.path(dir_aspm, "ss.par"),
overwrite = T, verbose = FALSE
)
```

### Starter File
Expand All @@ -81,10 +85,11 @@ starter <- SS_readstarter(file = file.path(dir_aspm, "starter.ss"), verbose = FA
starter$datfile <- "data_echo.ss_new"
starter$ctlfile <- "control.ss_new"
starter$init_values_src <- 1
SS_writestarter(starter,
dir = dir_aspm,
overwrite = TRUE,
verbose = FALSE)
SS_writestarter(starter,
dir = dir_aspm,
overwrite = TRUE,
verbose = FALSE
)
```

### Control File
Expand All @@ -105,30 +110,36 @@ Additionally, the likelihood for some components may need to be turned off as we
* initial F

```{r}
control <- SS_readctl_3.30(file = file.path(dir_aspm, "control.ss_new"),
datlist = file.path(dir_aspm, "data_echo.ss_new"),
verbose = FALSE)
control <- SS_readctl_3.30(
file = file.path(dir_aspm, "control.ss_new"),
datlist = file.path(dir_aspm, "data_echo.ss_new"),
verbose = FALSE
)
control$SR_parms
# Would need to run if PHASES were positive for "steep" and "sigmaR"
#control$SR_parms[c(2,3),"PHASE"] <- c(-4,-4)
# control$SR_parms[c(2,3),"PHASE"] <- c(-4,-4)
control$age_selex_parms
#would need to run if PHASES were positive for age selectivity
#control$age_selex_parms[,"PHASE"] <- control$age_selex_parms[,"PHASE"] * -1
control$size_selex_parms[,"PHASE"] <- control$size_selex_parms[,"PHASE"] * -1
# would need to run if PHASES were positive for age selectivity
# control$age_selex_parms[,"PHASE"] <- control$age_selex_parms[,"PHASE"] * -1
control$size_selex_parms[, "PHASE"] <- control$size_selex_parms[, "PHASE"] * -1
control$recdev_early_phase <- -4
control$recdev_phase <- -2
new_lambdas <- data.frame(like_comp = c(4,4,5,5,7,7,9,9,9,10),
fleet = c(1,2,1,2,1,2,1,2,3,1),
phase = rep(1, 10),
value = rep(0, 10),
sizefreq_method = rep(1, 10))
control$recdev_phase <- -2
new_lambdas <- data.frame(
like_comp = c(4, 4, 5, 5, 7, 7, 9, 9, 9, 10),
fleet = c(1, 2, 1, 2, 1, 2, 1, 2, 3, 1),
phase = rep(1, 10),
value = rep(0, 10),
sizefreq_method = rep(1, 10)
)
new_lambdas
control$lambdas <- new_lambdas
control$N_lambdas <- nrow(new_lambdas)
SS_writectl_3.30(control, outfile = file.path(dir_aspm, "control.ss_new"),
overwrite = TRUE, verbose = FALSE)
SS_writectl_3.30(control,
outfile = file.path(dir_aspm, "control.ss_new"),
overwrite = TRUE, verbose = FALSE
)
```

### Run ASPM
Expand All @@ -148,7 +159,9 @@ mods <- SSgetoutput(dirvec = c(
mods_sum <- SSsummarize(mods, verbose = FALSE)
SSplotComparisons(mods_sum, legendlabels = c("Ref", "ASPM"),
subplots = c(2,8,14), new = F)
SSplotComparisons(mods_sum,
legendlabels = c("Ref", "ASPM"),
subplots = c(2, 8, 14), new = F
)
```

Loading

0 comments on commit 986ae3f

Please sign in to comment.