Ex 10: Operating models and MSE

Source: vignettes/ex10_om_mse.Rmd

1. Background

This is the 10th WHAM example. We assume you already have wham installed and are relatively familiar with the package. If not, read the Introduction and Tutorial.

In this vignette we show how to:

  • Make a default wham input file without an ASAP3 dat file.
  • Make an operating model, simulate data, and fit models
  • Perform a simple management strategy evaluation (MSE) using an operating model with Beverton-Holt stock recruitment and a SCAA estimating model
  • Make some plots

2. Setup 

# devtools::install_github("timjmiller/wham", dependencies=TRUE)
library(wham)
library(ggplot2)
library(tidyr)
library(dplyr)

Create a directory for this analysis:

# choose a location to save output, otherwise will be saved in working directory
write.dir <- "choose/where/to/save/output"
dir.create(write.dir)
setwd(write.dir)

3. A simple operating model and two estimation models

Make a basic_info list of input components defining a simple default stock. We’ll then pass this to prepare_wham_input and fit_wham.

make_digifish <- function(years = 1975:2014) {
    digifish = list()
    digifish$ages = 1:10
    digifish$years = years
    na = length(digifish$ages)
    ny = length(digifish$years)

    digifish$n_fleets = 1
    digifish$catch_cv = matrix(0.1, ny, digifish$n_fleets)
    digifish$catch_Neff = matrix(200, ny, digifish$n_fleets)
    digifish$n_indices = 1
    digifish$index_cv = matrix(0.3, ny, digifish$n_indices)
    digifish$index_Neff = matrix(100, ny, digifish$n_indices)
    digifish$fracyr_indices = matrix(0.5, ny, digifish$n_indices)
    digifish$index_units = rep(1, length(digifish$n_indices)) #biomass
    digifish$index_paa_units = rep(2, length(digifish$n_indices)) #abundance
    digifish$maturity = t(matrix(1/(1 + exp(-1*(1:na - na/2))), na, ny))

    L = 100*(1-exp(-0.3*(1:na - 0)))
    W = exp(-11)*L^3
    nwaa = digifish$n_indices + digifish$n_fleets + 2
    digifish$waa = array(NA, dim = c(nwaa, ny, na))
    for(i in 1:nwaa) digifish$waa[i,,] = t(matrix(W, na, ny))

    digifish$fracyr_SSB = rep(0.25,ny)
    digifish$q = rep(0.3, digifish$n_indices)
    digifish$F = matrix(0.2,ny, digifish$n_fleets)

    digifish$selblock_pointer_fleets = t(matrix(1:digifish$n_fleets, digifish$n_fleets, ny))
    digifish$selblock_pointer_indices = t(matrix(digifish$n_fleets + 1:digifish$n_indices, digifish$n_indices, ny))
    return(digifish)
}
digifish = make_digifish()

Now define other components needed by prepare_wham_input (selectivity and M).

selectivity = list(model = c(rep("logistic", digifish$n_fleets),rep("logistic", digifish$n_indices)),
    initial_pars = rep(list(c(5,1)), digifish$n_fleets + digifish$n_indices)) #fleet, index

M = list(initial_means = rep(0.2, length(digifish$ages)))

Here we specify that recruitment deviations are independent random effects, with no stock-recruit relationship.

NAA_re = list(N1_pars = exp(10)*exp(-(0:(length(digifish$ages)-1))*M$initial_means[1]))
NAA_re$sigma = "rec" #random about mean
NAA_re$use_steepness = 0
NAA_re$recruit_model = 2 #random effects with a constant mean
NAA_re$recruit_pars = exp(10)

Now we can make the input list with prepare_wham_input

input = prepare_wham_input(basic_info = digifish, selectivity = selectivity, NAA_re = NAA_re, M = M)

We can then define an operating model (OM) by simulating data (and recruitment) from this input:

om = fit_wham(input, do.fit = FALSE)

