`vignettes/ex1_SNEMA_yellowtail_flounder.Rmd`

`ex1_SNEMA_yellowtail_flounder.Rmd`

In this vignette we walk through an example using the `wham`

(WHAM = the Woods Hole Assessment Model) package to run a state-space age-structured stock assessment model. WHAM is a generalization of code written for Miller et al. (2016), and in this example we apply WHAM to the same stock as in Miller et al. (2016), Southern New England / Mid-Atlantic Yellowtail Flounder. Here, we demonstrate the basic `wham`

workflow:

Load

`wham`

and data-
Specify several (slightly different) models:

m1: statistical catch-at-age (SCAA) model, but with recruitment estimated as random effects; multinomial age-compositions

m2: as m1, but with logistic normal age-compositions

m3: full state-space model (numbers at all ages are random effects), multinomial age-compositions

m4: full state-space model, logistic normal age-compositions

Fit models and check for convergence

Compare models by AIC and Mohn’s rho (retrospective analysis)

Review plots of input data, diagnostics, and results.

We assume you have already read the Introduction and installed `wham`

and its dependencies. If not, you should be able to just run `devtools::install_github("timjmiller/wham", dependencies=TRUE)`

.

Open R and load the `wham`

package:

library(wham)

For a clean, runnable `.R`

script, look at `ex1_SNEMA_yellowtail_flounder.R`

in the `example_scripts`

folder of the `wham`

package install:

wham.dir <- find.package("wham") file.path(wham.dir, "example_scripts") #> [1] "/tmp/RtmpMNbw9P/temp_libpath28652852460d/wham/example_scripts"

You can run this entire example script with:

write.dir <- "choose/where/to/save/output" # otherwise will be saved in working directory source(file.path(wham.dir, "example_scripts", "ex1_SNEMA_yellowtail_flounder.R"))

Let’s create a directory for this analysis:

# choose a location to save output write.dir <- "/path/to/save/output" # modify dir.create(write.dir) setwd(write.dir)

WHAM was built by modifying the ADMB-based ASAP model code (Legault and Restrepo 1999), and is designed to take an ASAP3 .dat file as input. We generally assume in `wham`

that you have an existing ASAP3 .dat file. If you are not familiar with ASAP3 input files, see the ASAP documentation. For this vignette, an example ASAP3 input file is provided.

Copy `ex1_SNEMAYT.dat`

to our analysis directory:

wham.dir <- find.package("wham") file.copy(from=file.path(wham.dir,"extdata","ex1_SNEMAYT.dat"), to=write.dir, overwrite=FALSE)

Confirm you are in the working directory and it has the `ex1_SNEMAYT.dat`

file:

list.files() #> [1] "CPI.csv" "ex1_SNEMAYT.dat" "ex1_test_results.rds" #> [4] "ex2_SNEMAYT.dat" "ex2_test_results.rds" "ex4_test_results.rds" #> [7] "ex5_summerflounder.dat" "ex5_test_results.rds" "ex6_test_results.rds" #> [10] "GSI.csv"

Read the ASAP3 .dat file into R:

asap3 <- read_asap3_dat("ex1_SNEMAYT.dat")

We use the `prepare_wham_input()`

function to specify the model name and any settings that differ from the ASAP3 file. Our first model will use:

- recruitment model: random about mean, no S-R function (
`recruit_model = 2`

) - recruitment deviations: independent random effects (
`NAA_re = list(sigma="rec", cor="iid")`

) - selectivity: age-specific (fix sel=1 for age 5 in fishery, age 4 in index1, and age 2 in index2)

input1 <- prepare_wham_input(asap3, recruit_model=2, model_name="Ex 1: SNEMA Yellowtail Flounder", selectivity=list(model=rep("age-specific",3), re=rep("none",3), initial_pars=list(c(0.5,0.5,0.5,0.5,1,0.5),c(0.5,0.5,0.5,1,0.5,0.5),c(0.5,1,0.5,0.5,0.5,0.5)), fix_pars=list(5,4,2)), NAA_re = list(sigma="rec", cor="iid"))

The stock-recruit model options in WHAM are:

- 1 = random walk,
- 2 = random about mean (default),
- 3 = Beverton-Holt, and
- 4 = Ricker

Note that the parameterization of age-specific selectivity is specific to SNEMA yellowtail flounder. We will use age-specific selectivity parameters for the first three selectivity blocks. Selectivity blocks define a selectivity configuration that can be used in one or more fleets or surveys over any set of years in the model. `prepare_wham_input()`

fixes age-specific parameters at zero if there are any age classes without any observations where the selectivity block is applied. The function configures all other age-specific parameters to be estimated. Generally there may be confounding of the selectivity parameters with either fully-selected fishing mortality (for fleets) or catchability (for surveys), so selectivity parameters for one or more ages would need to be fixed to allow convergence. An initial fit with all selectivity parameters freely estimated can be useful in determining which age(s) to fix selectivity at 1. In the code above, we have already determined the ages to fix selectivity at 1 (age 5 in fishery, age 4 in index 1, and age 2 in index 2). If you are interested in more details and options for selectivity in WHAM, see Example 4 and `prepare_wham_input()`

.

m1 <- fit_wham(input1, do.osa = F) # turn off OSA residuals to save time in ex

By default, `fit_wham()`

uses 3 extra Newton steps to reduce the absolute value of the gradient (`n.newton = 3`

) and estimates standard errors for derived parameters (`do.sdrep = TRUE`

). `fit_wham()`

also does a retrospective analysis with 7 peels by default (`do.retro = TRUE`

, `n.peels = 7`

). For more details, see `fit_wham()`

.

We need to check that `m1`

