2.1 Setting up the file structure

The ss3sim package is set up assuming there is an established base-case operating model (OM) and estimation model (EM) to work with. The ss3sim package comes with three generic built-in model setups that can be used as is, modified, or replaced. (For more details, see the vignette on modifying these model setups or the vignette on creating your own OM and EMs). Each OM and EM should be in its own folder. The OM folder should have the files:

yourOMmodel.ctl
yourOMmodel.dat
ss3.par
starter.ss
forecast.ss

The EM folder should have:

yourEMmodel.ctl
starter.ss
forecast.ss

In both cases, nothing more and nothing less. The names of the .ctl and .dat files are not important. The package functions will rename them after they are copied to appropriate folders. These files should be derived from .ss_new files but with file extensions as shown above. It is important to use .ss_new files so the files have consistent formatting. Many of the functions in this package depend on that formatting.

To obtain .ss_new files, open a command window in the OM and EM folders and type ss3_24o_safe to run the models. Once the model is finished running, it will generate a number of files, including .ss_new files. Remove the .ss_new file extension from the files needed and add the appropriate file extension.

2.2 Cases and scenarios

The high-level wrapper function run_ss3sim uses unique case identifiers (IDs) that combine to create unique scenarios. The types of cases are: data quality (D), estimation (E), fishing mortality (F), retrospective (R), and any other letter describing a time varying case (e.g. M for natural mortality, S for selectivity, or G for growth). These case IDs are followed by an alphanumeric stock or species identifier (e.g. cod). The different versions of each case are identified with numbers. For example, the base-case scenario for a cod stock might be: D0-E0-F0-M0-cod. The order of the cases does not matter, as long as they are used consistently.

ss3sim relies on a set of plain text files to control the simulation. These plain text files are read by get_caseval and turned into argument lists that are passed to run_ss3sim. The function create_argfiles creates template input files. It reads the various functions and parses the arguments and default values into plain text files. The default settings create these files:

1. E0-spp.txt (for estimation method cases)

2. F0-spp.txt (for fishing mortality cases)

3. index0-spp.txt (controlled by the D (data) case)

4. agecomp0-spp.txt (controlled by the D (data) case)

5. lcomp0-spp.txt (controlled by the D (data) case)

6. X0-spp.txt (for a time-varying parameter X, where X could be any time-varying case letter)