#simulate data from operating model
set.seed(8675309)
newdata = om$simulate(complete=TRUE)

Now put the simulated data in an input file with all the same configuration as the operating model.

temp = input
temp$data = newdata

Fit an estimation model that is the same as the operating model (self-test).

fit = fit_wham(temp, do.osa = FALSE, MakeADFun.silent = TRUE, retro.silent = TRUE)
fit$mohns_rho = mohns_rho(fit)
plot_wham_output(fit)

Set up a SCAA estimation model (recruitment as fixed effects). The only thing that is different is how the numbers at age are configured.

scaa_info = digifish
data_names = c("agg_catch", "catch_paa", "agg_indices","index_paa")
scaa_info[data_names] = newdata[data_names]

# recruitment as fixed effects
scaa_NAA_re = list(N1_pars = exp(10)*exp(-(0:(length(digifish$ages)-1))*M$initial_means[1]))
scaa_NAA_re$use_steepness = 0
scaa_NAA_re$recruit_model = 1

scaa_input = prepare_wham_input(basic_info = scaa_info, selectivity = selectivity, NAA_re = scaa_NAA_re, M = M, recruit_model = 1)

Fit the SCAA estimation model.

scaa_fit = fit_wham(scaa_input, do.osa = FALSE, MakeADFun.silent = TRUE, retro.silent = TRUE)
scaa_fit$mohns_rho = mohns_rho(scaa_fit)

4. A 2nd operating model with Beverton-Holt stock recruit relationship

Now make another operating model, bh_om, this time assuming a Beverton-Holt stock recruit relationship.

NAA_re = list(N1_pars = exp(10)*exp(-(0:(length(digifish$ages)-1))*M$initial_means[1]))
NAA_re$sigma = "rec" #random about mean
NAA_re$use_steepness = 1 #ok because M, WAA, etc are constant
NAA_re$recruit_model = 3 #Beverton-Holt
NAA_re$recruit_pars = c(0.5, exp(10))

# make input object for operating model
bh_input = prepare_wham_input(basic_info = digifish, selectivity = selectivity, NAA_re = NAA_re, M = M)

# make the operating model
bh_om = fit_wham(bh_input, do.fit = FALSE)

# simulate data and recruitment
set.seed(8675309)
sim_pop = bh_om$simulate(complete=TRUE)
temp = bh_input
temp$data = sim_pop

Fit an estimation model, bh_fit, that matches the operating model.

bh_fit = fit_wham(temp, do.osa = FALSE, MakeADFun.silent = TRUE, retro.silent = TRUE)
bh_fit$mohns_rho = mohns_rho(bh_fit)

Set up and fit a second estimation model with recruitment as fixed effects, i.e. with no stock recruit relationship, recruitment in each year is a separate parameter.

scaa_info = digifish
data_names = c("agg_catch", "catch_paa", "agg_indices","index_paa")
scaa_info[data_names] = sim_pop[data_names]

# recruitment as fixed effects
scaa_NAA_re = list(N1_pars = exp(10)*exp(-(0:(length(digifish$ages)-1))*M$initial_means[1]))
scaa_NAA_re$use_steepness = 0
scaa_NAA_re$recruit_model = 1 

scaa_input = prepare_wham_input(basic_info = scaa_info, selectivity = selectivity, NAA_re = scaa_NAA_re, M = M, recruit_model = 1)

# fit the scaa estimation model
scaa_fit = fit_wham(scaa_input, do.osa = FALSE, MakeADFun.silent = TRUE, retro.silent = TRUE)
scaa_fit$mohns_rho = mohns_rho(scaa_fit)

5. Closed-loop simulation study

Now we do a closed-loop simulation using the beverton-holt operating model. Catch advice is set at \(F_{40\%SPR}\) during the projection years.

First we use the prepare_projection function to set up an input that includes projection years. We set catch in the projection years to be nearly 0.