converged (`m1$opt$convergence`

should be 0, and the maximum absolute value of the gradient vector should be < 1e-06). Convergence issues may indicate that a model is misspecified or overparameterized. To help diagnose these problems, `fit_wham()`

includes a `do.check`

option to run Jim Thorson’s useful function TMBhelper::Check_Identifiable(). `do.check = FALSE`

by default. To turn on, set `do.check = TRUE`

. See `fit_wham()`

.

```
#> stats:nlminb thinks the model has converged: mod$opt$convergence == 0
#> Maximum gradient component: 1.04e-07
#> Max gradient parameter: log_F1
#> TMB:sdreport() was performed successfully for this model
```

`m2`

-`m4`

The second model, `m2`

, is like the first, but changes all the age composition likelihoods from multinomial to logistic normal with estimated observation error variances:

input2 = input1 input2$data$age_comp_model_indices = rep(7, input2$data$n_indices) input2$data$age_comp_model_fleets = rep(7, input2$data$n_fleets) input2$data$n_age_comp_pars_indices = rep(1, input2$data$n_indices) input2$data$n_age_comp_pars_fleets = rep(1, input2$data$n_fleets) input2$par$index_paa_pars = rep(0, input2$data$n_indices) input2$par$catch_paa_pars = rep(0, input2$data$n_fleets) input2$map = input2$map[!(names(input2$map) %in% c("index_paa_pars", "catch_paa_pars"))] m2 <- fit_wham(input2, do.osa = F) # turn off OSA residuals to save time in ex

Check that `m2`

converged:

```
#> stats:nlminb thinks the model has converged: mod$opt$convergence == 0
#> Maximum gradient component: 2.06e-12
#> Max gradient parameter: catch_paa_pars
#> TMB:sdreport() was performed successfully for this model
```

The third, `m3`

, is a full state-space model where numbers at all ages are random effects (`NAA_re$sigma = "rec+1"`

):

input3 <- prepare_wham_input(asap3, recruit_model=2, model_name="Ex 1: SNEMA Yellowtail Flounder", selectivity=list(model=rep("age-specific",3), re=rep("none",3), initial_pars=list(c(0.5,0.5,0.5,0.5,1,0.5),c(0.5,0.5,0.5,1,0.5,0.5),c(0.5,1,0.5,0.5,0.5,0.5)), fix_pars=list(5,4,2)), NAA_re = list(sigma="rec+1", cor="iid")) m3 <- fit_wham(input3, do.osa = F) # turn off OSA residuals to save time in ex

Check that `m3`

converged:

```
#> stats:nlminb thinks the model has converged: mod$opt$convergence == 0
#> Maximum gradient component: 2.43e-07
#> Max gradient parameter: logit_q
#> TMB:sdreport() was performed successfully for this model
```

The last, `m4`

, is like `m3`

, but again changes all the age composition likelihoods to logistic normal:

input4 = input3 input4$data$age_comp_model_indices = rep(7, input4$data$n_indices) input4$data$age_comp_model_fleets = rep(7, input4$data$n_fleets) input4$data$n_age_comp_pars_indices = rep(1, input4$data$n_indices) input4$data$n_age_comp_pars_fleets = rep(1, input4$data$n_fleets) input4$par$index_paa_pars = rep(0, input4$data$n_indices) input4$par$catch_paa_pars = rep(0, input4$data$n_fleets) input4$map = input4$map[!(names(input4$map) %in% c("index_paa_pars", "catch_paa_pars"))] m4 <- fit_wham(input4, do.osa = F) # turn off OSA residuals to save time in ex

Check that `m4`

converged:

```
#> stats:nlminb thinks the model has converged: mod$opt$convergence == 0
#> Maximum gradient component: 2.86e-06
#> Max gradient parameter: logit_q
#> TMB:sdreport() was performed successfully for this model
```

Store all models together in one (named) list:

mods <- list(m1=m1, m2=m2, m3=m3, m4=m4)

Since the retrospective analyses take a few minutes to run, you may want to save the output for later use:

save("mods", file="ex1_models.RData")

We can compare models using AIC and Mohn’s rho:

res <- compare_wham_models(mods, fname="model_comparison", sort=TRUE)

```
#> dAIC AIC rho_R rho_SSB rho_Fbar
#> m4 0.0 -1466.9 0.3610 0.0091 -0.0106
#> m2 294.2 -1172.7 3.1632 -0.0715 -0.0184
#> m3 5574.0 4107.1 0.1287 0.0304 -0.0162
#> m1 6313.4 4846.5 0.8207 0.1905 -0.2322
```

res$best

`#> [1] "m4"`

By default, `compare_wham_models()`

sorts the model comparison table with lowest (best) AIC at the top, and saves it as `model_comparison.csv`

. However, in this example the models with alternative likelihoods for the age composition observations are not comparable due to differences in how the observations are defined. Still, the models that treat the age composition observations in the same way can be compared to evaluate whether stochasticity in abundances at age provides better performance, i.e. `m1`

vs. `m3`

(both multinomial) and `m2`

vs. `m4`

(both logistic normal).

Let’s do projections for the best model, `m4`

, using the default settings (see `project_wham()`

):

m4_proj <- project_wham(model=m4)

There are 3 options for plotting WHAM output. The default (`out.type='html'`

) creates and opens an HTML file with plots organized into tabs (code modified from `r4ss::SS_html()`

):

plot_wham_output(mod=m4_proj, out.type='html')

Setting `out.type='pdf'`

saves the plots organized into 6 `.pdf`

files corresponding to the tabs in the `.html`

file (Diagnostics, Input Data, Results, Reference Points, Retrospective, and Misc). Setting `out.type='png'`

saves the plots as `.png`

files organized into 6 subdirectories.

plot_wham_output(mod=m4_proj, out.type='png')

Many plots are generated—here we display some examples: