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:

  1. Load wham and data

  2. 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

  3. Fit models and check for convergence

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

  5. Review plots of input data, diagnostics, and results.

1. Load data

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")

2. Specify model

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().

3. Fit model and check for convergence

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

Fit models 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")

4. Compare models

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).

Project the best model

Let’s do projections for the best model, m4, using the default settings (see project_wham()):

m4_proj <- project_wham(model=m4)

5. Plot input data, diagnostics, and results

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')
Example HTML file created by plot_wham_output(mod=m4, out.type='html').

Example HTML file created by plot_wham_output(mod=m4, 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:

Diagnostics

Index 2 4-panelLikelihood components

Conditional expected and posterior estimates of age-1 abundance (recruitment).Conditional expected and posterior estimates of age-5 abundance.

Fit to Index 1 age composition data.Residuals for Index 1 age composition.

One-step-ahead residuals Fleet 1 catch.One-step-ahead residuals Index 1 catch.One-step-ahead residuals Index 2 catch.

Input Data

CatchIndices of abundance

Results

SSB and FProportion of SSB at age

SSB-Rec on log-scaleSelectivity Fleet 1

Reference Points

Stock status (Kobe plot)Status relative to SPR

Retrospective

SSB with 7 peelsMohn’s rho (SSB) with 7 peels

Misc

Catch curves (Fleet 1 observed)Catch at age consistency (Fleet 1 observed)