om_input <- prepare_projection(bh_om, proj.opts=list(n.yrs=39, use.last.F=FALSE, use.avg.F=FALSE,
                use.FXSPR=FALSE, proj.F=NULL, avg.yrs=NULL, proj.catch = rep(1,39),
                cont.ecov=TRUE, use.last.ecov=FALSE, avg.ecov.yrs=NULL, proj.ecov=NULL))

Then we use fit_wham with do.fit = FALSE to set up a TMB model that can be used to simulate a population time series and associated catch and index observations with the stochastic assumptions made in the NAA_re options that generated bh_om.

temp <- fit_wham(om_input, n.newton=n.newton, do.sdrep=F, do.retro=F, do.osa=F, do.check=F, do.proj=F,
  MakeADFun.silent = TRUE, save.sdrep=FALSE, do.fit = F)

Now generate a time series of the population. A key component is to keep the same seed for each simulation of the MSE.

set.seed(8675309)
pop_om_base_period = temp$simulate(complete=TRUE)

Replace the realized abundance parameters in the operating model input, then reset the operating model.

om_input$par$log_NAA = pop_om_base_period$log_NAA

# reset the om
om <- fit_wham(om_input, n.newton=n.newton, do.sdrep=F, do.retro=F, do.osa=F, do.check=F, do.proj=F,
  MakeADFun.silent = TRUE, save.sdrep=FALSE, do.fit = F)

The "temp" version of the operating model (red) just has the initial recruitment value, exp(10). The "om" has the full recruitment time series (black), generated using the specified stock-recruit relationship and nearly 0 catch.

plot(om$years_full, log(om$rep$NAA[,1]), type = 'l', xlab = "Year", ylab = "log(Recruits (1000s))")
lines(om$years_full, log(temp$rep$NAA[,1]), col = "red")

Define how often an assessment (SCAA) and catch advice will be made (3 years)

assess.interval = 3 #years.step = integer()
# yearly catch advice
advice = rep(1,39)

The first advice model will be completed at the end of the base period. The SCAA model defined as above will be updated every 3 years (assess.interval) of the projection/evaluation period.

scaa_step = fit_wham(scaa_input, do.osa = FALSE, MakeADFun.silent = TRUE, retro.silent = TRUE)

# the SCAA model with actual projections to make catch advice using F40
scaa_step_proj <- project_wham(scaa_step, proj.opts=list(n.yrs=5, use.last.F=FALSE, use.avg.F=FALSE,
            use.FXSPR=TRUE, proj.F=NULL, proj.catch=NULL, avg.yrs=NULL, avg.rec.yrs = tail(scaa_step$years,30),
            cont.ecov=TRUE, use.last.ecov=FALSE, avg.ecov.yrs=NULL, proj.ecov=NULL), save.sdrep=FALSE, MakeADFun.silent=TRUE)

# store the OM recruitment time series at each step in the loop
NAAtrack = matrix(NA, 79, 13)

# make catcj advice for the first 3 years of the evaluation period
advice[1:assess.interval] = rep(mean(scaa_step_proj$rep$pred_catch[scaa_step_proj$input$data$n_years_model + 1:5]), assess.interval)

The feedback loop will use the estimation model catch advice to fish at \(F_{40}\) for the 3-year projection period. This will reduce the stock in the OM, which reduces the realized recruitment. This plot will show the (lower) OM recruitment time series updated at the assessment interval.