Illustration of input and output folder and file structure for an ss3sim simulation. Folders are shown in blue, input files in orange, and output files in grey. All input and output files are in plain text format. Case files (orange files at bottom left) combine cases (e.g. M0 for a given natural mortality trajectory) with species or life-history OMs and EMs (e.g. cod-like or sardine-like). Alternatively, a user can skip setting up case files and specify the simulation cases directly in R code (see the vignette section on using ss3sim_base()

After running create_argfiles(), look in your working directory for the template input files. (You can use getwd() to see your current working directory or setwd() to set a different working directory.) Change the case ID number (defaults to 0) and the species identifier to a three letter identifier. For example, you might use cod, sar, or fla for cod, sardine, or flatfish. An example filename would therefore be M1-sar.txt or lcomp2-fla.txt. The case D1 corresponds to the files index1-spp.txt, agecomp1-spp.txt, and lcomp0-spp.txt. The other case types have single argument files.

The first column in the text files denotes the argument to be passed to a function. The second argument denotes the value to be passed. See the help for a change function to see the arguments that need to be declared. For example, see ?change_f.

You can use any R syntax to declare argument values. For example, c(1, 2, 4), or seq(1, 100). Character objects do not need to be quoted as long as they are one word, but can be if you would like. However, be careful not to use the delimiter (set up as a semicolon) anywhere else in the file besides to denote columns. You can add comments after any # symbol.

Putting that all together, below is what an example F1-cod.txt file might look like:

 years; 1913:2012
 years_alter; 1913:2012
 fvals; c(rep(0, 25), rep(0.114, 75))

2.3 Case file names

Function names and associated case files. Case file names are given as the prefix. For example a case file for sample_index might be index0-cod.txt.

change_tv() can take any other case file name, as long as the first line is function_type; change_tv.

2.4 Bias adjustment

Bias adjustment helps assure that the estimated log-normally distributed recruitment deviations are mean-unbiased leading to mean-unbiased estimates of biomass (Methot and Taylor 2011). The high-level wrapper function run_ss3sim allows users to specify whether or not they would like to use the bias adjustment routine built into the package by setting the argument bias_adjust to TRUE or FALSE. If TRUE, the function runs bias_nsim replicates for each scenario (located in the subfolder bias within a scenario folder) and then averages the bias-adjustment parameter estimates. The default number of bias-adjustment replicates is five. If a bias-adjustment replicate fails to converge (checked by identifying whether the Hessian was invertible), the parameter estimates from that replicate are ignored. A minimum threshold of 80% converged bias-adjustment replicates (adjustable via the conv_crit argument) is used to ensure reliable parameter estimates. The mean bias adjustment runs are then used in the EM for all subsequent “iterations” (i.e. replicates using same models but different process and observation errors) within that scenario. We assume bias adjustment parameters applicable to all iterations within the scenario are more likely to be found by averaging across multiple sets of parameters.

The bias adjustment process creates several files in the bias folder (for each scenario) which can be examined by the user. AdjustBias.DAT contains the calculated adjustment parameters for each bias replicate, and AvgBias.DAT contains the average of these estimates which are used in all subsequent simulations for that scenario. Individual plots for each successful bias adjustment replicate are also created. If bias_adjust is set to FALSE, no bias adjustment replicates are executed and no bias adjustment is performed.

2.5 Output file structure

Internally, the function copy_ss3models creates a folder structure and copies the operating and estimation models. The folder structure looks like:

  D0-E0-F0-M0-cod/1/om
  D0-E0-F0-M0-cod/1/em
  D0-E0-F0-M0-cod/2/om
  D0-E0-F0-M0-cod/2/em
  ...

The integer values after the scenario ID represent different iterations; the total number of iterations run are specified by the user. If you are using bias adjustment (bias_adjust = TRUE) then there will be some additional folders. In that case the folders will look like:

  D0-E0-F0-M0-cod/bias/1/om
  D0-E0-F0-M0-cod/bias/1/em
  D0-E0-F0-M0-cod/bias/2/om
  D0-E0-F0-M0-cod/bias/2/em
  ...
  D0-E0-F0-M0-cod/1/om
  D0-E0-F0-M0-cod/1/em
  D0-E0-F0-M0-cod/2/om
  D0-E0-F0-M0-cod/2/em
  ...

Note that the OM and EM folders will be renamed om and em within each iteration, the OM and EM are checked to make sure they contain the minimal files (as listed above), the filenames will be all lowercase, the data file is renamed ss3.dat, the control files are renamed om.ctl or em.ctl, and the starter and control files are adjusted to reflect these new file names. The functions in this package assume you have set your working directory in R to be the base folder where you will store the scenario folders.

3.1 Creating the case files

We will base the simulation around the base-case files created for a cod-like species in the papers Ono et al. (2015) and Johnson et al. (2015), both of which used the ss3sim package. You can refer to these papers for details on how the models were set up.

If we were starting from scratch, we would run the function create_argfiles(), which would create a default set of configuration files in our R working directory. Instead we will start with the configuration files from those papers.

To investigate the effect of different levels of precision on the survey index of abundance, we will manipulate the argument sds_obs that gets passed to sample_index. This argument ultimately refers to the standard error on the log(index) value, as defined in SS3. In case 0, we will specify the standard deviation at 0.1 and in case 1 we will increase the standard deviation to 0.4. We can do this by including the line: sds_obs; 0.1 in the file D0-cod.txt and the line: sds_obs; 0.4 in the file D1-cod.txt.

The file index0-cod.txt will therefore look like:

fleets;  2
years;   list(seq(1974, 2012, by = 2))
sds_obs; list(0.1)

fleets refers to the fleet number we want to sample from (as defined in the model), 2 in codOM.dat refers to the survey fleet. We then run our survey from 1974 to 2012 with a standard error on the log(survey index) of 0.1.

We will not describe the length or age composition sampling here in detail, but you can refer to the help files ?sample_lcomp and ?sample_agecomp along with the case files included in the package data. Also see vignette section where we describe the theory behind the age- and length-composition sampling in ss3sim.

To investigate the effect of fixing versus estimating M, we will manipulate the argument natM_val that gets passed to change_e. The first entry of natM_val is the value which M is fixed or initialized at and the second refers the phase. In case 0, we will set the phase in which M is estimated to -1 (any negative phase number will tell SS3 to not estimate a particular parameter) and use the argument NA to fix M at the true value (i.e. do not modify it from what is in the EM .ctl file). In case 1, we will estimate M in phase 3 and initialize the estimation of M at 0.20. We can do this by including the line: natM_val; c(NA, -1) in the file E0-cod.txt and the line: natM_val; c(0.20, 3) in the file E1-cod.txt. This means that case E1 will initialize M at 0.20 and estimate M in the third phase. For example, the complete case file E0-cod.txt is:

natM_type;          1Parm
natM_n_breakpoints; NULL
natM_lorenzen;      NULL
natM_val;           c(NA,-1)
par_name;           LnQ_base_3_CPUE
par_int;            NA
par_phase;          -1
forecast_num;       0

We will set the fishing mortality \(F\) to a constant level at \(F_\mathrm{MSY}\) (\(F\) at maximum sustainable yield) from when the fishery starts (in year 25) until the end of our simulation (year 100). Therefore, this is our F0-cod.txt:

years;       1913:2012
years_alter; 1913:2012
fvals;       c(rep(0,25), rep(0.114,75))

Although we do not have any time-varying processes in this example simulation, we include M case files in which we describe M as stationary (M0-cod.txt):

function_type; change_tv
param;         NatM_p_1_Fem_GP_1
dev;           rep(0, 100)

Simply by changing the deviations from a vector of zeros to any vector pattern of deviations, we could add time-varying M. In a [later vignette section]((#time-varying) we describe how time-varying parameters are unique in ss3sim and how you can make any parameter described in SS3 time varying using this approach.

3.2 Checking the case files

It is a good idea to check that the case files are manipulating the SS3 model files as intended. One way we can do this is by running the change functions directly on the SS3 model files using the arguments specified within the case files and inspecting the modified SS3 files. Another way is to run a single iteration of your simulation through run_ss3sim and carefully inspect the OM and EM model files that it creates. For example, we could run:

run_ss3sim(iterations = 1, scenarios =
  c("D0-E0-F0-M0-cod",
    "D1-E0-F0-M0-cod",
    "D0-E1-F0-M0-cod",
    "D1-E1-F0-M0-cod"),
  case_files = list(F = "F", D = c("index", "lcomp", "agecomp"),
    E = "E", M = "M"),
  case_folder = case_folder, om_dir = om,
  em_dir = em)

And then check the SS3 files that are created.

3.3 Self testing: deterministic simulations to check the models for bias

Self testing is a crucial step of any simulation study. One way to accomplish this is to set up an EM that is similar to the OM and include minimal error.

We will run some simulations to check our EM for bias when process and observation error are minimized. To do this, we will start by setting up a 100 row (number of years) by 20 column (number of iterations) matrix of recruitment deviations, where all values are set to zero.

recdevs_det <- matrix(0, nrow = 100, ncol = 20)

Then we will set up case “estimation” files in which the initialized values of the recruitment deviation standard deviations, SR_sigmaR, are set to the nominal level of 0.001. We will name these files E100-cod.txt and E101-cod.txt. In the case files, the key element is setting par_name = SR_sigmaR (the SS3 parameter name) and par_int = 0.001 (the initial value). The arguments par_name and par_int are set up to handle vectors of parameter names and values. In this example, changing SR_sigmaR will be the second parameter in the vector. Again, you can look at all the case files in the package folder /inst/extdata/eg-cases/. As an example, here is the file E101-cod:

natM_type;          1Parm
natM_n_breakpoints; NULL
natM_lorenzen;      NULL
natM_val;           c(0.20,3)
par_name;           LnQ_base_3_CPUE,SR_sigmaR
par_int;            c(NA,0.001)
par_phase;          c(-1,-1)
forecast_num;       0

To minimize observation error on the survey index, we will create a data case 100 that has a standard deviation on the survey observation error of 0.001. Therefore, our case file index100-cod.txt will look like this:

fleets;  2
years;   list(seq(1974, 2012, by = 2))
sds_obs; list(0.001)

When we run the simulations, we will pass our deterministic recruitment deviations to the function run_ss3sim. Running 20 iterations should be enough to identify whether our models are performing as we expect. Note that by default 5 bias adjustment runs are performed since we specify that bias_adjust=TRUE.

run_ss3sim(iterations = 1:20,
  scenarios = c("D100-E100-F0-M0-cod", "D100-E101-F0-M0-cod"),
  case_files = list(F = "F", D = c("index", "lcomp", "agecomp"),
    E = "E", M = "M"),
  case_folder = case_folder, om_dir = om, em_dir = em,
  bias_adjust = TRUE, user_recdevs = recdevs_det)

We have written out the scenario names in full for clarity, but ss3sim also contains a convenience function expand_scenarios. With this function we could instead write:

x <- expand_scenarios(list(D = 100, E = 100:101, F = 0, M = 0),
  species = "cod")
run_ss3sim(iterations = 1:20, scenarios = x,
  case_folder = case_folder, om_dir = om, em_dir = em,
  bias_adjust = TRUE, user_recdevs = recdevs_det,
  case_files = list(F = "F", D = c("index", "lcomp", "agecomp"),
    E = "E", M = "M"))

Note that due to the way that SS3 is being used as an operating model, you may see an ADMB error in the console:

Error -- base = 0 in function prevariable& pow(const prevariable& v1,
  CGNU_DOUBLE u)

However, this is not a problem since ADMB is not used to optimize the OM — the error can safely be ignored.

3.4 Running the stochastic simulations

The package contains a set of normally-distributed recruitment deviations, with a mean of -0.5 and a standard deviation of 1 (bias-corrected standard-normal deviations). To use the pre-specified deviations remove the argument user_recdevs from run_ss3sim. Process error will now be unique for each iteration but consistent across scenarios for a specific iteration, to facilitate comparison between scenarios.

We will only run 100 iterations here to save time and keep the package size down, but in a real simulation testing study you may want to run many more iterations, perhaps testing how many iterations are required before the key results stabilize.

run_ss3sim(iterations = 1:100, scenarios =
  c("D0-E0-F0-M0-cod",
    "D1-E0-F0-M0-cod",
    "D0-E1-F0-M0-cod",
    "D1-E1-F0-M0-cod"),
  case_files = list(F = "F", D = c("index", "lcomp", "agecomp"),
    E = "E", M = "M"),
  case_folder = case_folder, om_dir = om,
  em_dir = em, bias_adjust = TRUE)

3.5 Reading in the output and plotting the results

The function get_results_all reads in a set of scenarios

and combines the output into two .csv files: ss3sim_scalar.csv and ss3sim_ts.csv. The scalar file contains values for which there is a single estimated value (e.g. MSY) and the ts file refers to values for which there are time series of estimates available (e.g. biomass for each year). The column names refer to the output from SS3, so it can be useful to consult the SS3 manual for details.

get_results_all(user_scenarios =
  c("D100-E100-F0-M0-cod",
    "D100-E101-F0-M0-cod",
    "D0-E0-F0-M0-cod",
    "D1-E0-F0-M0-cod",
    "D0-E1-F0-M0-cod",
    "D1-E1-F0-M0-cod"))

We will read in the .csv files:

scalar_dat <- read.csv("ss3sim_scalar.csv")
ts_dat <- read.csv("ss3sim_ts.csv")

Or if you would like to follow along with the rest of the vignette without running the simulations above, you can load a saved version of the output:

data("ts_dat", package = "ss3sim")
data("scalar_dat", package = "ss3sim")

First, we will calculate some useful values in new columns and separate the deterministic from the stochastic simulation runs:

scalar_dat <- transform(scalar_dat,
  steep = (SR_BH_steep_om - SR_BH_steep_em)/SR_BH_steep_om,
  logR0 = (SR_LN_R0_om - SR_LN_R0_em)/SR_LN_R0_om,
  depletion = (depletion_om - depletion_em)/depletion_om,
  SSB_MSY = (SSB_MSY_em - SSB_MSY_om)/SSB_MSY_om,
  SR_sigmaR = (SR_sigmaR_em - SR_sigmaR_om)/SR_sigmaR_om,
  NatM = (NatM_p_1_Fem_GP_1_em - NatM_p_1_Fem_GP_1_om)/
     NatM_p_1_Fem_GP_1_om)

ts_dat <- transform(ts_dat,
  SpawnBio = (SpawnBio_em - SpawnBio_om)/SpawnBio_om,
  Recruit_0 = (Recruit_0_em - Recruit_0_om)/Recruit_0_om)
ts_dat <- merge(ts_dat, scalar_dat[,c("scenario", "replicate",
    "max_grad")])

scalar_dat_det <- subset(scalar_dat, E %in% c("E100", "E101"))
scalar_dat_sto <- subset(scalar_dat, E %in% c("E0", "E1"))
ts_dat_det <- subset(ts_dat, E %in% c("E100", "E101"))
ts_dat_sto <- subset(ts_dat, E %in% c("E0", "E1"))

Now we will turn the scalar data into long-data format so we can make a multipanel plot with the R package ggplot2.

scalar_dat_long <- reshape2::melt(scalar_dat[,c("scenario", "D", "E",
  "replicate", "max_grad", "steep", "logR0", "depletion", "SSB_MSY",
  "SR_sigmaR", "NatM")], id.vars = c("scenario", "D", "E",
  "replicate", "max_grad"))
scalar_dat_long <- plyr::rename(scalar_dat_long,
  c("value" = "relative_error"))

Now we can create boxplots of the deterministic model runs.

library("ggplot2")
p <- ggplot(subset(scalar_dat_long, E %in% c("E100", "E101") &
       variable != "SR_sigmaR"), aes(D, relative_error)) +
     geom_boxplot() +
     geom_hline(aes(yintercept = 0), lty = 2) +
     facet_grid(variable~E) +
     theme_bw() + ylim(-0.4, 0.4)
print(p)

Relative error box plots for deterministic runs. In case E100, M is fixed at the true value; in E101 we estimate M. In case D100, the standard deviation on the survey index observation error is 0.001.

Let’s look at the relative error in estimates of spawning biomass. We will colour the time series according to the maximum gradient. Small values of the maximum gradient (approximately 0.001 or less) indicate that model convergence is likely. Larger values (greater than 1) indicate that model convergence is unlikely. Results of individual iterations are jittered around the vertical axis to aid in visualization. The following three blocks of code produce Figures~, , and .

p <- ggplot(ts_dat_sto, aes(x = year)) + xlab("Year") +
    theme_bw() + geom_line(aes(y = SpawnBio, group = replicate,
    colour = max_grad), alpha = 0.3, size = 0.15) + facet_grid(D~E) +
    scale_color_gradient(low = "gray", high = "red")
print(p)

Time series of relative error in spawning stock biomass.

p <- ggplot(ts_dat_sto, aes(year, SpawnBio_em, group = replicate)) +
  geom_line(alpha = 0.3, aes(colour = max_grad)) + facet_grid(D~E) +
  scale_color_gradient(low = "darkgrey", high = "red") + theme_bw()
print(p)

Spawning stock biomass time series.

p <- ggplot(subset(scalar_dat_long, E %in% c("E0", "E1")),
       aes(D, relative_error)) +
     geom_boxplot() + geom_hline(aes(yintercept = 0), lty = 2) +
     facet_grid(variable~E) +
     geom_jitter(aes(colour = max_grad),
       position = position_jitter(height = 0, width = 0.1),
       alpha = 0.4, size = 1.5) +
     scale_color_gradient(low = "darkgrey", high = "red") +
     theme_bw()
print(p)

Relative error box plots for stochastic runs. In case E0, M is fixed at the true value; in E1 we estimate M. In case D1, the standard deviation on the survey index observation error is 0.4. In case D0, the standard deviation is quartered representing an increase in survey sampling effort.

9.1 Indices of abundance

The function sample_index facilitates generating relative indices of abundance (CPUE for fishery fleets). It samples from the biomass trends for the different fleets to simulate the collection of CPUE and survey index data. The user can specify which fleets are sampled in which years and with what amount of noise. Different fleets will “see” different biomass trends due to differences in selectivity, so each index is associated with the same fleet between the OM and EM. Since \(q=1\) for the OM, the indices are actually absolute indices of (spawning) biomass.

In practice, sampling from the abundance indices is relatively straightforward. The OM .dat file contains the annual biomass values for each fleet, and the sample_index function uses these true values as the mean of a distribution. The function uses a bias-corrected log-normal distribution with expected values given by the OM biomass and a user-provided standard deviation term that controls the level of variance.

More specifically, let

\[\begin{aligned} B_y&=\text{the true (OM) biomass in year $y$}\newline \sigma_y&=\text{the standard deviation provided by the user}\newline X\sim N(0, \sigma_y^2)&=\text{a normal random variable}\end{aligned}\]

Then the sampled value, \(B_y^\text{obs}\) is

\[B_y^\text{obs}=B_y e^{X-\sigma_y^2/2}\]

which has expected value \(\mathrm{E}\left[B_y^\text{obs}\right]=B_y\) due to the bias adjustment term \(\sigma_y^2/2\) (SR_sigmaR in SS3). This process generates log-normal values centered at the true value. It is possible for the user to specify the amount of uncertainty (e.g. to mimic the amount of survey effort), but it is not possible to induce bias in this process (as currently implemented).

9.2 Age and length compositions

The functions sample_agecomp and sample_lcomp sample from the true age and length compositions using either a multinomial or Dirichlet distribution. The user can further specify which fleets are sampled in which years and with what amount of noise (via sample size).

The following calculations demonstrate how ss3sim handles compositional data. They are shown for age compositions but apply equally to length compositions. The multinomial distribution \(\mathbf{m} \sim \text{MN}\left(\mathbf{p}, n, A\right )\) is defined as

\[\begin{aligned} \mathbf{m} &= m_1, \ldots, m_A\\ &=\text{the number of fish observed in bin $a$}\\ \mathbf{p}&= p_1, \ldots, p_A\\ &=\text{the true proportion of fish in bin $a$}\\ n&=\text{sample size} \\ A&=\text{number of age bins}\end{aligned}\]

In the case of SS3, actual proportions are input as data so \(\mathbf{m}/n\) is the distribution used instead of \(\mathbf{m}\). Thus, the variance of the estimated proportion for age bin \(a\) is \(\mathrm{Var}\left[M_a/n\right]=p_aq_a/n\). Note that \(m_a/n\) can only take on values in the set \(\{0, 1/n, 2/n, \ldots, 1\}\). With sufficiently large sample size this set approximates the real interval \([0,1]\), but the key point being that only a finite set of rational values of \(m_a/n\) are possible.

The multinomial, as described above, is based on assumptions of ideal sampling and is likely unrealistic, particularly with large sample sizes. One strategy to add realism is to allow for overdispersion.

9.3 Mean length at age

This function exists in the package, but has not been thoroughly tested or used in a simulation before. We caution users using this function to be extra cautious, and treat it as a beta version. Contact the developers for more information.

9.4 Conditional age at length

Conditional age-at-length (CAAL) data is an alternative to age compositions when there are paired age and length observations on the same fish. SS3 has this capability, and it has many appealing attributes. For instance, using CAAL data, instead of both age and length compositions of the same sampled fish, avoids including the same data twice. CAAL data are also expected to be more informative about growth than marginal age composition data although to date few studies have examined this data type (He et al. 2015, Monnahan et al. 2015). However, their implementation in ss3sim is significantly more complicated than simple age compositions (which are assumed independent of lengths).

A CAAL matrix, as implemented in SS3, is inputted the same way and in the same block in the .dat file as age compositions. Each row of this matrix corresponds to a length bin (for a year and fleet), while the columns remain ages and the cells are the number of fish. Thus, each row represents the observed age distribution for a length bin, conditioned on the fish lengths that were observed in the length compositions. Obviously for older or younger fish the age distribution will be truncated, and often many rows will be empty because no fish of that length bin were observed (they can be dropped from SS3 in this case).

Imagine sampling \(N\) fish for a year and fleet, taking lengths, and binning them to create the length distributions (numbers of fish in each length bin). From there, several strategies are possible for sampling ages from those fish: (1) age all fish, (2) take random subset of fish independent of length bin, or (3) take a fixed number of fish from each length bin. ss3sim can currently only handle option (2) above, and by extension (1). Future versions could include (3) but this is not currently implemented.

In any case, we need a hierarchy of sample sizes to implement this data. Consider a row of the CAAL matrix. The sample size for that row comes from the length compositions and represents the maximum number of fish that could be aged. If not all fish are aged, then a new sample size must be drawn that is strictly less than or equal to it. This will then be used to draw ages randomly from the expected values. If we consider all rows for a fleet and year (one for each length bin), then the sum of those will be the sample size for the CAAL data. If all fish are aged ((2) above), then no subsampling is done. However, if the CAAL sample size is less than the length, we need to be careful to not age more fish in a length bin than were in that length bin in the first place. We accomplish this in the code by doing sampling without replacement for vectors of length bins equal to the number of fish in them. This ensures realistic sampling. If the option (3) above were implemented, a different strategy would need to be implemented. For instance, if the user wants 10 fish from each length bin, but only 5 fish were observed, what to do?

A further complication is that when Dirichlet sampling is used for length compositions the number of fish observed will be real-valued, not whole fish. One cannot simply multiply by the length composition sample size to get whole numbers since they are real, and rounding or truncating would be unsatisfactory. Currently the function simply draws a multinomial sample from the length compositions of specified size (`Nsamp’). However, this does not gaurantee that fewer fish are aged than lengthed. If you are specifying a small number of fish to age relative to length, then this might be OK. However, we discourage the use of Dirichlet length samples when using CAAL data as currently implemented, unless you really understand the consequences of this issue.

We conclude by noting that the CAAL data are inherently tied to the length compositions (e.g., as if a trip measured lengths and ages for all fish), in contrast to the age compositions which were generated independently of the length data (e.g., as if one trip measured only ages and a second only lengths). This is an important distinction.

9.5 Empirical weight at age

This data type is very different than all the others. With this setting in SS3 all length-based processes are bypassed by assigning a weight to each age in each year, from which spawning biomass/fecundity can be determined. These values are not data in the sense they have a likelihood, but are generated from samples.

This sampling function was implemented for a specific study (Kuriyama et al. 2015), and may not be satisfactory for other studies. For this reason we will not explain the specifics of the function here and instead suggest contacting the developers if you wish to use empirical weight at age in your simulation.

13.1 For Unix (OS X and Linux)

To check if SS3 is in your path: open a Terminal window and type which ss3_24o_safe and hit enter. If you get nothing returned, then SS is not in your path. The easiest way to fix this is to move the SS3 binary to a folder that’s already in your path. To find existing path folders type echo $PATH in the terminal and hit enter. Now move the SS3 binary to one of these folders. For example, in a Terminal window type:

sudo cp ~/Downloads/SS3 /usr/bin/

You will need to use sudo and enter your password after to have permission to move a file to a folder like /usr/bin/.

Also note that you may need to add executable permissions to the SS3 binary after downloading it. You can do that by switching to the folder where you placed the binary (cd /usr/bin/ if you followed the instructions above), and running the command:

sudo chmod +x ss3_24o_safe

Check that SS3 is now executable and in your path:

which ss3_24o_safe

If you followed the instructions above, you will see the following line returned:

/usr/bin/ss3_24o_safe

If you have previously modified your path to add a non-standard location for the SS3 binary, you may need to also tell R about the new path. The path that R sees may not include additional paths that you have added through a configuration file like .bash_profile. If needed, you can add to the path that R sees by including a line like this in your .Rprofile file. (This is an invisible file in your home directory.)

Sys.setenv(PATH=paste(Sys.getenv("PATH"),"/my/folder",sep=":"))

13.2 For Windows

To check if SS is in your path for Windows 7 and 8: open a DOS prompt and type ss3_24o_safe -? and hit enter. If you get a line like ss3_24o_safe is not recognized... then SS3 is not in your path. To add the SS3 binary file to your path, follow these steps:

1. Find the correct version of the ss3_24o_safe.exe binary on your computer

2. Record the folder location. E.g. C:/SS3/

3. Click on the start menu and type environment

4. Choose Edit environment variables for your account under Control Panel

5. Click on PATH if it exists, create it if does not exist

6. Choose PATH and click edit

7. In the Edit User Variable window add to the end of the Variable value section a semicolon and the SS3 folder location you recorded earlier. E.g. ;C:/SS3. Do not overwrite what was previously in the PATH variable.

8. Restart your computer

9. Go back to the DOS prompt and try typing ss3_24o_safe -? and hitting return again.