4 Options for run_SSMSE
Many inputs are possible for the run_SSMSE() option. Here, we will describe some of the options available. For detailed documentation, type ?SSMSE::run_SSMSE() into the R console.
4.1 Scenarios in SSMSE
Note that multiple scenarios can be called in run_SSMSE(), often through vector inputs to run_SSMSE(). Below, we will describe inputs needed to run 1 scenario.
4.2 Operating model
The operating model (OM) for SSMSE should be a Stock Synthesis model. This could be any fitted Stock Synthesis model, except for models that use Empirical Weight at age.
There are 3 built-in models that comes with SSMSE: cod, growth_timevary, and SR_env_block. To use a built in model, specify the name of the model as part of the OM_name_vec. For example, to use the cod model in run_SSMSE for 1 scenario, set OM_name_vec = "cod". Otherwise, the path to the OM model should be specified in OM_in_dir_vec.
4.3 The Management Strategy and (if applicable) Estimation model (EM)
The management strategy (and EM) can be specified in one of two ways:
- Using an SS3 model and its forcast file to project catch by fleet
- Using a custom management strategy via a function in R
In theory, any management strategy should work, as long as it can take the data file produced by the OM as output and provide back to SSMSE future catches by fleet.
4.3.1 Specify the Management Strategy in a SS3 model
An SS3 model can be set up as the EM for the MSE. To use this option, specify "EM" as part of MS_vec. As with the OM, the built-in cod model could be used; just specify "cod" in the EM_name_vec. To use any other SS3 model as the EM, specify the path in EM_in_dir_vec.Note that models that use empirical weight at age can not yet be used as estimation models. Also, the .ss_new input files are used for an OM, while the original input files are used for the EM.
Future catches will be determined from the forecasting file settings of the SS3 model. SSMSE will change the number of forecast years to match the number of years between assessments, but other specifications need to be made by the user.
4.3.2 Using a custom management strategy/procedure
Users can outline a custom managment strategy as an R function to use. As long as the correct inputs and outputs are used, any estimation method and management procedure can be used. For example, here is a simple function that just sets future catches as half the sampled catches in a specified year:
constant_catch_MS <- function(OM_dat, nyrs_assess, catch_yr = 100,
frac_catch = 0.5, ...) { # need to include ... to allow function to work
# set catch the same as the previous year (sampled catch).
# catch is in the same units as the operating model, in this case it is in
# biomass.
catch <- data.frame(
year = (OM_dat$endyr + 1):(OM_dat$endyr + nyrs_assess), # the years to project the model forward
seas = 1, # hard coded from looking at model
fleet = 1, # hard coded from looking at model
catch = OM_dat$catch[OM_dat$catch$year == catch_yr, "catch"]*frac_catch,
catch_se = 0.05) # hard coded from looking at model
catch_bio <- catch # catch in biomass. In this case, catch is in biomass for both. Could also be left as NULL
catch_F <- NULL # catch in terms of F, can be left as NULL.
discards <- NULL # discards can be left as NULL if there are no discards
catch_list <- list(catch = catch,
catch_bio = catch_bio,
catch_F = catch_F,
discards = discards)
}Let’s assume this function is saved in a file within the working directory named constant_catch_MS.R.
This function can then be used in a call to run_SSMSE():
# define sample structure
datfile <- system.file("extdata", "models", "cod", "ss3.dat", package = "SSMSE")
sample_struct <- create_sample_struct(dat = datfile, nyrs = 6) # note warning
sample_struct$lencomp <- NULL # don't use length sampling
# run the SSMSE routine
run_result_custom <- run_SSMSE(
scen_name_vec = "constant-catch",
out_dir_scen_vec = "my_results",
iter_vec = 1,
OM_name_vec = "cod",
OM_in_dir_vec = NULL,
MS_vec = "constant_catch_MS", # use custom fun
custom_MS_source = "constant_catch_MS.R",
use_SS_boot_vec = TRUE,
nyrs_vec = 6,
nyrs_assess_vec = 3,
run_EM_last_yr = FALSE,
run_parallel = FALSE,
sample_struct_list = list(sample_struct),
seed = 12345)4.4 Sampling options
Currently, the only available sampling option is to use the bootstrapping module within SS3 itself. This means specifying use_SS_boot_vec = TRUE. Details on how sampling is done using the bootstrapping module in SS3 is available in the “Bootstrap Data Files” section of the SS3 user manual.
Users also need to specify how and which data types should be sampled for each future year in the simulation in sample_struct_list. sample_struct_list is a list of lists. The first level is a list for each scenario; then, for the scenario, there is a list of dataframes, each of which specifying which years and fleets of data should be sampled in the future as well as which standard errors or sample sizes should be used.
The helper function create_sample_struct() can be used to help users generate the list of dataframes for a scenario. See an example of this function’s use in the simple example or by typing ?SSMSE::create_sample_struct() into the R console.
Currently, users can sample data for catches (treated as fixed if sample_catch = FALSE for the scenario), CPUE (i.e., indices of abundance), length and age composition, conditional length at age compositions, mean size at age, and mean size. It is not yet possible to sample data for generalized size compositions, tagging data, and morph compositions.
Users can specify the sampling during the historical period of the model through the sample_struct_hist input to run_SSMSE. This is an optional input and has the same structure as sample_struct_list, but the years will be before the original end year of the OM.
4.5 Projecting the operating model
By default, SSMSE will simply extend the structure of the OM into the future using the same parameter values as in the last year of the model. However, the user may want to allow changes in parameter values to occur as the OM is extended forward in time. Users can input values into the future_om_list to accomplish this. This input is a list of lists. For the first level of lists, each represents a separate change to make. Within the first level, there are 4 list components that outline the details of the change to make. There are 2 main choices: 1) to specify model changes by telling SSMSE how values should be sampled or 2) to input your own custom values for future values of a given parameter.
4.5.1 The structure of future_om_list
For example, here is an example list containing just one future change for the model. It shows that the model should be changed to 4.5 in year 103 and afterwards:
[[1]]
[[1]]$pars
[1] "SizeSel_P_3_Fishery(1)"
[[1]]$scen
[1] "replicate" "scen2"
[[1]]$pattern
[1] "model_change"
[[1]]$input
first_yr_averaging last_yr_averaging last_yr_orig_val first_yr_final_val
1 NA NA 102 103
ts_param method value
1 mean absolute 4.5[1] 1[1] 4Note there is just one change specified here. For the change, there are four list items that are required for any specified change.
- The first list item is named “pars”. It contains a vector of parameter name(s) to apply the change to. The names should be the same as the names in
r4ss::SS_read_pars()or can include the following:
"rec_devs"- For specifying recruitment deviations"impl_error- For specifying implementation error in transforming a management procedure’s specified future catch into realized catch."all"- Apply the change to all parameters in the model, with the exception of “SR_sigmaR”, “SR_regime”, “SR_autocorr”, and “impl_error” parameters. Changing the first 3 stock recruitment related parameters would conflict with changes in recruitment deviations. Since implementation error is not a parameter specified in the control file and SSMSE does not rely on the implementation error parameter created through the forecasting file, including the implementation error in “all” did not seem appropriate.
In this case, we just want to apply the change to the SizeSel_P_3_Fishery(1) parameter of the model.
The second item is “scen”, which contains a vector of information about how to apply the changes within and across scenarios. The first value is an option to specify how the change should be applied among scenarios, either “randomize” (use a different value for each iteration of each scenario) or “replicate” (use the same set of values across scenarios for the same number iteration, but each value will be different across iterations within the same scenario). Note that the same values will be applied across iterations and scenarios if there isn’t any stochasticity in the values being drawn (e.g., standard deviation set at 0), regardless if “randomize” or “replicate” is chosen. In the example above, there is no stochasticity, so specifying randomize or replicate does not matter. Following the first value are the names of the scenarios to apply the change to. If the change should be applied to all of the scenarios, “all” can be used in place of naming every scenario. In this example, The change will only be applied to the scenario named “scen2”, so the input for “scen” is
c("randomize", "scen2").The “pattern” is a vector of character inputs. The first value should be either “model_change”, meaning that SSMSE will calculate the change values, or “custom” which allows the user to directly put in the values that they want in the model. In this case, “model_change” is used. A second input can also be added when the “model_change” pattern is used, which specifies the distribution to pull from. This could be “normal” or “lognormal”, but if no input is provided, a normal distribution will be used for sampling.
The fourth item, named “input”, is a dataframe, which contains different column name values depending on if the “model_change” pattern is used or if the “custom” pattern is used. For the model_change options, a dataframe in input specifies a change for the parameter of a distribution. Because we are using the model change option, we have the following columns:
first_yr_averaging: The first year to average historical values from the model, if using for the change. This should beNAif historical averaging will not be used.last_yr_averaging: The last year to average historical values from the model, if using for the change. This should beNAif historical averaging will not be used.last_yr_orig_val: The last year of the future deviations with the original model value. This value will not be changed, but the following year’s value will be.first_yr_final_val: The first year where the final value as specified will be reached. If no changes are made, the final value will continue forward in the OM through all projected years. For step changes, thefirst_yr_final_valis just 1 year after thelast_yr_orig_val. However, if a more gradual adjustment toward the final value is desired, this could be any number of years after the original value.- ts_param: The sampling parameter to modify. Options are “mean”, “ar_1_phi” (if using autocorrelation), “sd”, and “cv”. If not included in the dataframe, “mean” is assumed to be the same as in the previous model year, sd is assumed to be 0, and we assume no autocorrelation (as specified through an autoregressive AR1 process). Phi should be between -1 and 1 (this creates a stationary process) and if phi = 1, this creates a random walk (as random walk is a special case of an AR1 model). Note that both sd and cv cannot be specified within the same input dataframe.
- method: How to apply the change relative to the historical (if specified) or previous year’s value (if first_yr_averaging and last_yr_averaging for that row are NA). Options are “multiplicative”, “additive”, and “absolute”. This example uses “absolute”, meaning that the number in “value” is directly used.
- value. The value of the parameter change. This may be NA if historical averaging is used.
4.5.2 Example of future_om_list denoting a gradual change
Suppose we wanted to change the value of “SizeSel_P_3_Fishery(1)” in scen2 over 4 years after year 102 to arrive at a value of 4.5 in year 106. We can do this by using my_future_om_list, but changing the value in the first row of first_yr_final_val to 106:
$pars
[1] "SizeSel_P_3_Fishery(1)"
$scen
[1] "replicate" "scen2"
$pattern
[1] "model_change"
$input
first_yr_averaging last_yr_averaging last_yr_orig_val first_yr_final_val
1 NA NA 102 106
ts_param method value
1 mean absolute 4.54.5.3 Example of future_om_list with random deviations
Suppose we now wanted the value of “SizeSel_P_3_Fishery(1)” to change randomly according to a normal distribution with a standard deviation of 0.1 around a mean of 4.5 from 104 onwards. This can be done by adding a line specifying a change in standard deviation (which for now, has been assumed to be 0) to the data frame:
new_vals <- data.frame(first_yr_averaging = NA,
last_yr_averaging = NA,
last_yr_orig_val = 103,
first_yr_final_val = 104,
ts_param = "sd",
method = "absolute",
value = 0.1)
my_future_om_list[[1]][["input"]] <- rbind(my_future_om_list[[1]][["input"]],
new_vals)
my_future_om_list[[1]]$pars
[1] "SizeSel_P_3_Fishery(1)"
$scen
[1] "replicate" "scen2"
$pattern
[1] "model_change"
$input
first_yr_averaging last_yr_averaging last_yr_orig_val first_yr_final_val
1 NA NA 102 106
2 NA NA 103 104
ts_param method value
1 mean absolute 4.5
2 sd absolute 0.1Note that the last_yr_orig_val and first_yr_final_val are different than the line for the mean, which is allowed.
4.5.4 Example of using historical values for determining parameter values
This example applies random annual deviations to all parameters for scenarios scen2 and scen3.
future_om_list_2 <- vector(mode = "list", length = 1)
future_om_list_2 <- lapply(future_om_list_2, function (x) x <- vector(mode = "list", length = 4))
names(future_om_list_2[[1]]) <- c("pars", "scen", "pattern", "input")
future_om_list_2[[1]][["pars"]] <- "all"
future_om_list_2[[1]][["scen"]] <- c("randomize", "scen2", "scen3")
future_om_list_2[[1]][["pattern"]] <- "model_change" # defaults to using normal dist
future_om_list_2[[1]][["input"]] <- data.frame(first_yr_averaging = c(1, 1),
last_yr_averaging = c(100, 100),
last_yr_orig_val = c(100, 100),
first_yr_final_val = c(101, 101),
ts_param = c("cv", "mean"),
method = c("absolute", "multiplier"),
value = c(0.1, 1))Note that the choice of year range for historical values does not matter unless the parameter is already time-varying in the original operating model or has become time varying through a previous change. Otherwise, the base model parameter will be applied. If no historical years are included, then the base parameter value will be the basis of comparison for relative changes (i.e., method = multiplier or additive).
4.5.5 Example using “custom” pattern instead of “model_change”
custom_future_om_list <- create_future_om_list(example_type = "custom",
list_length = 1)
custom_future_om_list[[1]]
[[1]]$pars
[1] "impl_error"
[[1]]$scen
[1] "randomize" "all"
[[1]]$pattern
[1] "custom"
[[1]]$input
par scen iter yr value
1 impl_error scen1 1 101 1.05
2 impl_error scen1 2 101 1.05
3 impl_error scen1 3 101 1.05
4 impl_error scen1 4 101 1.05
5 impl_error scen1 5 101 1.05
6 impl_error scen1 1 102 1.05
7 impl_error scen1 2 102 1.05
8 impl_error scen1 3 102 1.05
9 impl_error scen1 4 102 1.05
10 impl_error scen1 5 102 1.05
11 impl_error scen1 1 103 1.05
12 impl_error scen1 2 103 1.05
13 impl_error scen1 3 103 1.05
14 impl_error scen1 4 103 1.05
15 impl_error scen1 5 103 1.05
16 impl_error scen1 1 104 1.05
17 impl_error scen1 2 104 1.05
18 impl_error scen1 3 104 1.05
19 impl_error scen1 4 104 1.05
20 impl_error scen1 5 104 1.05
21 impl_error scen1 1 105 1.05
22 impl_error scen1 2 105 1.05
23 impl_error scen1 3 105 1.05
24 impl_error scen1 4 105 1.05
25 impl_error scen1 5 105 1.05
26 impl_error scen1 1 106 1.05
27 impl_error scen1 2 106 1.05
28 impl_error scen1 3 106 1.05
29 impl_error scen1 4 106 1.05
30 impl_error scen1 5 106 1.05
31 impl_error scen2 1 101 1.10
32 impl_error scen2 2 101 1.10
33 impl_error scen2 3 101 1.10
34 impl_error scen2 4 101 1.10
35 impl_error scen2 5 101 1.10
36 impl_error scen2 1 102 1.10
37 impl_error scen2 2 102 1.10
38 impl_error scen2 3 102 1.10
39 impl_error scen2 4 102 1.10
40 impl_error scen2 5 102 1.10
41 impl_error scen2 1 103 1.10
42 impl_error scen2 2 103 1.10
43 impl_error scen2 3 103 1.10
44 impl_error scen2 4 103 1.10
45 impl_error scen2 5 103 1.10
46 impl_error scen2 1 104 1.10
47 impl_error scen2 2 104 1.10
48 impl_error scen2 3 104 1.10
49 impl_error scen2 4 104 1.10
50 impl_error scen2 5 104 1.10
51 impl_error scen2 1 105 1.10
52 impl_error scen2 2 105 1.10
53 impl_error scen2 3 105 1.10
54 impl_error scen2 4 105 1.10
55 impl_error scen2 5 105 1.10
56 impl_error scen2 1 106 1.10
57 impl_error scen2 2 106 1.10
58 impl_error scen2 3 106 1.10
59 impl_error scen2 4 106 1.10
60 impl_error scen2 5 106 1.10
61 impl_error scen3 1 101 1.10
62 impl_error scen3 2 101 1.10
63 impl_error scen3 3 101 1.10
64 impl_error scen3 4 101 1.10
65 impl_error scen3 5 101 1.10
66 impl_error scen3 1 102 1.10
67 impl_error scen3 2 102 1.10
68 impl_error scen3 3 102 1.10
69 impl_error scen3 4 102 1.10
70 impl_error scen3 5 102 1.10
71 impl_error scen3 1 103 1.10
72 impl_error scen3 2 103 1.10
73 impl_error scen3 3 103 1.10
74 impl_error scen3 4 103 1.10
75 impl_error scen3 5 103 1.10
76 impl_error scen3 1 104 1.10
77 impl_error scen3 2 104 1.10
78 impl_error scen3 3 104 1.10
79 impl_error scen3 4 104 1.10
80 impl_error scen3 5 104 1.10
81 impl_error scen3 1 105 1.10
82 impl_error scen3 2 105 1.10
83 impl_error scen3 3 105 1.10
84 impl_error scen3 4 105 1.10
85 impl_error scen3 5 105 1.10
86 impl_error scen3 1 106 1.10
87 impl_error scen3 2 106 1.10
88 impl_error scen3 3 106 1.10
89 impl_error scen3 4 106 1.10
90 impl_error scen3 5 106 1.10The inputs required for pattern = custom are similar to using pattern = model_change, but there is no need to specify a distribution as the second vector input of “pattern” and the column names in input are different. Also that the change is “randomized”, which indicates that a value for each scenario and each iteration are necessary, whereas if the first value of “scen” was “replicate”, separate values for each scenario should not be specified, but rather a single value for “all” should be used. The columns in the dataframe are:
- par: the parameter name or “all” if it should be applied to all parameters in the “pars” input
- scen: which scenario the value should apply to, or “all” if “replicate” is used as the scenario input
- iter: which number iteration the value should apply to. These are absolute iteration numbers.
- yr: the year the value applies to. These should be after the original end year of the OM.
- value: the value to apply to the model
4.5.6 How is the operating modified to accomodate changes as specified in the future_OM_list?
All changes are made by converting the parameter(s) with changes to be time varing by using additive annual parameter deviations. Because this is an operating model, Stock Synthesis is not run with estimation, so to get the changes into the OM, values (drawn or calculated in model_changes or specified by the user using “custom”) are directly input as parameter deviations into the ss.par file.
If an OM already contains time varying parameters, these parameters will also be converted to additive parameter deviations before applying the changes specified in the future_OM_list.
4.5.7 Example of specifying recruitment deviations
This shows some code for putting recruitment deviations with a mean of 0 and the same standard deviation as the historical recruitment deviations in years 1 to 100:
template_mod_change <- create_future_om_list(list_length = 1)
rec_dev_specify <- template_mod_change[[1]]
rec_dev_specify$pars <- "rec_devs"
rec_dev_specify$scen <- c("replicate", "all")
rec_dev_specify$input$first_yr_averaging <- 1
rec_dev_specify$input$last_yr_averaging <- 100
rec_dev_specify$input$last_yr_orig_val <- 100
rec_dev_specify$input$first_yr_final_val <- 101
rec_dev_specify$input$ts_param <- "sd"
rec_dev_specify$input$value <- NA
rec_dev_list <- list(rec_dev_specify)
rec_dev_list[[1]]
[[1]]$pars
[1] "rec_devs"
[[1]]$scen
[1] "replicate" "all"
[[1]]$pattern
[1] "model_change"
[[1]]$input
first_yr_averaging last_yr_averaging last_yr_orig_val first_yr_final_val
1 1 100 100 101
ts_param method value
1 sd absolute NA