for(y in seq(3,39,3)) {
    # set up projection of operating model
    om_input <- prepare_projection(bh_om, proj.opts=list(n.yrs=39, use.last.F=FALSE, use.avg.F=FALSE,
                use.FXSPR=FALSE, proj.F=NULL, proj.catch=advice, avg.yrs=NULL,
                cont.ecov=TRUE, use.last.ecov=FALSE, avg.ecov.yrs=NULL, proj.ecov=NULL))
    temp <- fit_wham(om_input, n.newton=n.newton, do.sdrep=F, do.retro=F, do.osa=F, do.check=F, do.proj=F,
    MakeADFun.silent = TRUE, save.sdrep=FALSE, do.fit = F)
    set.seed(8675309)
    updated_sim = temp$simulate(complete=TRUE)
    om_input$par$log_NAA = updated_sim$log_NAA
    NAAtrack[,y/3] = updated_sim$NAA[,1]
    lines(om$years_full, log(updated_sim$NAA[,1]), col = viridisLite::viridis(n= 39)[y], lty = 2)

    # increase terminal year
    scaa_info_step = make_digifish(min(digifish$years):(max(digifish$years)+y))
    scaa_info_step$agg_catch = rbind(updated_sim$agg_catch, updated_sim$agg_catch_proj[1:y,,drop=F])
    scaa_info_step$agg_indices = rbind(updated_sim$agg_indices, updated_sim$agg_indices_proj[1:y,,drop=F])
    scaa_info_step$catch_paa = abind::abind(updated_sim$catch_paa, updated_sim$catch_paa_proj[,1:y,,drop=F], along = 2)
    scaa_info_step$index_paa = abind::abind(updated_sim$index_paa, updated_sim$index_paa_proj[,1:y,,drop=F], along = 2)

    step_input = prepare_wham_input(basic_info = scaa_info_step, selectivity = selectivity, NAA_re = scaa_NAA_re, M = M)
    scaa_step = fit_wham(step_input, do.osa = FALSE, MakeADFun.silent = TRUE, retro.silent = TRUE)
    scaa_step_proj <- project_wham(scaa_step, proj.opts=list(n.yrs=5, use.last.F=FALSE, use.avg.F=FALSE,
                use.FXSPR=TRUE, proj.F=NULL, proj.catch=NULL, avg.yrs=NULL, avg.rec.yrs = tail(scaa_step$years,30),
                cont.ecov=TRUE, use.last.ecov=FALSE, avg.ecov.yrs=NULL, proj.ecov=NULL), save.sdrep=FALSE, MakeADFun.silent=TRUE)

    advice[y + 1:assess.interval] = rep(mean(scaa_step_proj$rep$pred_catch[scaa_step_proj$input$data$n_years_model + 1:5]), assess.interval)
}

Now compare the original OM simulated pop with \(F\) = 0 (purple) to the estimation model re-fit at 3-year time intervals providing catch advice equal to \(F_{40}\) (teal).

par(mfrow = c(2,2))
pal = viridisLite::viridis(n=3)
plot(temp$years_full, log(pop_om_base_period$NAA[,1]), type = 'l', xlab = "Year", ylab = "log(recruits)", col = pal[1])
lines(temp$years_full, log(NAAtrack[,13]), col = pal[2])
abline(v = max(temp$years), lty=2)

plot(temp$years_full, pop_om_base_period$SSB, type = 'l', xlab = "Year", ylab = "SSB", ylim = c(0,450000), col  = pal[1])
lines(temp$years_full, updated_sim$SSB, col = pal[2])
abline(v = max(temp$years), lty=2)

plot(temp$years_full, pop_om_base_period$pred_catch, type = 'l', xlab = "Year", ylab = "Catch (mt)", ylim = c(0,80000), col = pal[1])
lines(temp$years_full, updated_sim$pred_catch, col = pal[2])
points(temp$years_full, scaa_info_step$agg_catch, col = pal[2])
abline(v = max(temp$years), lty=2)

y.max = max(1.1*updated_sim$FAA_tot[,10]/exp(updated_sim$log_FMSY), na.rm=T)
plot(temp$years_full, updated_sim$FAA_tot[,10]/exp(updated_sim$log_FMSY), type = 'l', xlab = "Year", ylim=c(0,y.max), ylab = "F/FMSY", col = pal[2])
abline(v = max(temp$years), lty=2)
abline(h = 1, col = 'red', lty=2)