Supplemental Data: Operating Model Equations and Components
Population Dynamics
The operating model uses the following equations to mimic the dynamics of an age structured population that is partitioned into two regional subpopulations. Parameter descriptions and values are shown in Table S.1.
Equilibrium number of age a fish in region r in the absence of fishing:
(1)
where is the equilibrium unfished number of age-0 recruits, is the equilibrium recruitment fraction to region r, with and , and M is the annual instantaneous rate of natural mortality, assumed constant over time and age.
Equilibrium spawning biomass (both regions) in the absence of fishing:
(2)
whereis the weight-at-age, based on the length-weight relationship
(3)
and length-at-age equation
(4)
andis the maturity-at-age, calculated as a logistic function of length-at-age.
Number of age a fish in region r in year t:
(5)
where is the number of age-0 recruits in year t, is the recruitment fraction to region r in year t, and is the annual instantaneous rate of fishing mortality at age a in region r during year t. The numbers-at-age for the initial year were generated to reflect randomness in recruitment and F during years prior to the simulated period.
Number of recruits in year t:
(6)
where is the spawning stock biomass in year t and is a normal random deviate with zero mean and standard deviation . Curvature in the relationship R(SSB) is controlled by the steepness parameter h.
Fishing mortality at age a in region r during year t:
(7)
where is the selection coefficient for age a fish, Et[X] denotes the expected value with respect to time tof random variable X, the expected values for the two sets of values follow annual deterministic patterns specified in the experimental design (illustrated in Figure S.1), and is a normal random deviate with zero mean and standard deviation .
Selection at age a:
(8)
Spawning stock biomass in year t:
(9)
where .
Relative spawning stock biomass in year t (Depletion):
(10)
Fraction of age-0 fish recruiting to region r in year t:
is driven by environmental variation plus randomness. For region one is
(11)
where is a normal random deviate with zero mean and standard deviation and for regiontwo is. If the experimental factor EnvirP is set at the gradual level, denoted as , there is gradual cyclical change in and
(12)
where
(13)
If the EnvirP factor is at the constant level (EC) the expected recruitment fraction values are constant over time and the value for region one is equal to RecFracinit. If EnvirP has the abrupt setting (EA), the expected recruitment fraction for region one is initially equal to 60% () and changes as a step function in year 12 to 30%, causing a reversal in which region receives the most recruitment. The different patterns of variation are illustrated in Figure S.1.
Catch (numbers) of age a fish in region r in year t:
(14)
Yield (catch weight) of age a fish in region r in year t:
(15)
Data Generation
The following equations are used to generate various types of data that are input to stock synthesis for analysis or compared with estimates output by stock synthesis. One dataset perfectly reflects the random variations in recruitment and exploitation that drive the population dynamics. Another dataset includes additional random variation due to sampling (observation error).
Variables with No Observation Error
Fishery catch (numbers) per unit effort index for region r in year t:
(16)
This equation follows from a rearrangement of the basic relationships between (a) catch (in numbers), fishing mortality and numerical abundance () and (b) fishing mortality and effort (), and the assumption that catchability (Q) is the same in both regions.
Survey catch (numbers) of age a fish from region r in year t:
(17)
where the survey selectivity coefficients, , are calculated using the same formula as the fishery selectivity coefficients but with a reduced age at 50% selection (Table S.1).
Survey biomass index for region r in year t:
(18)
where the Survey Numbers at age a for region r correspond to the population numbers at age for that region (eq. 5) but reduced by the survey selection coefficients.
Environmental index in year t that drives the recruitment fraction to region one.
(19)
Data Inputs to Stock Synthesis that Include Observation Error
The data series for fishery CPUE and Survey Biomass include lognormal observation error.
Fishery catch per unit effort index of abundance in year t:
(20)
where is a normal random deviate with zero mean and standard deviation 0.5. For experimental treatments that ignored spatial structure the regional observations were added together.
Survey biomass index for region r in year t:
(21)
where is a normal random deviate with zero mean and standard deviation 0.25. For experimental treatments that ignored spatial structure the regional observations were added together.
Fishery age composition observations for region r in year t:
The fishery age composition data include multinomial error with 400 observations of age from each region (no age reading error) and predicted proportions at age a for region r in year t given by
(22)
For experimental treatments that ignore spatial structure the observed regional proportions at age ain year twere combined using
(23)
where is the observed proportion of age a fish in the sample from region r in year t and is the total catch in numbers from region r, equal to . Here we make the unrealistic assumption that the Ca,r values are known without error.
Survey age composition observations for region r in year t:
The survey age composition data include multinomial error with 200 observations of age from each region (no age reading error) and predicted proportions at age a for region r in year t given by
(24)
For experimental treatments that ignore spatial structure the observed regional proportions at age ain year twere combined using
(25)
where is the observed proportion of age a fish in the sample from region r in year t and is the total survey catch in numbers from region r, equal to .
Supplemental Data: Evaluating Convergence
To provide some assurance that the stock synthesis parameter estimates corresponded to maximum likelihood estimates we employed the algorithm described using pseudo-code in Figure S.3. For each iteration of each experimental treatment and estimation model configuration we conducted sets of runs that started with initial parameter values that were “jittered” by adding to the true parameter values the product of three values: 0.1, a standard normal random deviate, and the difference between each parameter’s lower and upper bound. For each estimated parameter that was not a recruitment deviation the lower bound was set at 0.25 times the true parameter value and the upper bound was set at four times the true parameter value. The jitter for the recruitment deviations was the product of a standard normal random deviate times 0.4.
The algorithm has an initial loop that evaluates five sets of jittered stock synthesis results and determines if these runs produced valid likelihood values (rather than ones that were “not a number”) and had invertible Hessian matrices. For those runs that are acceptable, the algorithm saves the results from the run that produces the smallest negative log-likelihood (NLL). At least one (and up to five) additional jittered stock synthesis runs are then conducted, evaluated for acceptability (valid NLL and an invertible Hessian matrix) and compared with the best results obtained from the first set of five. The algorithm returns either with (a) the best results from the first set of five (if the NLL from the sixth run is greater than the best NLL from the first set of five) or with results from the next run that produces a valid NLL and invertible Hessian and has a smaller NLL than the best NLL from the first set of five, or (b) a flag indicating non-convergence.
For the vast majority of the 5400 experimental treatments the algorithm returned after doing six stock synthesis runs, but a small fraction of the treatments required one or more additional runs. On average, each treatment used 6.008 runs and 5.96 of the runs produced valid NLL values and invertible Hessians. For 3570 of the 5400 treatments (66.1%) the different sets of parameters produced by each run resulted in all the NLL values being equal to the minimum NLL value found across the runs for that treatment. However, there were multiple local minimum NLL values for about one third of the treatments. There were two distinct NLL values for 1377 treatments (25.5%), three distinct values for 409 treatments (7.6%), four distinct values for 40 treatments (0.7%), and five distinct values for 4 treatments (0.1%).
Figure Legends
Figure S.1. Illustration of regional exploitation and environmental patterns that drive the spatial distribution of the recruits. The solid lines indicate the patterns for the expected values. Each replicate includes random variation. The three graphs on the left represent the three levels of the environmental pattern (EnvirP) influencing the expected fractional recruitment to regionone over time: (1) Constant(EC), where there is no systematic change in the environmental pattern over time, and therefore no systematic change in the fractional recruitment to the regions; (2)Gradual (EG), where there is a sinusoidal environmental pattern driving the distribution of recruits; and (3) Abrupt(EA), where the expected fractional recruitment to regionone suddenly drops to a lower level in year 12 of the time series. The three graphs on the right represent the three exploitation histories (ExpRate) experienced by the two regions: (1) Constant (FC), where there are no time trends in expected fishing mortality in each region; (2) Mixed (FM), where the expected fishing mortality rate in one region increases over time, while the expected fishing rate in the other region decreases over time at the same rate; and (3) Increasing (FI), where the expected fishing mortality rate increases in both regions.The two shades of lines indicate the two regions.
Figure S.2. Summary of bias (median relative error, MRE) for estimates of SSB0, SSBCurrent, and Depletion from each of the nine operating model scenarios and six estimation model configurations. For each replicate the relative error of Depletion is related algebraically to the ratio of the relative error of SSBCurrent to the MRE of SSB0. When the MRE of SSB0 is less than the RE of SSBCurrent the RE of Depletion is positive. The distinct zigzag patterns for each type of relative error within each panel illustrate clear differences between MCs with and without survey data. Configurations that use survey data are less biased than those without survey data. When comparing two-region MCs (MC2.S.nE)with corresponding one-region MCs (MC1.S)using the correct spatial structure does not always improve assessment performance. The changes in bias are also inconstant across operating model scenarios. The factors EnvInd, ExpRate, and EnvirP produce complicated patterns indicating interactions among these factors. The solid horizontal line at a median value of zero indicates no bias in results.
Figure S.3. Summary of accuracy (median absolute RE, MARE) for estimates of SSB0, SSBCurrent and Depletion from each of the nine operating model scenarios and six estimation model configurations. Configurations that use survey data are more accurate than those without survey data. Accuracy of SSB0is improved by the environmental index for all scenarios except FC:EAand FC:EG. The largest effect of ExpRateon the accuracy of SSBCurrentis seen when EnvirPis gradual (EG) and ExpRateswitches between constant (FC) and mixed (FM) under MC1.nS. When considering the overall effects of these factors on the accuracy of SSB0, ExpRatehas more of an effect than EnvirP.The solid horizontal line at a median value of zero indicates no bias in results.
Figure S.4. Pseudo-code for the algorithm used to evaluate whether the stock synthesis runs had converged on the maximum likelihood estimates.
# Initialize some constants
Min.Iter.to.Test = 5
Max.Iter.to.Test = 10
OK.Hessian.Tolerance = 0.1
Hessian.Status = Bad
Loop1.Best.Available = FALSE
Use.Loop1.Best = FALSE
Loop1.Best.Iter = 0
Best.Iteration = 0
# NLL is the abbreviation for negative log-likelihood
Loop1.Best.NLL = 999999999
Best.NLL = 999999999
# begin first loop
forIter in 1 to Min.Iter.to.Test
do a jittered (10%) Stock Synthesis run
get SS3 output variables (e.g., NLL, SB0, R0, SB25, DEPL, SBMSY, MSY)
test for valid NLL and inverted Hessian and set OK.Hessian.YN
ifOK.Hession.YN = "Y"
if NLL < Best.NLL
save the SS3 output files as best-files
Loop1.Best.Available = TRUE
Loop1.Best.NLL = NLL
Loop1.Best.Iteration = Iter
Best.NLL = NLL
Best.Iteration = Iter
end if
end if
save the results( Iter, NLL, OK.Hessian.YN, SB0, R0, SB25, DEPL, SBMSY, MSY )
end for # first loop
# begin second and final loop
Init.Iter = Min.Iter.to.Test + 1
forIter in Init.Iter to Max.Iter.to.Test
do a jittered (10%) Stock Synthesis run
get SS3 output variables
if valid NLL # try again if invalid NLL
if NLL < Best.NLL
Best.NLL = NLL
Best.Iteration = Iter
end if
test for inverted Hessian and set OK.Hessian.YN
ifOK.Hession.YN = "Y"
# case 1: Likelihood within Tolerance of smallest NLL
case abs( NLL - Best.NLL ) = OK.Hessian.Tolerance
Hess.Msg = "OK Hessian and best likelihood or within tolerance of best likelihood"
OK.Hessian.YN[iter] = "Y<"
Hess.status <- Hess.OK # set the return value
break
# case 2: use saved Loop 1 Stock Synthesis results that had OK Hessian and best NLL
case Loop1.Best.Available
Use.Loop1.Best = TRUE
Hess.Msg ="OK Hessian but gave results from Loop.1 as the best likelihood" )
Hess.status <- Hess.OK # set the return value
break
# case 3: not smallest NLL, try again
case NLL > Best.NLL
Hess.Msg = "OK Hessian but not best likelihood"
end case
end if # OK.Hession.YN = "Y"
end if # valid NLL
save the results( Iter, NLL, OK.Hessian.YN,SB0, R0, SB25, DEPL, SBMSY, MSY )
end for # second loop
# special case 1: bad Hessian after last iteration but OK results available from Loop 1
ifHess.status is Bad and Loop1.Best.Available
Hess.msg = "Bad Hessian but gave results from Loop.1 best likelihood from iteration ", Loop1.best.iter
OK.Hess.YN[Loop1.best.iter] <- "Y<"
replace current SS3 output files with the best-files saved from Loop 1
Hess.status <- Hess.OK # set the return value
end if
# special case 2: OK Hessian at break from Loop 2 but better OK results available from Loop 1
if Use.Loop1.Best
Hess.msg = "OK Hessian but gave results from Loop.1 best likelihood from iteration ", Loop1.best.iter
OK.Hess.YN[Loop1.best.iter] <- "Y<"
replace current SS3 output files with the best-files saved from Loop 1
end if
write the following
Hess.msg
"Stopped after iteration", iter
"Best.like =", best.like, "at iter =", best.iter )
save the results( Iter, NLL, OK.Hessian.YN,SB0, R0, SB25, DEPL, SBMSY, MSY )
return(Hess.status )
Table S.1. Biological and fishery parameters used in the operating model. Values were derived from past stock assessments for black rockfish. The same values are used for each of the two regions, except for the region-one initial fractional recruitment and the regional initial F values.
Parameter / ValueNatural mortality (M, y-1) / 0.16
Steepness (h) / 0.6
Average maximum length ( cm) / 44.2
Growth coefficient (k, y-1) / 0.2266
Weight-length relationship alpha(α, kg/cm3) / 1.68E-05
Weight-length relationship beta (β) / 3
Length at 50% maturity (cm) / 39.53
Maturity versus length slope coefficient (cm-1) / 0.4103
Number of age-0 fish at unfished equilibrium (R0) / 2000
Log-scale standard deviation (SD) for annual recruitment (σR) / 0.5
Log-scale SD for exploitation (σF) / 0.2
Log-scale SD for recruitment distribution (σRecD) / 0.2
Initial fractional recruitment for region one (RecFracinit) / 60%
Initial Ft for region one (y-1) / 0.11
Initial Ft for region two (y-1) / 0.055
Fishery Logistic Selectivity (Sa)
Age at 50% selection (, y) / 6
Widtha of the selection curve (b, y) / 1.5
Catchability coefficient (QF, y-1) / 0.0001
Survey Logistic Selectivity (SurvSa)
Age at 50% selection (y) / 4
Widtha of the selection curve (, y) / 1.5
Catchability coefficient (QS) / 0.2
aWidth of the interval between 5% and 95% selection.
Table S.2. Linear mixed effect models fit to the relative error values for both SSB0 and SSBCurrent. Factors EnvirP, ExpRate and MC were treated as fixed effects and Replicate was treated as a random effect. The table includes the degrees of freedom (df) for each model, the ΔAIC values used to select the best fitting model for each variable, and the negative log likelihood values (LogLik). The bolded values indicate the models selected for RE(SSB0) or RE(SSBCurrent).
SSB0 / SSBCurrentModel / df / ΔAIC / LogLik / ΔAIC / LogLik
EnvirP + ExpRate + MC + EnvirP*ExpRate + EnvirP*MC + ExpRate*MC + EnvirP*ExpRate*MC / 56 / 31.63 / 11075.63 / 20.80 / 5248.32
EnvirP + ExpRate + MC + EnvirP*MC + ExpRate*MC + EnvirP*ExpRate / 36 / 60.45 / 11071.22 / 43.55 / 5246.94
EnvirP+ExpRate+MC+EnvirP*MC+ ExpRate*MC / 32 / 28.82 / 11083.04 / 22.76 / 5253.34
EnvirP+ ExpRate+MC+EnvirP*MC+EnvirP*ExpRate / 26 / 587.54 / 10797.68 / 367.78 / 5074.82
EnvirP+ExpRate+MC+EnvirP*ExpRate+ExpRate*MC / 26 / 31.63 / 11075.63 / 20.80 / 5248.32
EnvirP + ExpRate + MC + ExpRate*MC / 22 / 0.00 / 11087.44 / 0.00 / 5254.71
EnvirP+ExpRate+MC+EnvirP*MC / 22 / 555.91 / 10809.49 / 346.99 / 5081.22
EnvirP+ExpRate+MC+EnvirP*ExpRate / 16 / 547.50 / 10807.69 / 338.93 / 5079.25
EnvirP+ExpRate+MC / 12 / 515.87 / 10819.51 / 318.14 / 5085.65
EnvirP / 5 / 2552.07 / 9794.41 / 1453.71 / 4510.86
ExpRate / 5 / 2452.05 / 9844.42 / 1378.66 / 4548.39
MC / 8 / 576.51 / 10785.19 / 357.78 / 5061.83