x 4 w ‘
. ., ,
:..;..:..m. we? x.
4.1.... 3%....uus.
f... x: g. 1... WC.
12‘ \cal
.{q‘fivi}! .

my... a has

$1.5...

3 ‘ “in. .

. a. .
Divan
4mm». «fa

3L»: .

355.1 3
(. 3|
‘1..-

«muwiew

huiflwto. .:

r. \!
a.

z... .1
, Java?
. ‘ V . . H, V fury?“
99...? . ,
‘4‘. . ‘ .
4

 

 

 

4.\.......
.. .51... C

Q0137

 

LIBRARY
Michigan State
University

 

 

 

This is to certify that the
dissertation entitled

MODEL SELECTION AND DATA WEIGHTING METHODS
FOR STATISTICAL CATCH-AT-AGE ANALYSIS:
APPLICATION TO 1836 TREATY WATER STOCK

Doctoral

 

ASSESSMENTS

presented by

Brian C. Linton

has been accepted towards fulfillment
of the requirements for the
degree in Department of Fisheries and
Wildlife and Program in
Ecology, Evolutionary Biology,
and Behavior

 

/\
LI W (9,. tom

 

:7 Major Professor’s Signature
’4‘“; I 1 9x007.
U .

 

Date

MSU is an affinnative-action, equal-opportunity employer

_.- -.n._

m.-.—.--.-.--.-.—.--.-

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATEDUE

DATEDUE

DAIEDUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/07 p‘IClRC/DateDuetindd-pfi

 

MODEL SELECTION AND DATA WEIGHTING METHODS FOR STATISTICAL
CATCH-AT-AGE ANALYSIS: APPLICATION TO 1836 TREATY WATER STOCK
ASSESSMENTS

By

Brian C. Linton

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Fisheries and Wildlife and Program in Ecology, Evolutionary Biology,
and Behavior

2007

ABSTRACT
MODEL SELECTION AND DATA WEIGHTIN G METHODS FOR STATISTICAL
CATCH-AT-AGE ANALYSIS: APPLICATION TO 1836 TREATY WATER STOCK
ASSESSMENTS
By
Brian C. Linton
Recommended harvest limits for lake trout Salvelinus namaycush and lake
Whitefish Coregonus clupeaformis stocks in the 1836 Treaty Waters of the Great Lakes
are based on statistical catch-at-age analysis (SCAA). The assessment models and
methods are similar to those used to assess fish stocks in many of the word’s major
fisheries. My objective was to evaluate these methods with an eye towards suggesting
improvements both for 1836 treaty waters and more generally. My results provide
general guidance to stock assessment scientists with regard to data weighting and
selecting among alternative assessment models. As a first step, I performed an analysis
of the Lake Huron lake Whitefish models’ sensitivity to changes in “known” inputs and
model structure, selected as examples of basic type of assessment used throughout treaty
waters for lake Whitefish and lake trout. All of the Lake Huron lake Whitefish models
were sensitive to changes in the methods used to estimate recruitment and time-varying
selectivity, as well as to changes in their objective functions, and this indicated that
further study of these aspects of the assessment methods was warranted.
Specifically with regard to the objective function, the assessment models were
sensitive to changes in pre-specified variances associated with process and observation

errors, which are used to weight the different data sources. This result is consistent with

concerns expressed more broadly in the literature. I evaluated alternative approaches for

estimating log catchability (process error) and log total catch (observation error) standard
deviations within SCAA using Monte Carlo simulations: an ad hoc approach that tunes
the model predicted log total catch standard deviation to match a prior value, and a
Bayesian approach using either strongly or weakly informative priors for log catchability
standard deviation. When process error variance is large relative to observation error
(likely for many fisheries), reliable estimates of log catchability and log total catch
standard deviations can be obtained in SCAA using a Bayesian approach with only a
weakly informative prior on log catchability standard deviation.

The sensitivity of the Lake Huron Whitefish models to the method used to model
time-varying selectivity is also consistent with indications in the broader literature that
SCAA assessments can be sensitive to misspecification of selectivity. I therefore
evaluated four approaches for modeling time-varying selectivity within SCAA using
Monte Carlo simulations: double logistic functions with one, two and all four of the
function parameters varying over time, as well as age-specific selectivity parameters that
all varied over time. None of these estimation methods out performed the others in all
cases. In addition, I compared model selection methods to identify good (i.e., accurately
matching the true fish population) estimation models. Degree of retrospectivity, the best
selection method, was based on a retrospective analysis of bias in model parameter
estimates as the data time series for estimation is sequentially shortened. I recommend
this method of model section when considering different time-varying selectivity

estimation approaches in SCAA.

ACKNOWLEDGMENTS

I would like to thank everyone who made this dissertation possible. My advisor
and graduate committee chairman, Dr. James Bence, who gave me the opportunity to
work on this project and taught me all about fish population dynamics and stock
assessment. The rest of my graduate committee including: Dr. Michael Jones, Dr. Daniel
Hayes, and Dr. Guilherme Rosa for their constructive criticism concerning my research
and this dissertation. The staff and students of the Quantitative Fisheries Center provided
both technical and moral support for this project. In particular, Dr. Michael Wilberg
helped me get my first Monte Carlo simulation study up and running. The Modeling
Subcommittee of the ‘1836 Treaty Waters provided the initial impetus for this work, as
well as providing technical knowledge on the lake trout and lake Whitefish stock
assessments for the 1836 Treaty Waters. The Michigan Department of Natural Resources
and US. Fish and Wildlife Service provided the funding for this project.

I would also like to thank my family and friends for their mental, emotional, and
spiritual support over the last five years. Most notably, I would like to thank my wife,
Jennifer Linton, who I had the good fortune to meet and marry during my time here. In
addition, I thank my parents, Clifford and Irene Linton, and my brothers, Eric and John
Linton, for supporting me in this venture, as in every undertaking I have attempted in life.
Special thanks also go to all my friends at Martin Luther Chapel and the Graduate
Intervarsity Christian Fellowship chapter. Lastly, and in the chief place, I thank God for

blessing me with the desire and the ability needed to pursue this course. Thank you all.

iv

TABLE OF CONTENTS

LIST OF TABLES ............................................................................................................. vi
LIST OF FIGURES ............................................................................................................ x
CHAPTER 1
SENSITIVITY ANALYSIS OF LAKE WHITEFISH STOCK ASSESSMENT MODELS
USED IN THE 1836 TREATY WATERS OF LAKE HURON ......................................... 1
Introduction .......................................................................................................................... 1
Methods ................................................................................................................................ 4
Stock Assessment Model ......................................................................................... 5
Projection Model ...................................................................................................... 9
Sensitivity Analysis ............................................................................................... 10
Results ................................................................................................................................ 15
Discussion .......................................................................................................................... 24
References .......................................................................................................................... 29
CHAPTER 2
EVALUATING METHODS FOR ESTIMATING PROCESS AND OBSERVATION
ERRORS IN STATISTICAL CATCH-AT- AGE ANALYSIS ......................................... 49
Introduction ........................................................................................................................ 49
Methods .............................................................................................................................. 53
Data Generating Model .......................................................................................... 54
Estimation Models. ................................................................................................ 57
Ad Hoc Estimation Model ...................................................................................... 58
Bayesian Estimation Models .................................................................................. 60
Estimation Model Evaluation ................................................................................ 62
Results ................................................................................................................................ 63
Discussion .......................................................................................................................... 66
References .......................................................................................................................... 73
CHAPTER 3
EVALUATING AND SELECTING METHODS FOR ESTIMATING TIME-VARYING
SELECTIVITY IN STATISTICAL CATCH-AT-AGE ANALYSIS ............................... 85
Introduction ........................................................................................................................ 85
Methods .............................................................................................................................. 88
Data Generating Model .......................................................................................... 89
Estimation Models ................................................................................................. 93
Model Selection Methods ...................................................................................... 99
Results .............................................................................................................................. 103
Estimation Models ............................................................................................... 104
Model Selection ................................................................................................... 106
Discussion ........................................................................................................................ 108
References ........................................................................................................................ 11 3

LIST OF TABLES

Table 1.1. Predicted values for fully selected trap-net mortality (TN F), fully selected
gill-net mortality (GN F), population biomass (lbs), SSBR for unfished
population, SSBR at reference mortality schedule, SSBR ratio, projected
TAC/HRG (lbs), and the negative log-likelihood values from the unmodified lake
Whitefish stock assessment models for the 1836 treaty-ceded waters of Lake
Huron in 2001 ........................................................................................................ 31

Table 1.2. Percent difference from baseline values for fully selected trap-net mortality
(TN F), fully selected gill-net mortality (GN F), population biomass, SSBR for
unfished population, SSBR at reference mortality schedule, SSBR ratio, and
projected TAC/HRG from the lake Whitefish stock assessment models for the
1836 treaty-ceded waters of Lake Huron in 2001, when data input values were
increased (+), decreased (-), and set to specific values. Some changes led the
models to fail to converge (fc) ............................................................................... 32

Table 1.3. Percent difference from baseline values for fully selected trap-net mortality
(TN F), fully selected gill-net mortality (GN F), population biomass, SSBR for
unfished population, SSBR at reference mortality schedule, SSBR ratio, and
projected TAC/HRG from the lake Whitefish stock assessment models for the
1836 treaty-ceded waters of Lake Huron in 2001, when data input values were
increased (+), decreased (-), and set to specific values. For maturity schedule, an
increase (+) means maturity values were shifted up to the next oldest age, while a
decrease (-) means maturity values were shifted down to the next youngest age.
Some changes led the models to fail to converge (fc) ........................................... 33

Table 1.4. Percent difference from baseline values for fully selected trap-net mortality
(TN F), fully selected gill-net mortality (GN F), population biomass, SSBR for
unfished population, SSBR at reference mortality schedule, SSBR ratio, and
projected TAC/HRG from the lake Whitefish stock assessment models for the
1836 treaty-ceded waters of Lake Huron in 2001, when data input values were
increased (+), decreased (—), and set to specific values .......................................... 35

Table 1.5. Percent difference from baseline values for fully selected trap-net mortality
(TN F), fully selected gill-net mortality (GN F), population biomass, SSBR for
unfished population, SSBR at reference mortality schedule, SSBR ratio, and
projected TAC/I-IRG from the lake Whitefish stock assessment models for the
1836 treaty-ceded waters of Lake Huron in 2001, when natural mortality
parameters were increased (+) and decreased (-) ................................................... 36

Table 1.6. Percent difference from baseline values for fully selected trap-net mortality
(TN F), fully selected gill-net mortality (GN F), population biomass, SSBR for
unfished population, SSBR at reference mortality schedule, SSBR ratio, and
projected TAC/HRG from the lake Whitefish stock assessment models for the

vi

1836 treaty-ceded waters of Lake Huron in 2001, when catchability parameters
were increased (+) and decreased (—) ..................................................................... 37

Table 1.7. Percent difference from baseline values for fully selected trap-net mortality

(TN F), fully selected gill-net mortality (GN F), population biomass, SSBR for
unfished population, SSBR at reference mortality schedule, SSBR ratio, and
projected TAC/HRG from the lake Whitefish stock assessment models for the
1836 treaty-ceded waters of Lake Huron in 2001, when data input values were
increased (+) and decreased (-) .............................................................................. 38

Table 1.8. Percent difference from baseline values for fully selected trap-net mortality

(TN F), fully selected gill-net mortality (GN F), population biomass, SSBR for
unfished population, SSBR at reference mortality schedule, SSBR ratio, and
projected TAC/HRG from the lake Whitefish stock assessment models for the
1836 treaty—ceded waters of Lake Huron in 2001, when recruitment parameters
were increased (+) and decreased (-). Some changes led the models to fail to
converge (fc) ......................................................................................................... 39

Table 1.9. Percent difference from baseline values for fully selected trap—net mortality

(TN F), fully selected gill-net mortality (GN F), population biomass, SSBR for
unfished population, SSBR at reference mortality schedule, SSBR ratio, and
projected TAC/HRG from the lake Whitefish stock assessment models for the
1836 treaty-ceded waters of Lake Huron in 2001, when selectivity parameters
were increased (+) and decreased (-). The selectivity fimction parameters were
the first inflection point (p1), first slope (p2), second inflection point (p3), and
second slope (p4). Some changes led the models to fail to converge (fc) ............ 40

Table 1.10. Percent difference from baseline values for fully selected trap-net mortality

(TN F), fully selected gill-net mortality (GN F), population biomass, SSBR for
unfished population, SSBR at reference mortality schedule, SSBR ratio, and
projected TAC/HRG from the lake Whitefish stock assessment models for the
1836 treaty-ceded waters of Lake Huron in 2001, when selectivity parameters
were increased (+) and decreased (-). The selectivity function parameters were
the first inflection point (p1 ), first slope (p2), second inflection point (p3), and
second slope (p4). Some changes led the models to fail to converge (fc) ............ 42

Table 1.11. Percent difference from baseline values for fully selected trap-net mortality

(TN F), fully selected gill-net mortality (GN F), population biomass, SSBR for
unfished population, SSBR at reference mortality schedule, SSBR ratio, and
projected TAC/HRG from the lake Whitefish stock assessment models for the
1836 treaty-ceded waters of Lake Huron in 2001, when likelihood emphasis
factors were increased (+) and decreased (-). Some changes led the models to fail
to converge (fc) ..................................................................................................... 43

Table 1.12. Percent difference from baseline values for fully selected trap-net mortality

(TN F), fully selected gill-net mortality (GN F), population biomass, SSBR for

vii

unfished population, SSBR at reference mortality schedule, SSBR ratio, and
projected TAC/HRG from the lake Whitefish stock assessment models for the
1836 treaty-ceded waters of Lake Huron in 2001, when model structure was
modified. Some changes led the models to fail to converge (fc) .......................... 45

Table 1.13. Predicted values for fully selected trap-net mortality (TN F), fully selected
gill-net mortality (GN F), population biomass (lbs), SSBR for unfished
population, SSBR at reference mortality schedule, SSBR ratio, projected
TAC/HRG (lbs), and the negative log-likelihood values from the lake Whitefish
stock assessment models for the 1836 treaty-ceded waters of Lake Huron in 2001,
when changes improved model fit. The selectivity function parameter p4 was the
second slope. Selectivity was abbreviated sel., and decrease was abbreviated

dec .......................................................................................................................... 47
Table 2.1. Symbols and descriptions of variables used in data generating and estimation

models .................................................................................................................... 75
Table 2.2. Data generating and estimation model equations ............................................ 77
Table 2.3. Posterior probability density equations for estimation models ........................ 78

Table 2.4. Values of variables used in data generating model to create simulated data
sets .......................................................................................................................... 79

Table 3.1. Symbols and descriptions of variables used in data generating and estimation

models .................................................................................................................. 1 16
Table 3.2. Data generating and estimation model equations .......................................... 120
Table 3.3. Posterior probability density equations for estimation models ...................... 121

Table 3.4. Values of quantities used in data generating model to create simulation data
sets ........................................................................................................................ 122

Table 3.5. Values used to define prior probability densities in estimation models ........ 123

Table 3.6. Median relative errors (MRE), median absolute relative errors (MARE), and
number of replicates (N) for estimates of final population biomass and
exploitation rate produced by the time-varying selectivity estimation models:
double logistic functions with one (DLl), two (DL2), and four (DL4) time-
varying parameters, and time-varying age-specific selectivity parameters
(ASP) ................................................................................................................... 124

Table 3.7. Median relative errors (MRE), median absolute relative errors (MARE), and

number of replicates (N) for estimates of final population biomass and
exploitation rate chosen by the model selection methods: root mean square error

viii

(RMSE), deviance information criterion (DIC), and degree of retrospectivity
(DR) .................................................................................................................... 125

ix

LIST OF FIGURES

Figure 1.1. 1836 treaty-ceded waters and lake Whitefish management units in lakes
Huron, Michigan and Superior .............................................................................. 48

Figure 2.1. Box plots representing relative error distributions for estimates of log total
catch standard deviation across different levels of catchability and total catch
variance. The bars represent median relative errors. The boxes, whiskers, and
circles represent 25th and 75th, 10th and 90th, and 5th and 95th percentiles of the
distributions, respectively ...................................................................................... 80

Figure 2.2. Differences in median absolute relative errors (MARE) between informative
Bayesian approach and ad hoc approach across different levels of catchability and
total catch variance. Symbols represent informative Bayesian approach MARE
values minus ad hoc approach MARE values ....................................................... 81

Figure 2.3. Differences in median absolute relative errors (MARE) between the objective
Bayesian approach and the ad hoc approach across different levels of catchability
and total catch variance. Symbols represent objective Bayesian approach MARE
values minus ad hoc approach MARE values ....................................................... 82

Figure 2.4. Box plots representing relative error distributions for estimates of log
catchability standard deviation across different levels of catchability and total
catch variance. The bars represent median relative errors. The boxes, whiskers,
and circles represent 25th and 75th, 10th and 90th, and 5th and 95th percentiles of
the distributions, respectively ................................................................................ 83

Figure 2.5. Box plots representing relative error distributions for estimates of total
abundance in the last year of analysis across different levels of catchability and
total catch variance. The bars represent median relative errors. The boxes,

whiskers, and circles represent 25th and 75th, 10th and 90th, and 5th and 95th
percentiles of the distributions, respectively .......................................................... 84

Figure 3.1. Box plots representing relative error distributions for estimates of population
biomass in the last year of analysis across different data generating models. The
data generating and estimation models include double logistic functions with one
(DLl), two (DL2), and four (DL4) time-varying parameters, and time-varying
age-specific selectivity parameters (ASP). The bars represent median relative
errors. The boxes, whiskers, and circles represent 25th and 75th, 10th and 90th,
and 5th and 95th percentiles of the distributions, respectively ............................ 126

Figure 3.2. Box plots representing relative error distributions for estimates of
exploitation rate in the last year of analysis across different data generating
models. The data generating and estimation models include double logistic
functions with one (DLl), two (DL2), and four (DL4) time-varying parameters,

and time-varying age-specific selectivity parameters (ASP). The bars represent
median relative errors. The boxes, whiskers, and circles represent 25th and 75th,
10th and 90th, and 5th and 95th percentiles of the distributions, respectively ....127

Figure 3.3. The percentage of model nms when the model selection methods chose the
best or nearly best estimation model based on estimates of final population
biomass. The model selection methods include root mean square error (RMSE),
deviance information criterion (DIC), and degree of retrospectivity (DR). The
best or nearly best estimation model(s) is defined as the model(s) producing A)
the lowest final population biomass relative error, B) within 5% of the lowest
final population biomass relative error, and C) within 10% of the lowest final
population biomass relative error ......................................................................... 129

Figure 3.4. The percentage of model runs when the model selection methods chose the
best or nearly best estimation model based on estimates of final exploitation rate.
The model selection methods include root mean square error (RMSE), deviance
information criterion (DIC), and degree of retrospectivity (DR). The best or
nearly best estimation model(s) is defined as the model(s) producing A) the
lowest final exploitation rate relative error, B) within 5% of the lowest final
exploitation rate relative error, and C) within 10% of the lowest final exploitation
rate relative error .................................................................................................. 131

Figure 3.5. Box plots representing relative error distributions for estimates of population
biomass in the last year of analysis chosen by model selection methods across
different data generating models. The data generating models include double
logistic functions with one time-varying parameter (DLl) and time-varying age-
specific selectivity parameters (ASP). The model selection methods include root
mean square error (RMSE), deviance information criterion (DIC), and degree of
retrospectivity (DR). The bars represent median relative errors. The boxes,
whiskers, and circles represent 25th and 75th, 10th and 90th, and 5th and 95th
percentiles of the distributions, respectively ........................................................ 132

Figure 3.6. Box plots representing relative error distributions for estimates of
exploitation rate in the last year of analysis chosen by model selection methods
across different data generating models. The data generating models include
double logistic functions with one time-varying parameter (DLl) and time-
varying age-specific selectivity parameters (ASP). The model selection methods
include root mean square error (RMSE), deviance information criterion (DIC),
and degree of retrospectivity (DR). The bars represent median relative errors.
The boxes, whiskers, and circles represent 25th and 75th, 10th and 90th, and 5th
and 95th percentiles of the distributions, respectively ......................................... 133

xi

CHAPTER 1
SENSITIVITY ANALYSIS OF LAKE WHITEFISH STOCK ASSESSMENT MODELS

USED IN THE 1836 TREATY WATERS OF LAKE HURON

Introduction

In 1836, Native American Bands in the region to become the state of Michigan
signed a treaty with the US. government which reserved their right to fish in the
Michigan waters of lakes Huron, Michigan, and Superior. These fishing rights were
reaffirmed by the US. federal courts in 1979. The federal district court later approved
the fishery regulations created by the Chippewa/Ottawa Treaty Fishery Management
Authority (COTFMA) in 1982, while mandating that total allowable catches (TACs) or
harvest regulating guidelines (HRGs) be established for important fish species in order to
prevent over-fishing. Federal, state, and tribal biologists worked together to estimate
TACs for lake Whitefish Coregonus clupeaformis during 1979-1982. During this period,
the stock assessment methods used in the treaty waters were evolving and constrained by
limited data. Where possible stock sizes were estimated by application of a simple age-
structured model. Although there was no formal harvest policy, TACs were generally set
near the estimated maximum sustainable yield if the stock size was near the associated
biomass and to lower values when stock sizes were lower (e. g., AHWG 1979).

The 1985 Consent Decree laid out a 15 year agreement between federal, state and
tribal agencies for the allocation of fishery harvest between the parties. The Technical
Fisheries Review Committee (TFRC) was created by the decree to assess stocks of

important fish species. As part of this mandate, the TF RC recommended TAC/HRGs for

lake Whitefish stocks within the ceded territory to federal, state and tribal governments.
Stock assessments produced for the TFRC were generally based on simple age-structured
models (Clark and Smith 1984). The 1985 decree did not specify a harvest policy, but
based on TTWG (1984) the TFRC adopted a policy to limit total mortality to specified
levels less than 70%.

The 2000 Consent Decree was a new 20 year agreement, which set guidelines for
the management of important fish species, as well as allocating fishery harvest. As part
of the new decree, the Technical Fisheries Committee (TF C) was formed, which serves
many of the same functions as did the TFRC under the previous decree. Also at this
time, COTFMA was reorganized as the Chippewa/Ottawa Resource Authority (CORA).
Unlike the previous decree, a reference mortality rate for lake Whitefish of 65% was
specified, which partially defines a harvest policy. New methods for conducting lake
Whitefish stock assessments and projecting TAC/HRGs were developed during the
negotiation period for the 2000 Consent Decree by an interagency modeling group. The
decree specifies that a newly formalized Modeling Subcommittee (MSC) of the TFC
should build upon the work of the interagency modeling group to continue the lake
Whitefish stock assessment program.

The new stock assessment methods employed statistical catch-at-age models,
which were created for each lake Whitefish stock by the interagency modeling group and
further developed by the MSC. These stock assessment models used catch-at-age and
effort data from the commercial fisheries to estimate population abundances, mortality
rates, fishery harvests, and other population parameters of interest. Estimated quantities

from the assessment models were used to project each stock’s abundance and mortality

rates into the future, and then TAC/HRGs were calculated from these projections and a
reference mortality rate.

The 2000 Consent Decree established requirements governing the calculation of
TAC/HRGs. The reference level of total annual mortality (65%) specified for lake
Whitefish plays a different role depending on whether the yield from a particular
management unit is allocated entirely to the tribes (tribal unit) or partially allocated to the
state (shared unit). For shared units, 65% total mortality is treated as an upper limit and
TACs are established so as to allocate the yield between the parties as specified in the
decree. State and tribal management agencies are responsible for separately
implementing management actions (e.g., limits on entry to the fishery, gear restrictions,
size limits, and trip limits) to constrain fishery yield at or below levels specified by
TACs. If state or tribal fishery harvest exceeded their TAC/HRG by 25% or more, either
in a single year or over the course of five years, then that party’s TAC in the following
year is reduced by the amount that the previous TAC was exceeded. For tribal units, 65%
total mortality is viewed as an upper target level, and management actions by the tribes
are intended to prevent this level from being exceeded on average.

One of the complications of applying a reference mortality rate to the results of
the new age-structured assessment models is that these models account for the fact that
fishing mortality varies with age. The MSC chose a conservative solution to the problem
for lake Whitefish by further defining the reference mortality rate. First, for the reference
mortality rate, the maximum total mortality across all ages was not to exceed the
specified value of 65% (for most units). In addition, the spawning stock biomass per

recruit (SSBR) at this mortality schedule was required to be at least 20% of the SSBR for

the unfished stock. If the SSBR was below the 20% threshold, then the maximum total
mortality was reduced until the resulting SSBR was at least 20% of the unfished SSBR.

Due to the rapid development and implementation of the stock assessment
models, not all of the approaches used in the models have been fully evaluated. For
example, there were numerous methods for modeling each of the biological processes
represented within the models from which the MSC analysts could select. Once a
particular method for modeling a process was chosen, reasonable parameter starting
values and bounds on what values those parameters could take also had to be selected by
the analysts. It was unknown how much these choices affected stock assessment results.
Therefore, my objective was to further evaluate the stock assessment models for lake
Whitefish in the 1836 treaty waters of Lake Huron, with a view toward suggesting
possible improvements. This objective linked to a broader goal for my work, to form the
basis for advice that is broadly applicable in the field of fishery stock assessment. As a
first step to achieve this objective, I performed a general analysis of the models’

sensitivity to changes in “known” inputs and model structure.

Methods

The 1836 treaty waters of Lake Huron were divided into five lake Whitefish
management units, each thought to contain a distinct lake Whitefish stock (Figure 1.1).
Separate stock assessment models were developed for each of the lake Whitefish
management units. When the models were originally developed, it was assumed that the

net movement of lake Whitefish between management units was nil.

Stock Assessment Model

Here I provide an overview of the stock assessment models’ general structure.
Ebener et a1. (2005) provides a detailed description of the models. All of the stock
assessment models consisted of two basic submodels, a population submodel and an
observation submodel. The population submodel described the population dynamics of
the stock in terms of abundance-at-age excluding the first year and the first age in

subsequent years:

‘Za
Na+l,y+l = Na,ye ,y ,

where Na,y was the number of fish in age a and year y and lay was the total

instantaneous mortality rate in age a and year y. Numbers-at-age in the first year were
estimated as a vector of relative population variation parameters (i.e. a vector of
deviations that must sum to zero). A population scaling parameter then converted these
deviations to numbers-at-age. Numbers of fish in the first age of each year also were
estimated as a series of scaled deviations using the same population scaling parameter,

but were penalized for deviating too greatly from a Ricker stock-recruitment function:

N _e‘.BGy—a0-l ’

aGy-(

00,)! = a0 -1)

where N was the number of fish in the first age a0 and year y, G y_(a0 _1) was the

a 0 , y
number of eggs produced aO-l years prior to year y, a was the productivity parameter,

and 6 was the density dependent parameter. The number of eggs was calculated within
the submodel, based on a constant weight-specific fecundity. The productivity and

density dependent parameters were estimated within the submodel. Numbers-at-age were

converted to biomass using observed mean weight-at-age data. Total mortality consisted

of four component parts:

Za,)’ = M +ML,a,y +FG,a,y + FT,a,y’
where M was the natural mortality rate, M L, a, y was the sea lamprey induced mortality rate
in age a and year y, F G, a, y was the gill net fishing mortality rate in age a and year y, and

F130,), was the trap net fishing mortality rate in age a and year y. Natural mortality was

assumed to be constant for all ages and years, and was estimated as a model parameter.
Pauly’s equation (Pauly 1980) was used to calculate an initial value for the natural
mortality parameter to provide a reasonable starting point for parameter estimation. Sea
lamprey mortality was calculated externally to the model based on observed sea lamprey
wounding rates. Fishing mortality was estimated by relaxing the assumptions of the fully

separable fishing mortality model and allowing gear selectivity to vary with time:
Fi,a,y = Si,a,yqiEi,y§i,y ’
where Si,a,y was the gear selectivity of age a fish in fishery i and year y, q; was the

catchability in fishery i, E; y was the observed fishing effort in fishery i and year y, and

(,3), was the deviation in fishing mortality from direct proportionality to observed fishing

effort in fishery i and year y. Selectivity was estimated with a double logistic function of
age, and one of the parameters of the function was allowed to change with time according
to a quadratic function. This allowed age-specific selectivity to change gradually over

time. An adjustment factor was applied to the observed gill net effort in order to account

for changes in the number of meshes deep that were set through time.

The observation submodel predicted catch-at-age for the gill net and trap net

fisheries. Catch-at-age was predicted using Baranov’s catch equation:

 

Fi,a,y —Z
Cifiay = Zay Na-y(l—e thy),

where Ci.a.y was the number of age a fish caught in fishery i during year y, and all of the

other parameters were estimated in the population submodel. Predicted catch-at-age was
converted to a total annual catch and a proportion of catch-at-age for each fishery. An
underreporting factor, representing the proportion of the actual catch that was reported,
was applied to the total catch in order to account for underreporting and discards in the
fisheries. The underreporting factor was obtained by comparing reported fishery landings
to actual sales.

The parameter values providing the best fit were found using Bayesian methods
(i.e., prior densities were assigned to all parameters). In particular, best fit parameter
estimates maximized the joint posterior density, and for numerical reasons this was done
by finding parameter values that minimized the weighted sum of the negative log
likelihoods and the negative log prior densities. Separate likelihood components were
calculated for gill net total catch, gill net proportion of catch-at—age, trap net total catch,
and trap net proportion of catch-at-age. Total annual catch was assumed to follow a
lognorrnal distribution, with the negative log likelihood (ignoring some additive

constants) given by:

where a,- was the standard deviation for log-scale observed total catch in fishery i, C ,3 y

A

was observed total numbers of fish caught in fishery i and year y, C was predicted

13y
total numbers of fish caught in fishery i and year y, and n was the total number of years
included in the model. Observed catch was reported as weight of fish harvested, which
was converted to numbers of fish using the observed mean weight of a harvested fish.

Proportion of catch-at-age was assumed to follow a multinomial distribution, with the

negative log likelihood (ignoring some additive constants) expressed as:
n m A
L<PJ= ‘ZNEM Elm... mam),
y=l a=l

where N E, ,3 y was the effective number of fish used to calculate the age composition in

fishery i and year y (Fournier and Archibald 1982), Pi.y.a was the observed proportion of

catch-at-age a in fishery i and year y, 13, a was the predicted proportion of catch-at-age

,y
a in fishery i and year y, n was the total number of years included in the model, and m
was the total number of ages included in the model. In addition to the likelihood
components, the joint posterior density included terms related to prior densities for the
model parameters. First, deviations of predicted recruitrnents from the Ricker stock-
recruitrnent function were assumed to follow a lognorrnal distribution. Second, deviation
of predicted natural mortality from the prior natural mortality value (i.e. the Pauly’s
equation value) was assumed to follow a lognorrnal distribution. Third, deviations in the
fishing mortality from direct proportionality to observed fishing effort were assumed to
follow a lognorrnal distribution. The log of all remaining model parameters were

assigned proper uniform prior densities, which follows common practice with the intent

of being weakly informative. Therefore, prior densities of the log of the remaining
parameters were constants for all parameter values.

Each likelihood component, the prior density for deviations between recruitment
and the stock-recruitment function predictions, the prior density for natural mortality, and
the prior density for deviations in the fishing mortality from direct proportionality to
observed fishing effort were weighted by an emphasis factor as described by Methot
(1990). If all likelihood components, prior densities, and their associated standard
deviations or effective sample sizes were correctly specified, then the emphasis factors
should all be 1.0. If there was a misspecification in the objective function, then the
emphasis factors provide a simple way for analysts to adjust how closely the model

attempts to fit observed and predicted data for each likelihood component.

Projection Model

Recommended yields for a reference (sometimes called target) mortality rate were
then calculated using stock assessment model output in a projection model. The stock
assessment model output included estimated numbers-at-age, estimated total mortality,
estimated natural mortality, and assumed sea lamprey mortality, all from the last year of
the model, as well as, estimated trap net and gill net mortality rates that were averaged
over the last three years of the model, and estimated average recruitment (over the last ten
years) of the model. Along with the stock assessment model output, observed weight-at-
age in the fisheries, observed mean proportion of females in the population, observed
maturity schedules represented as year and age-specific proportions, and observed time of
year of spawning represented as a proportion of the year were also used in the projection

model for SSBR calculations.

The projection model took the abundance-at-age estimates from the beginning of
the last year of the stock assessment model, projected abundance to the beginning of the
year for which recommended yields were desired, then projected yields for the trap net
and gill net fisheries. Trap net and gill net fishery multiplier parameters were used to
adjust age-specific fishing mortality rates by the same proportion for each age. The
values of the two multipliers were set so as to achieve the reference mortality rate, while
maintaining a desired allocation between trap net and gill net yield. There were two steps
to determining the appropriate values for the multipliers, which corresponded to how the
reference mortality was defined. First, the multipliers were adjusted so that the
maximum total annual mortality for any age did not exceed the reference (typically 65%).
Second, the ratio of SSBR at this mortality schedule to SSBR without fishing was
calculated (hereafter the SSBR ratio). If this ratio was less than 0.2, then the multipliers

were decreased until the SSBR ratio equaled 0.2.

Sensitivity Analysis

Sensitivity analysis quantifies the effect of changes made to a model’s input
values and underlying assumptions on the model’s output (Morgan and Henrion 1990).
My sensitivity analysis tested changes to the stock assessment models’ input quantities
and model structure (i.e., underlying model assumptions). Changes to observed input
data represented possible changes in data collection (e.g., collecting more or less data),
while changes to input values based on expert judgment (e.g., parameter starting values)
represented a changes made by the analyst during the model fitting process. Changes in

model structure were based on alternative modeling procedures suggested in the

10

literature. The WFH—03 management unit model was not included in the following
analysis, due to a lack of convergence to a satisfactory solution.

The stock assessment models were tested for their sensitivity to changes in input
values. The observed mean weight-at-age of harvested fish was varied for all ages at
once by i10% of the original values. The year- and age-specific maturity schedule was
varied by reassigning maturity values from each age to the next oldest age (e. g., maturity
values for age 4 fish became the maturity values for age 5 fish). Then the first age was
given a maturity value of zero. Similarly, the maturity schedule was varied by
reassigning maturity values from each age to the next youngest age (e.g., maturity values
for age 4 fish became the maturity values for age 3 fish), and setting the maturity in the
last age equal to 1.00. Fecundity was adjusted by making it a linear function of average
weight-at-age at the time of spawning. The gill net adjustment factors for number of
meshes deep set through time were set equal to 1.00 to test the overall effect of the
adjustments. The gill net adjustment factors also were varied using the following
formula, which assumed the trend in the factors over time was alternatively more and less

extreme than originally thought:
xy = 7c +chO,y -)—C),
where xy was the new adjustment factor in year y, 7: was the average of the original

adjustment factors across all years, x0,y was the original adjustment factor in year y, and

scalar c alternatively equaled 0.8 to represent a less extreme trend and 1.2 to represent a
more extreme trend. Adjustment factors were included in the original models to account
for underreporting in each year of the fisheries. The underreporting factors were set to

1.00 for one fishery at a time to test the overall effect of the adjustments. The

11

underreporting factors were increased and decreased by a value of 0.2 for one fishery at a
time. The proportion of females in the population was set to 0.5. The proportion of
females also was increased and decreased by a value of 0.2. The time of year of
spawning was increased and decreased by a value of 0.2. Bounds for each model
parameter, which limited the range of values a given parameter could take, were
increased one at a time by decreasing the lower bound by 20% of the original value and
increasing the upper bound by 20% of the original value. Bounds for each model
parameter were decreased one at a time by increasing the lower bound by 20% of the
original value and decreasing the upper bound by 20% of the original value. Starting
values for each model parameter were increased and decreased one at a time by 20% of
the original values. Natural mortality was altered by fixing the parameter to the starting
value and by increasing and decreasing the starting value by 20% of the original starting
value.

The stock assessment models were tested for their sensitivity to changes in model
structure. Recruitment in each year was estimated as a free parameter without any
penalty for deviating from stock-recruitment model predictions. Also, a Beverton-Holt
stock-recruitment function, rather than a Ricker stock-recruitment function, was used to

predict recruitment (Beverton and Holt 1957):

N _ aGy—ao—l
“W _ 1+ [30

 

y-ao—l.

Rather than using deviations between observed and predicted numbers of fish caught in
the objective function, deviations between observed and predicted biomass of fish caught
were used. The predicted numbers of fish caught were converted to mass of fish caught,

using the mass—at-age of a harvested fish, comparing them assuming a lognormal

12

distribution. The likelihood component emphasis factors were doubled and halved one at
a time. Gamma likelihood components were substituted for all lognormal likelihood
components, keeping the same coefficient of variation:

L(C,' ) = —¢i i In[ (Shy ] - (gay :I ,

y=1 13y Chy

 

 

where ¢,~ was the inverse of the squared coefficient of variation for observed harvest in

fishery i (Cadigan and Myers 2001) and the other variables were the same as in the
lognormal likelihood component. Dirichlet likelihood components were substituted for
all multinomial likelihood components, with fixed parameters setting the effective sample

size equal 100:
n m A m A m A
L(P.-)= Z 1nr[2r,-P.-,y,a]— Zlnrlrle,y,a)+ Xvi-Pry. —1)1nP.-,y,. .
y=l a=l a=l 0:]

where 7,- represented the effective sample size for fishery i, F was the gamma function,

and the other variables were the same as in the multinomial likelihood component.

Each stock assessment model was rerun for each of the changes tested. In order
to better specify the standard deviation around the stock recruitment relationship, an
initial recruitment standard deviation was input into the model. The standard deviation of
predicted recruitment was then calculated at the conclusion of model fitting. The
predicted recruitment standard deviation then replaced the former input standard
deviation, and the model was renm leading to a new predicted recruitment standard
deviation. This process was repeated 50 times with the goal of getting the ratio between
input recruitment standard deviation and predicted recruitment standard variation as close

to unity as possible. After the 50 runs, the model was considered to have converged to a

13

satisfactory solution if: 1.) the ratio of recruitment standard deviations was between 0.98

and 1.02; 2.) the maximum gradient component, which measures the maximum amount

of change in parameter estimates during model fitting, was less than 1 x 10-2; and 3.) the

Hessian matrix, which is used to calculate standard deviations for the parameter
estimates, was positive definite.

The sensitivity of the stock assessment models to change was monitored by
tracking several of the models’ output quantities. The output quantities of interest
included: the estimated fully selected gill net and trap net fishing mortality rates averaged
for the last three years of the assessment, estimated population biomass averaged for the
last three years of the assessment, estimated SSBR of the unfished population, predicted
SSBR at reference mortality levels, estimated SSBR ratio, and the estimated yield
calculated for reference mortality rates for the projected population. Model sensitivity
was calculated as the percent difference of the test quantity of the adjusted model from
the baseline value of the test quantity of the original model (Table 1.1):

a-eo
0

x100,

 

D%=

where D% was the percent difference, 60 was the baseline value of the test quantity,

andt9‘ was the value of the test quantity from the adjusted model. I considered a model to

be sensitive to a change if that change produced a 10% or greater change in one of the

output quantities.

14

Results

All of the stock assessment models were sensitive to changes in the input values.
As expected, increasing observed mean weight-at-age of a harvested fish led to an
increase in the projected TAC/HRG in all of the projection models (Table 1.2).
Likewise, decreasing mean weight-at-age led to a decrease in the projected TAC/HRG in
all of the projection models. These effects upon projected TAC/HRGs were greater for
the gill net fishery in WFH-Ol and WFH-04, and were greater for the trap net fishery in
WFH-02 and WFH-OS. Changes in mean weight-at-age of a harvested fish in the gill net
fishery had no effect upon the projected TAC/HRG in the WFH-OS model, due to the
small size of the gill net fishery in that management unit.

Surprisingly, setting gill net effort adjustment factors for number of meshes deep
set through time to one, increasing gill net effort adjustment factors, and decreasing gill
net effort adjustment factors increased the projected TAC/HRG by 34.7% and changed
the remaining test quantities to a lesser degree (0. 146%), except for SSBR of the
unfished population which was unaffected, in the WFH-02 model (Table 1.2). All of
these changes to gill net effort adjustment factors had slight effects (OJ-3.3%) on all of
the test quantities, except for SSBR of the unfished population, in all of the other models;
though no clear patterns were apparent. The WFH-Ol model failed to converge when gill
net effort adjustment factors were set equal to one.

As anticipated, shifting the maturity schedule later by one age led to substantial
decreases (21 .2-48.2%) in the SSBR of the unfished population and SSBR at the
reference mortality schedule, with a greater decrease in SSBR at the target schedule,

because the fish were maturing later after more mortality had occurred, and mortality was

15

higher for the reference schedule (Table 1.3). The greater decrease in SSBR at reference
mortality schedule led to a decrease in the SSBR ratio. Likewise, shifting the maturity
schedule earlier by one age led to substantial increases (13.6-55.4%) in the SSBR of the
unfished population, SSBR at the reference mortality schedule, and the SSBR ratio in all
of the models, due to the resulting increase in spawning biomass. Unexpectedly, shifting
the maturity schedule later by one age increased the projected TAC/HRG in WFH-02 and
WFH-04 by 38.5% and 14.4% respectively. Changes in the maturity schedule also had
modest influence on fully selected gill net and trap net mortality, biomass, and projected
TAC/HRG (01-59%) in all of the models. There was some influence because maturity
schedule values are used to calculate the number of eggs produced for the stock-
recruitrnent function, and this affects the objective function. The WFH-02 model failed
to converge when the maturity schedule was shifted earlier by one age.

Setting the average proportion of females in the population equal to 0.5 led to an
increase (21.8-31.6%) in SSBR of the unfished population and SSBR at the reference
mortality schedule in all of the models, except the WFH-04 model which failed to
converge (Table 1.3). As expected, increasing the proportion of females led to an
increase (47 .9-52.6%) in SSBR of the unfished population and SSBR at the reference
mortality schedule in all of the models, because the spawning stock was considered to be
the mature females within the population. Decreasing the proportion of females led to a
decrease (48.9-52.6%) in SSBR of the unfished population and SSBR at the reference
mortality schedule in all of the models, due to the resulting decrease in spawning stock.

All of the adjustments made to the proportion of females led to slight changes (OJ-5.4%)

l6

in the fully selected fishing mortalities, biomass, and SSBR ratio; and larger changes
(23.8-3 7.2%) in the projected TAC/HRG for the WFH-02 model.

As expected, increasing trap net fishery underreporting adjustment factors led to
increases in fully selected trap net mortality (12.8-18.1%) in order to account for the
increased trap net harvest, except for a 1.3% decrease in trap net mortality in the WFH-OS
model, and decreases in fully selected gill net mortality (7.6-22.3%; Table 1.4).
Decreasing trap net fishery underreporting adjustment factors led to decreases in fully
selected trap net mortality (10.7-13.4%) due to the lower trap net harvest, except for a
1.1% increase in trap net mortality in the WFH-OS model, and increases in fully selected
gill net mortality (5.5-22.6%). Likewise, increasing gill net fishery underreporting
adjustment factors led to increases in fully selected gill net mortality (6.8-29.8%) and
decreases in fully selected trap net mortality (13.4-16.7%) due to increased gill net
harvest, except for a 1.8% increase in trap net mortality in the WFH-OS model.
Decreasing gill net fishery underreporting adjustment factors led to decreases in fully
selected gill net mortality (6.1-18.7%) and increases in fully selected trap net mortality
(103-14.5%) due to decreased gill net harvest, except for a 1.2% decrease in trap net
mortality in the WFH-OS model. The small, but unforeseen, changes in firlly selected
fishing mortality rates (< 2%) in the WFH-OS model appeared to be due to the small gill
net fishery, which effectively makes WFH-OS a one (trap net) fishery system. It appears
the WFH-OS assessment model accounted for adjustments in observed trap net harvest by
making large changes to the biomass and small changes to fishing mortality. Likewise,
changes in gill net harvest led to only small adjustments of the biomass and trap net

fishing mortality because of the small size of the fishery. Changes in the fishery

17

underreporting adjustment factors also affected biomass (0.5-28.5%), SSBR at the target
mortality schedule (0-4.4%), and SSBR ratio (0-4.4%), though no patterns were apparent.
Changes in the fishery underreporting adjustment factors had no effect on the SSBR of
the unfished population.

As anticipated, increasing the time of year of spawning led to a decrease (2.7-
12.4%) in SSBR of the unfished population, SSBRat the reference mortality schedule,
and SSBR ratio because fewer fish survived to spawn later in the year (Table 1.3).
Increasing the time of year of spawning also led to an increase (0.3-36.2%) in the
projected TAC/HRG for all of the models because the spawning stock was exposed to the
fisheries for a longer period of time before spawning. Decreasing the time of spawning
led to an increase (2.8-14.5%) in SSBR of the unfished population, SSBR at the reference
mortality schedule, and SSBR ratio because more fish would survive to spawn earlier in
the year. Decreasing the time of spawning also led to a decrease (0.4-1.1%) in the
projected TAC/HRG for all of the models because the spawning stock was exposed to the
fisheries for a shorter period of time, except in the WFH-02 model which had an
unexpected increase in the TAC/I-IRG of 33.0%. Adjustments to the time of spawning
led to slight changes (0-5.1%) in the fully selected fishing mortalities and biomass with
no clear pattern in all of the models. These slight changes appeared because time of
spawning is used to calculate the number of eggs produced for the stock-recruitment
function, which influenced the objective function.

Both increasing and decreasing the parameter bounds for natural mortality led the
WFH-O2 model to converge to the same solution, different from the original one, where

fully selected trap net mortality decreased by 0.2%, fully selected gill net mortality

18

increased by 1.2%, biomass increased by 4.6%, SSBR of the unfished population
remained unchanged, SSBR and SSBR ratio decreased by 0.1%, and projected
TAC/l-IRG increased by 34.7% (Table 1.5). This alternate solution was very similar to
the one reached by the model for this unit when changes were made to the gill net effort
adjustment factors. Increasing and decreasing the natural mortality parameter’s starting
value for this unit also led to the same solution described above. Decreasing natural
mortality’s starting value led to 0.1% changes in fully selected trap net mortality,
biomass, SSBR and SSBR ratio, and projected TAC/I-IRG in the WFH-04 model. The
WFH-Ol and WFH-05 models were unaffected by changes to natural mortality.

Surprisingly, decreasing trap net catchability bounds and increasing and
decreasing gill net catchability bounds in the WFH-02 model led to the same alternate
solution described above, where projected TAC/HRG increases by 34.7% while all the
other test quantities, except for SSBR of the unfished population, changed from 01-46%
(Table 1.6). Increasing and decreasing the trap net catchability starting value and
increasing the gill net starting value again led to the same alternate solution for WFH-02
with the 34.7% increase in projected TAC/HRG. None of the other models showed any
sensitivity to changes in the catchability parameters.

Increasing and decreasing the population scaling parameter’s bounds, and
decreasing the population scaling parameter’s starting value led to the state with the
34.7% increase in projected TAC/HRG in the WFH-02 model (Table 1.7). Decreasing
the relative population variation parameters’ bounds led to a 31.9% increase in firlly
selected gill net mortality, a 25.2% increase in the projected TAC/HRG, and smaller

changes (22-80%) in fully selected trap net mortality, biomass, and SSBR and SSBR

19

ratio in the WFH-02 model. The other models were unaffected by changes to the
population scaling parameter and relative population variation parameters.

Increasing the bounds of the Ricker stock-recruitment function’s productivity
parameter led to a 37.2% increase in projected TAC/HRG and smaller changes (0.4-
5.4%) in the fully selected fishing mortalities and biomass in the WFH-02 model (Table
1.8). Decreasing the bounds of the Ricker function’s productivity parameter led to
changes (0.2-10.8%) in all of the test quantities, except for SSBR for the unfished
population, for the WP H-01 and WFH-OS models. The WFH-02 model failed to
converge when both the Ricker function’s productivity parameter’s bounds and starting
value were decreased. Increasing the bounds of the Ricker function’s density dependence
parameter, increasing the starting value of the productivity parameter, and increasing and
decreasing the starting value of the density dependence parameter led to the state where
the projected TAC/HRG increases by 34.7% in the WFH-O2 model. Increasing the
bounds of the Ricker function’s density dependence parameter led to a 20.6-43.4%
increase in the projected TAC/HRG and smaller changes (0.3-15.2%) in the remaining
test quantities, except for SSBR of the unfished population, for the WFH-Ol and WFH-04
models.

Increasing the bounds of the gill net selectivity function’s first inflection point and
decreasing the bounds of the gill net selectivity function’s first slope parameter in the
WFH-02 model led to the same state noted earlier with the 34.7% increase in the
projected TAC/HRG (Table 1.9). Decreasing the bounds of the gill net selectivity
function’s first inflection point led to changes (3.0-15.4%) in all of the test quantities,

except SSBR of the unfished population, in the WFH-02 and WFH-04 models. The

20

WFH-02 model failed to converge when the bounds on the gill net selectivity frmction’s
second inflection point were widened. Decreasing the bounds of the trap net selectivity
function’s first inflection point led to changes (5.3-18.1%) in all of the other test
quantities, except for SSBR of the unfished population, in the WFH-04 model (Table
1.10). The WFH-Ol model failed to converge when the starting value for the gill net
selectivity function’s second slope parameter was increased. Increasing and decreasing
the starting values for the gill net selectivity function’s first and second inflection points,
decreasing the starting values for the gill net selectivity function’s first and second slope
parameters, decreasing the starting value for the trap net selectivity function’s first
inflection point, increasing and decreasing the starting value for the trap net selectivity
function’s second inflection point, and decreasing the starting value for the trap net
selectivity function’s second slope parameter led to the alternate state with the 34.7%
increase in the projected TAC/HRG in the WFH-02 model. The WFH-02 model failed to
converge when the starting value for the trap net selectivity function’s first inflection
point was increased. The WFH-04 model failed to converge when the starting value for
the gill net selectivity function’s first inflection point was decreased. Adjustments to the
starting values for the gill net selectivity function’s parameters led to 0.2-131.0% changes
in fully selected trap net mortality, 1.4-1,277.2% changes in fully selected gill net
mortality, 1.6-50.9% changes in biomass, 0.7-14.4% changes in SSBR and SSBR ratio,
2.4—87.7% changes in projected TAC/HRG, and no change to SSBR for the unfished
population for the WFH-05 model. Increasing the starting value for the trap net
selectivity function’s second inflection point, and increasing and decreasing the starting

value for the trap net selectivity function’s second slope parameter led to a 0.1% change

21

in firlly selected trap net mortality for the WFH-Ol model. Decreasing the starting value
for the trap net selectivity function’s second inflection point led to changes (OJ-0.4%) in
all of the test quantities, except for SSBR of the unfished population, for the WFH-Ol
model.

Increasing and decreasing the likelihood emphasis factor for natural mortality led
to the alternate solution with a 34.7% increase in the projected TAC/HRG for the WFH-
02 model (Table 1.11). The WFH-Ol model failed to converge when the trap net catch
and age composition emphasis factors were increased. The WFH-O4 model failed to
converge when the trap net catch emphasis factor was increased, and when the trap net
and gill net age composition emphasis factors were decreased. All the remaining
adjustments to the likelihood emphasis factors led to positive and negative changes (0.1-
62.0%) that showed no pattern in all of the test quantities, except SSBR of the unfished
population, for all of the models.

All of the stock assessment models also were sensitive to changes in model
structure. Holding natural mortality constant at its starting value in the WFH-02 model
led to the state with the 34.7% increase in projected TAC/HRG (Table 1.12). Modeling
fecundity as a linear function of weight led to changes (0-38.8%) in all of the test
quantities, except for SSBR of the unfished population, for the WFH-Ol , WFH-02, and
WFH-OS models, because fecundity was used to calculate the number of eggs produced
(stock size) for the stock-recruitment function (Table 1.6). The WFH-04 model failed to
converge when fecundity was modeled as a linear function of weight.

Estimating each year’s recruitment as a free parameter led to changes (0. 1 -39.8%)

in all of the test quantities, except for SSBR of the unfished population, for the WFH-02

22

and WFH-OS (Table 1.12). The WFH-Ol and WFH-04 models failed to converge when
recruitment was estimated as free parameters. Estimating recruitment using a Beverton-
Holt stock-recruitment model led to (0.3-54.3%) changes in all of the test quantities,
except for SSBR of the unfished population, in all of the models.

Fitting mass, instead of numbers, of fish caught in the objective firnction led to
changes (3.7-39.6%) in all of the test quantities, except for SSBR of the unfished
population, for the WFH-02, WFH-04, and WFH-OS models (Table 1.12). The WFH-Ol
model failed to converge when the mass of fish caught was used in the objective function.

The use of the gamma likelihood function in place of the lognormal likelihood
function led to small changes (01-07%) in all of the test quantities, except for SSBR of
the unfished population, in the WFH-Ol, WFH-04, and WFH-OS models (Table 1.12).
The WFH-02 model failed to converge when the gamma likelihood firnction was used.
The use of the Dirchlet likelihood function in place of the multinomial likelihood
function led to changes (0-25.9%) in all of the test quantities, except for SSBR of the
unfished population, for all of the models.

Most of the adjustments made to the models led to negative log-likelihood values
that were the same as, or higher than, the original likelihood values, which means that the
model fit was not improved. In particular, the alternate solution often arrived at by the
WFH-02 model had a higher likelihood value (4,340.5) than the original model (4,337.6).
There were, however, several changes that led to a decrease in the negative log-likelihood
value, which means that the changes produced parameter estimates that fit the data better
than the original parameter estimates. In particular a better fit was obtained after

decreasing the bounds of the Ricker recruitment function’s density dependence parameter

23

in the WFH-Ol and WFH-04 models and after decreasing the starting value of the gill net
selectivity function’s second slope parameter in the WFH-OS model (Table 1.13). These
instances of better model fit could be due to random chance given the large number of
model changes explored. Likelihood values could not be directly compared to determine
better model fit in cases where the model structure was changed or when the likelihood

emphasis factors were adjusted, because these changes altered the objective function.

Discussion

I performed a simple sensitivity analysis of the stock assessment models for lake
Whitefish in the 1836 treaty waters of Lake Huron to changes in input quantities and
model structure. The changes I tested could be divided into two alternate categories that
affect the way in which the results are interpreted. First, changes to the observed data
and model structure led to changes in the objective function (negative log-likelihood) and
thus altered the optimal solution (i.e., the best-fit parameter estimates) from the optimal
solution of the baseline model. In this case, changes in the output quantities represent the
model seeking the new optimal solution. Second, changes to the parameter starting
values and parameter bounds did not alter the optimal solution from the optimal solution
of the baseline model. In this case, changes in the output quantities mean that the model
has become trapped at a local minimum for the objective function or found the true
global minimum for the objective function depending upon whether the likelihood value
is greater than or less than, respectively, the baseline model’s likelihood value.

A simple sensitivity analysis, like the one conducted here, can be useful for
identifying models that are unstable and highly sensitive to change. The WFH-02 stock

assessment model appears to be such a sensitive model. Thirty-five of the 111 changes

24

tested led the model to converge to an alternate solution, which provided a poorer fit
between observed and predicted values than the original model. The alternative solution
was similar to the original model’s solution except for a large increase in the projected
TAC/HRG, due to a change in the estimated selectivity patterns. It appears that the
WFH-OZ model can easily become trapped at a local minimum for the objective function,
which leads to this alternate solution, rather than finding the global minimum. Besides
identifying unstable models, sensitivity analysis can also provide clues for analysts as
they seek the best fit for an unstable model. Sensitivity analysis can reveal to which
parameters of the model important outputs are most sensitive to change. It is critical for
analysts to try a wide range of starting values and bounds for those parameters in order to
help ensure that the global minimum for the objective function is found each time the
model is updated. Failure to find the global minimum can lead to dangerous
management decisions] (e.g., setting harvest limits based on an overestimated projected
yield, as in WFH-02). To this end, I have created a program for the MSC to automate my
sensitivity analysis, using AD Model Builder software (ADMB 2002). This program will
allow analysts to more easily evaluate the sensitivity of the stock assessment models
whenever the models are updated.

All of the stock assessment models were sensitive to changes in the stock-
recruitrnent function’s parameter bounds. Decreasing the density dependence parameter
bounds led to better fit parameter estimates for the WFH-Ol and WFH-04 models, which
significantly reduced the projected TAC/HRG in both models. It appears particularly
important to do sensitivity analysis using a range of starting values and bounds for

recruitment parameters each time the models are updated.

25

Sensitivity analysis also can reveal patterns of sensitivity across models, which
may point to assumptions about underlying basic model structure (i.e. the way various
biological, fishery, and data producing processes are included in the models) that should
receive more attention. I found that the Lake Huron stock assessment models were
sensitive to my assumptions embodied in stock-recruitment functions, gear selectivity,
and assumed probability distributions used to define the likelihood functions. Of
particular importance, the WFH-02 and WFH-04 model test quantities underwent similar
changes, which resulted in increased TAC/HRGs, when the Beverton-Holt recruitment
function was employed. A number of authors have considered the consequences of
assuming different stock-recruitment relationships to the management advice stemming
from those assumed relationships (Myers et a1. 1994; Barrowman and Myers 2000).
Other authors have discussed the relative merits of estimating stock-recruitment
parameters inside stock assessment models versus outside them (Maunder and Deriso
2003). The issue here is somewhat different than is addressed in that work since I was
only considering short-term projections. My concern here is more on whether including
stock-recruitment functions as a form of a “prior” influences and potentially improves
estimates and resulting short-term management advice, given a harvest policy exists. A
simulation study, where either freely estimated recruitment or priors based on different
recruitment functions were used would allow for a more detailed analysis of how
different approaches to estimation of recruitment fare.

I did not explicitly consider alternative approaches to estimating selectivity.
However, all of the models showed sensitivity to changes in gear selectivity starting

values and parameter bounds. While sometimes the changes were small, in other cases

26

changes were pronounced (e.g., WFH-02 and WFH-OS). These results reinforce concerns
that have arisen about the general suitability of the current double logistic selectivity
function during the development of the stock assessment models. As a result of problems
encountered during the original models’ development, all the selectivity parameters are
estimated in only one of the Lake Huron lake Whitefish assessment models (WFH-Ol). In
each of the other models, some of the selectivity parameters must be held constant in
order for the models to even converge on a solution. Thus, issues clearly go beyond
simply finding the best starting values and parameter bounds. Reduced or constrained
versions of the double logistic are not the only alternatives. For example, logistic curves
(Punt et al. 2001), double logistic curves (Methot 1990), gamma-type fimctions (Deriso et
al. 1985), and polynomials (Fournier 1983) have all been used to model selectivity.
Kimura (1990) and Radomski et al. (2005) found that use of an inappropriate selectivity
function can greatly increase the error in modeling results. This is another area where a
simulation study could be used to help evaluate the current and alternative approaches to
modeling selectivity. Alternatively, an empirical selectivity experiment could be used to
determine the actual gear selectivity, but this would need to be done for both gill nets and
trap nets.

All of the models showed some sensitivity to changes in the negative log
likelihood function, both when likelihood emphasis factors were altered and when
alternate distributional assumptions were made. In theory, if the assumed standard
deviations for the natural mortality, catch, and effort data and the assumed maximum
effective sample sizes for the age composition data are correct, then all of the likelihood

emphasis factors should be set to one. Methot (1990) warns that sensitivity to changes in

27

the emphasis factors indicates that there is some inconsistency between the data sources
or that some process has not been modeled correctly. Sensitivity analysis of the models
to the likelihood emphasis factors should be tested whenever the emphasis factors are
adjusted during model updates in order to report this sensitivity along with model results.
Replacing the lognormal likelihood function with the gamma likelihood function led to
only small changes in the test quantities. Cadigan and Meyers (2001) found similar
results when comparing the two likelihood functions, although they emphasized that the
gamma likelihood function is more robust to invalid distributional assumptions than the
lognormal. Williams and Quinn (2000a, 2000b) successfully used the Dirchlet likelihood
function to represent age composition data for Pacific herring, where sample sizes were
large. Replacing the multinomial likelihood function with the Dirchlet likelihood
function led to some changes in the test quantities, particularly the TAC/HRGs, in all of
the models. Again I believe a simulation study could be used to evaluate the robustness
of assessments based on these alternative distributions, and to evaluate potential
approaches to selecting between them.

In conclusion, sensitivity analysis provides a useful tool for analysts applying
stock assessment models. Running a sensitivity analysis whenever models are updated
with new data can reveal unstable models which are highly sensitive to change.
Furthermore, such analysis can identify particular parameters or assumptions that
generally have a large influence on outputs of interest. This can help focus attention on
these aspects of the assessment models. Such attention could come in the form of using
simulations to evaluate performance of alternative modeling approaches or collecting

new kinds of data to distinguish among modeling choices.

28

References

ADMB. 2002. AD Model Builder, version 6.0.4 edition. Otter Research Ltd., Nanaimo,
British Columbia.

AHWG (Ad Hoc Working Group). 1979. Second report of the Ad Hoc Working Group
to assess stocks of lake Whitefish, bloater chub, and lake trout in treaty waters of
the upper Great Lakes—State of Michigan: Lake Huron. Memeo. Rep. 38 pp.
with appendix.

Barrowman, N. J ., and R. A. Myers. 2000. Still more spawner-recruitrnent curves: the
hockey stick and its generalizations. Canadian Journal of Fisheries and Aquatic
Sciences 57:665-676.

Beverton, R. J. H., and S. J. Holt. 1957. On the Dynamics of Exploited Fish
Populations. Chapman and Hall, London.

Cadigan, N. G., and R. A. Myers. 2001. A comparison of gamma and lognormal
maximum likelihood estimators in a sequential population analysis. Canadian
Journal of Fisheries and Aquatic Sciences 58:560-567.

Clark, R. D., Jr., and K. D. Smith. 1984. Quota calculations with SAP.1 (Stock
Assessment Package). MDNR Institute for Fisheries Research, Ann Arbor.

Deriso, R. B., T. J. Quinn 11, and P. R. Neal. 1985. Catch-age analysis with auxiliary
information. Canadian Journal of Fisheries and Aquatic Sciences 42:815-824.

Ebener, M. P., J. R. Bence, K. Newman, and P. Schneeberger. 2005. Application of
statistical catch-at-age models to assess lake Whitefish stocks in the 1836 treaty-
ceded waters of the upper Great Lakes. Pages 271-309 in L. C. Mohr and T. F.
Nalepa, editors. Proceedings of a workshop on the dynamics of lake Whitefish
(Coregonus clupeaformis) and the amphipod Diporeia spp. in the Great Lakes.
Great Lakes Fishery Commission, Technical Report 66, Ann Arbor.

Fournier, D. 1983. An analysis of the Hecate Strait Pacific cod fishery using an age-
structured model incorporating density-dependent effects. Canadian Journal of

Fisheries and Aquatic Sciences 40:1233-1243.

Fournier, D., and C. P. Archibald. 1982. A general theory for analyzing catch at age
data. Canadian Journal of Fisheries and Aquatic Sciences 39:1195-1207.

Kimura, D. K. 1990. Approaches to age-structured separable sequential population
analysis. Canadian Journal of Fisheries and Aquatic Sciences 47:2364-2374.

29

Maunder, M. N., and R. B. Deriso. 2003. Estimation of recruitment in catch-at-age
models. Canadian Journal of Fisheries and Aquatic Sciences 60:1204-1216.

Methot, R. D. 1990. Synthesis model: an adaptable framework for analysisi of diverse
stock assessment data. Pages 259-277 in L. Low, editor. Proceedings of the
symposium on applications of stock assessment techniques to gadids.
International North Pacific Fisheries Commission, Bulletin 50, Vancouver.

Morgan, M. G., and M. Henrion. 1990. Uncertainty: a Guide to Dealing with
Uncertainty in Quantitative Risk and Policy Analysis. Cambridge University
Press, Cambridge.

Myers, R. A., A. A. Rosenberg, P. M. Mace, N. J. Barrowman, and V. R. Restrepo.
1994. In search of thresholds for recruitment overfishing. ICES Journal of
Marine Science 51:191-205.

Pauly, D. 1980. On the interrelationships between natural mortality, growth-parameters,

and mean environmental-temperature in 175 fish stocks. Journal du Conseil
39: 175-192.

Punt, A. E., D. C. Smith, R. B. Thomson, M. Haddon, X. He, and J. M. Lyle. 2001.
Stock assessment of the blue grenadier Macruronus novaezelandiae resource off
south-eastern Australia. Marine and Freshwater Resources 52:701-717.

Radomski, P., J. R. Bence, and T. J. Quinn 11. 2005. Comparison of virtual population
analysis and statistical kill-at-age analysis for a recreational kill-dominated
fishery. Canadian Journal of Fisheries and Aquatic Sciences 62:436-452.

Williams, E. H., and T. J. Quinn. 2000a. Pacific herring, Clupea pallasi, recruitment in
the Bering Sea and north-east Pacific Ocean, 1: relationships between different
populations. Fisheries Oceanography 9:285-299.

Williams, E. H., and T. J. Quinn. 2000b. Pacific herring, Clupea pallasi, recruitment in

the Bering Sea and north-east Pacific Ocean, 11: relationships to environmental
variables and implications for forecasting. Fisheries Oceanography 9:300-315.

30

 

m0$wqm $0.30 00.0 2.0 3.0 02.0550 m00.0 mm.0 3-3m?
vw.vmo.m 03.mmm 0m.0 2.0 N00 0m _._0m.m mmd 0N0 00.05?
383.0 ~0de 00.0 0N0 3.0 0_w.000.m mad 000 No.2:
3.03.0 09.03” 5.0 2.0 wmd 0m~.3m.m SWO 50.0 00-3;

 

00258.5 GmmOaa. 008 mmmm mmmm mom Mmmm 027.95 3955 m 20 m 7:. 6002
800530 00002030

 

.Som

5 SEE 803 me 3063 0o0oo-bmob 0mm~ 20 .80 26008 808383 xoowm 508005 83 0000080: 20 80¢ mos—g 002280702
0339: 20 05“ A30 OMEGA; 088 .800 6:8 mmmm 630008 328.85 85808 “a Mmmm dove—anon 0050:: new mmmm A30
mmeEoB cove—smog Am 200 0.2558 Soffiw 00828 30.0 .E 2.5 3:388 8?me 088—8 300 new 835» 08306.5 .2 033.

31

 

m0 m0- m0 m0- m0 m0- 0.0 0.0 0.0 0.0.. 0.0- 0.0 m0- m0 8080020008000 20

 

2 2 2 30 am 2. 8- 3 u 0882.0“ :80 20

9o 8 S. 80 od 8 ed as ed as 8. ed as 30 00892003 20

98. SN 30 30 as 8. as ed cs 8. ed as as as 08-3-2003 E
3.53

no no- 8 Qo ed od 8 as no me. 2- no mo- 2 0589.? E8 20

I No No 9o 3 8. e. _- 3 u asses? :80 20

9.2- 0.2 o0 0.0 30 od 3 od ed as 00 od 80 3 ”093-3003 20

2. S 3. as as 8. 3 od od 30 9o 8 ca 30 0093.283 E

3-83

0.: Sm :0- :0- :0- :0- 3 as 2. 3. S 2 3- to. 088%.. not... 20

Sm 2- no- 8. 2. 2 :0 3 n cause? 8.8 78

mm. 3 ed ed as 8. 30 9o od od 30 0.0 30 8. 0090-200; 20

02. m2 ca 0.0 as 80 cs 30 ed 80 9o 8 od od 08.3.80; E
8-53

2 3- cs 3- 80 :0- od od 80 ed to- no 3 S- amass? :00... 20

um om om om om om om o; H 0:083:30 fiche ZO

v.2- 4.2 as 30 ca 3 ca 8 8 cs 30 ed as 30 0090-200; 20

3.- 3. cs 3 ca 80 8 0.0 3 ca 30 od 8 as 002.8003 E
5-53

- + - + - + - + .. + - + - + 08000 .8 :00888Q
088.3 as. 8% 80m .00 00mm 8580 8820 a 20 a E

 

800880 80020310
.800 090380 9 =8 8 20008 05 00. 000880 080m .8003 2.08% 3 8m 08 .3 00000800 A+V 000888
083 00200, 808 800 8:3 . Son 8 883 0005 .8 8083 0008-388 0mm: 05 .80 20008 80808000 x080 000823 000: 05 80.0
050.3“. 0008.88 000 .008 mmmm 6800000 3:388 00:08.80 8 mmmm 80003000 00:00:: 80 Mmmm £00805 800—800 .E
260 3:888 88:8 0080—00 300 .00 2.5 3:888 80.080 0080—00 3.3 80 82? 08—008 80¢ 0088000 80080 .NA 0300.

 

32

 

 

33

00 00 0.8 0.2 0.0 0.0 00 no n 2080 02082.0

3- 3 0.2 4.:- 3... 0.9.- 0.2 000. 40- v0 00 00- ? 00- 2820. 00302
8-003

2- 00 2: 4.0- n: 4.2- ...m 3- 3- no no 00- 00 00- 083880 08E

00 0.0 00 00 0.8. 000 0.8. 020 00 00 00 00 0.0 00 2080 80.880

0,0 00 0.0 0.0 0.0 0.0 00 m0 M 20800 80800.5

3; 00 30 0mm- 4.8 me.- 000 0.00. 00 E 3- m. _- Q- 00- 2.62.... 0.5.2
20.003

00m 0.0m 0m 0...- S E- 00 E- 3. 8 2 Z S- 3. 0503830 2.8

NR 0.8 00 00- 3...- at. 3:.- 004 4.0 S 2 2 40- m0 0080 80820

0.00 to- 0. 8. 000 E E 3 no u 0080 ass-58$

a 0.00 a mom. a S...- 0 SN- 0 am e S e no- 23808 00302
8-003

40. no 00 2- 0.2 3 7 S. 3.- do- 3 00 00 :0 00 055.880 08E

00 00 0.0 00 000- 0.8 000- 000 00 00 00 00 0.0 00 28.0 0288.0

00 00 0. 0 m 0. 0 m 00 0.0 00 no u 2.80 02:88.0

00- .2 3.0 3m- 0mm 00..- m.m~ 3.0- 3- 00 S :0- 2 S- 28008 00:02
8-003

- + - + - + - + - + - + - + 000000 00 8000005
0008.8 8a 0000 0000 .00 008 020.000 .8620 0 78 0 E

 

80.00000 00002030

 

A000 0w8>000 00 :00 00 20008

05 00— 008000 080m .000 08088.0 080 08 00 8500 000.000 083 0020.» 3:808 00008 3 0000800 0 0:03 .000 080—0 08: 05 8
a: 03:00 0.83 830., 38508 8008 A5 000088 :0 6300000 38008 8 m .820..- 008000 00 80 0:0 .3 00000800 A5 0000088
0003 0030> 808 S00 0003 ._00N 8 8.82 0000 .8 80003 00000-3000. 022 05 8.0 £0008 808000000 00000 0000003 0002 05 800
Gym-<92. 0008.88 000 .008 Mmmm 60000000 300088 00:08.00.— 00 Mmmm 60003000 000008. 8.0 .mmmm 000805 8003000 .E
200 30388 80-5w 0080—00 300 .E 2.5 300088 000.008 0080000 300 8.0 820> 0000000 88-0 000080-00 80080 .m._ 030.2.

 

0.0-
06

5o
od

m0
od

0%..
od

m.:
mdm-

0.0-
mdm

0.0
mdm-

m0-
mdm

0.0- 06 md We-
od od od od

0.0
od

5o-
0.0

90:30am .00 085.
208$ 000.8005

00.0080 2 0300

34

 

_.N- .0.— Nd- 0o No- 0.0 o.o od .0.— NA mdm Nam—- w; NT w0000000000>00020

 

3 3 3 0.0 00 0.0- 0.0. 3 n 00020: 02:2 20

0.00 2:- 3 0.0 S 0.0 00 0.0 000 3:- 0.00- 0.00 S- 2 02022 032 E

:- 20- 00- 0.0 0.0- 0.0 2 3 n 02020: 02:2 E-
3.53

000 0.2- .3 00 No 0.0 0.0 0.0 2: 0.0 0- 00 0.0- 0.2- 0.: 02082 02:2 20

0.0- S- 3- S. 0.0 0.0- 0.0 3 u 00022 00:2 20

0.0 3- 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0- 0.0- 0.0 2: 02- 00020: 02:2 E

0.0- no- 00- 0.0 0.0- 00 3.- 3 n 00020: 02:2 2..-
2.00.0

0.00 00- $- 2 2- S 0.0 0.0 0.2 3- 00 S- 0.:- 2: 00022 00:2 20

EN 00 to 0.0 0.0 No- 0.0 3 n 20:20: 02:2 20

0.00 0.00 0.0 0.0- 0.0 0.0- 0.0 0.0 0.2 S- 0.0- E 0.2 2:- 0005.0: 032 E

0. _ m wo- 0.0- 0.0 0.0 00 0.0- 00 n 002000: 02:2 7:-
8-003

0: E- 0.0 0.0- 0.0 0.0- 0.0 0.0 3: I _- 00 0.0. S:- o.: 00:20: 02,2 20

$- 00- 00- 0.0 0.0. 0 :- 00 3 u 00020: 02:2 20

0: 0.0- 3.- mm :- 00 0.0 00 E $- 3- S 0.2 S _- 05020: 02:2 E

00- S- 3- 0.0 0 0- 0.0 0.0- 3 u 02020: 02:2 E
8-00.0

- + - + - + - + - + - + - + 0w0000 00 000000009
00065 20: 0000 0000 .00 0000 0200:: 00820 0 20 0 7:

 

00000000 00002020

 

.0022, 00.00000 00 000 000 .3 000000000 A+v 00000005

0003 000_0> 00000 800 000.5 . Sow 00 000000 00—04 00 000003 00000-5000 02: 000 00-0 200000 0000000000 00000 0000003 00.00 05 800
OMEU<H 0800.800 000 .0000 mmmm 6300000 300008 000000000 00 mmmm 00003000 0000000 000 Mmmm 000805 00003000 A0
20V b00008 000-_=w 0000200 0:00 A0 2.5 b00006 000-000 0000200 300 .000 0020> 0000000 0000 000000.000 000000 .0; 0300-

35

 

 

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0003 000000 00000-0000 000002
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 000000 0000000000 000002
00-20?
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0- 0.0 0000.» 000000 00000000 000002
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 000000 000000000 0000002
00-00005
mdm 0.0m 0.0- 0.0- 0.0- 0.0- 0.0 0.0 0.0 0.0 NA N0 Nd- md- 00_0> 000000 000000000 00002
0.0m 0.0m 0.0- 0.0- 0.0- 0.0- 0.0 0.0 0.0 0.0 m0 m0 No- Nd- 000000 00000000 00002
8.00005
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 00_0> 0000 000000000 00002
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 000000 000000000 000002
5.30?
- + - + - + - + - + - + - + 0w0000 00 000000000
OED/0.0. 0000 ~53 mmmm 000 0mmm 000000D 0000005 0 20 0 Z0-

 

00000000 00000—0>m

A-v 000000000 000 A+v 000000000 0003 0000000000

000000000 000000 0003 . Son 00 0000: 0x04 .00 000003 00000-000000 03— 0.0 00-0 £00000 000000000000 00000 00000005 0000— 000 0000
00EU<H 0080.800 000 .0000 0mmm 6300000 000000000 000000.000 00 Mmmm 00000000 0000000 00.0 Mmmm 0000005 000000000 A0
20v 000000000 000-Ew 00000—00 0:00 A0 2.5 00000000 00-00-00 00000—00 0:00 00-0 000—00, 0000000 .000 0000000000 0000000 .m._ 0300-

36

 

06 ed od o.o ed 06 ed od od od o.o od ed 06 egg—32823252030 ZO

 

o.o od 3. ed 3. S0 3 ed 3 9o 3 od 3 ed €58 £8228 20

3 od o.o 3 9o 3 9o 3 od 30 3 9o o.o od 02.? Ea: £5238 7:

0.0 0.0 06 od 06 0.0 od o.o 0.6 0.0 o.o 0.6 o.o o.o £0555 5253.330 E
nonmmk/

od o.o od od ed 9o 3 od 3 od od ed 3 ed 83> :38 £5228 20

0.0 Cd o6 0.0 0.6 Cd o.o od od 0.0 Cd o.o od od $08905 5259.330 20

3 od ed 3 3 9o 3 ed 3 9o 3 od 3 ed 02? :58 £5228 8.

ed 3 3 3 ed 2 3 2 od 3 3 ad 3 0.0 $58 £5228 5
v0.3.2},

06 hfim 06 Ad- Cd Ad- od 06 od 06 o.o NA 0.0 Nd- os3> 3252 3253.330 20

him 54% Ad- Ad- Ad- Ad- o.o Cd 0.? 0.? NA NA Nd- Nd- $0555 5253.330 20

5.va ficvm Ad- Ad- Ad- Ad- od 06 9v 0.? NA NA N6. N6. 033> 3222 52532030 7C.

him 06 Ad- o.o Ad- o.o od od 9v 0.0 NA od N.o.. od £0555 b2532030 7C.
Noummk/

ad ad od 06 ed od 3. ed 3. 3 od 3 9o o.o 28> 8:5 £5238 20

od 3 Qo od ed 9o 3 3 3 ed 3 ed od 3 258 £5228 20

ed od 3 3 ed 3 ed 3 od 3 8. 3 ed 3 2;? :38 £5228 E

od 3 ed 3 3 9o 3 3 od od .3 3 od 3 2:58 €38.28 E
29.3.23

- + - + - + - + - + - + - + owqfio mo 8388809
CEUAC. mmmm mmmm .32 mmmm 20028.53 wwwEomm ”A 20 ”A 75.

 

$20830 =282§m

 

.3 338020 28 A+v 393.58 803 88083.80

325288 80:? . Sam 8 :98: 8:3 .3 8895 3304338 32 05 .80 2288 80833.3 0.8% 5282.5 83 2: 88m
0EU<H 338.88 28 .028 Mmmm 633:8 ©2808 3896on E .mmmm down—smog gamma: .52 mmmm £3803 corn—smog AHA
ZOV .3388 027:3 33028 53A .E 2.5 52.8.88 008-98 vogue—om 38% no.“ 329, 082093 80¢ 3:80.23 838m 6A 038.

37

 

:.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: 88:8 8:288 8:28

:.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: 28> 8:8 888 8:288

:.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: 888 880.... 8:288
8-8;»

:.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: 8:88 8:288 8:28

:.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: 28> 8:8 888 8:288

:.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: 888 888 8:288
8.8,;

an: :.: a: :.: N: :.: :.: :.: :.: :.: :.: :.: :.: :.: 88:8 8:288 8:28

5.: :.: _.:- :.: _.:- :.: :.: :.: :.: :.: 8 :.: :.:- :.: 28> 8:8 88:... 8:288

S: 5.: S. :.:- 3. 3. :.: :.: :.: 8. S S 8:- :.:- 888 888 8:288
8.53

:.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: 88:8 8:288 8:28

:.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: :.: 28> 8:8 882.. 8:288

 

0.0 0.0 0.0 0.0 0.0 0.0 o.o 0.0 0.0 0.0 0.0 0.0 0.0 o.o 308002 80—03 80320808
3.2%»
- + .. + - + - + - + - + - + 03020 ,3 8008.309
03254:. 0:8 mmmm mmmm .808 Mmmm 803283 38805 ”A 20 .2 7E.

 

30238000 coumagm

 

.3 030.0800 08: TV 0030808

0:03 3020.» :88 80A: 8023 Jean 8 80.82 0023 .«0 30:03 0000043808 on: 05 80.: 20008 8083030 2003 320033 003 05 808
OMEU<H 00803.8 08: .038 Mmmw 630053 £3.58 00808000: 00 mmmm 8028800 0038.30 :00 Mmmm £3805 8020—0008 Am
200 AHA—8.88 008-2% 00000.3 5:0.“ Am 25 328.808 008.980 00000—3 53% .5: 302.8, 08—33 80.2 008080.22: 8080a SA 030,—.

38

 

dd dd dd dd dd dd dd dd dd dd dd dd dd dd 0=_0> 20222880823002

 

dd dd dd dd dd dd dd dd dd dd dd dd dd dd :00002.n\ 80828002
dd dd dd dd dd dd dd dd dd dd dd dd dd dd 020:, 2028:»: 80888002
:.:- dd d2 dd d.: dd dd dd 06- dd md: dd dd: dd :00002 8 80888002
moi-:3
dd dd dd dd dd dd dd dd dd dd dd dd dd dd 030> 20888 80888002
ddN- dd 0N dd 0N dd dd dd :.N- dd N.: dd md- dd :000020\ 80888002
dd dd dd dd dd dd dd dd dd dd dd dd dd dd 0:20> 202228 80888002
dd dd dd dd dd dd dd dd dd dd dd dd dd dd :000028 80828002
3-38:3
0.0m 0.0m :d- :.d- :d- :.d- dd dd 0.0 0.0 N.: N.: Nd- Nd- 020> 202828 80828002
dd 0.3. dd :d- dd :.d- dd dd dd 0.0 dd N.: dd Nd- :00002.n\ 80888002
0.: 0.0m 0.: :d- 0.: :d- 0: dd 0.: 0.0 0.: N. 2 0: Nd- 020> 2022:: :0 80888000:
0: N.0m 0: dd 0.: dd 0: dd 0.: 8m 0.: :.: 0.: 0d- :00002 c 80888002
No-22?
dd dd dd dd dd dd dd dd dd dd dd dd dd dd 020> 888880888002
0.?- dd 0N- dd 0.N- dd dd dd N.m :- dd wd dd dd dd :000008 80828002
dd dd dd dd dd dd dd dd dd dd dd dd dd dd 00:0> 202828 80828002
Nd- dd 0d- dd dd- dd dd dd 0d- dd Nd dd dd dd :000028 8088882
522-23
- + - + - + - + - + - + - + 0w0020 :0 00202085
02:25.3: 0::0: 22mm mmmm :02 Mmmm 0028-83 80805 n: 20 n: 2.:-

 

:0:::800O 0020203.:
.05 0w:0>000 0: 20: 0: 20008 02: 00: 88020 080m .3 00:00:000 000 A+: 00:00:00: 0:03 8808800
8082800: 0023 . :ddN 0: 008m 020-: :0 80:03 000004308: 0mm: 02: :0: 20008 808::0::0 200:: 2:20:83 020: 02: 80::
Ova-<05:- 00:00.:0:0 000 .020: mmmm 680020: 820808 0000:0.:0: :0 mmmm 80200800 0028-80 :0: mmmm .::0805 00208000 A”:
20: 320808 :00-28 008020: >20”: .0: 2b 520808 80.00:: 008020: 2:20“: :0: :020> 0:20:02 80:: 00080.28 8080:: d: 030.:-

 

39

 

od od od o.o o.o od 06 ed od od ad ad 0.0 o.o egg EEENQAom 20

ed od o.o od cd o6 o.o od o.o od o6 o c od od mugon mm ._om 20
um o.o om ed om ed on ed on od ow od um 06 2:? REE Nu ._0m 20
_.N ed @5- ed 0&- ed od o.o 5o od md- o.o wan- ed 3559 Nat—om ZO

vo-Em?

new od #0. od fid- o.o o.o od 06 ed NA od Nd- od o=_m>_mE::Q._om 20
0.3” scum fio- _.o- fie- fio- od od of 06 NA N4 N6- Nd- oin>fi£fim§g8 20

c o om o o om o.o o.“ o.o um 06 om od om od om $358 mm ._0m 20
bém o.o _.o- od fio- o.o o.o c6 9v o.o NA o.o Nd- o.o 25g 3:5 mm ._om 20
5.3“ od fio- od _.o- od od od 04» o.o .0..— o.o Nd- o.o mngop Nu. .3m 20
5% 5.3“ fio- _.o- _.o- _.o- od od 06 we NA N; No- Nd- ofig BEE Nu .6m 20
:1 5.3” m.w- fio- Ow- fio- 06 ed v.2- 9v o6. N; 0.0 Nd- mnason 330m 26
mo-TEB

od om o.o om o.o om o.o om od ow od ow o.o om 03:; REE 1:8 26
ed _.o od od o.o od 06 ed od o.o od o.o od o.o 0513325 338va

 

o.o od 06 ed o.o od o.o o.o o.o o.o o.o o.o od ed 8569 Md. .3m 20
ed od o.o o.o od o.o od od od od od od od o.o 2:? REE mm .3m 20
ed od o.o od 06 o.o o.o od 06 ed od od o.o o.o WEBB mm .3m 20
o.o od od od o.o od o.o od od ed o.o od od o.o 33> BEE E- .Em 20
ed od o.o o.o o.o o.o o.o od od o.o o.o cd ed 06 £525 Nah—om ZG
Sim},

- + - + - + - + - + - + - + omSEo mo coumtomoa

GMEU<H oufi mmmm mmmm mom mmmm 8595 3285 m 20 m Zu-

 

momESBO cocaigm
.Aob 09268 8 2mm 8 20on of to. mowSfio oEom .QS ago—m vacuum was Ammo “Eon couoocfi
288m ANS ego—m ram .33 «Eon couoomfi SE 05 203 £808an 8:83 33628 2:. .3 c8353 28 at @8885 v.53
880.53g b33028 5:3 . Sow E 553 8:3 mo £895 Ego-bag cm: 05 8.“ 2258 Eon—magma x8; ammBEB 83 05 89a
OED-«H 380.85 Ea “0:8 mmmm 633:8 £3.88 8:96on 3 mmmm dorm—smog Baa-«S. .8.“ mmmm £8805 5:233 .E
29 mum—888 “oz-5w vogue—om ba .3 2.5 £3.88 «oz-m8» 3828 $3 8m 32? 05—83 88m oofizobmc “58$ d; 033-

 

4O

 

wfim 5.;

nd ed od fin o. T
2 o.o o.o adm- «.m-
85883 a ma 38833 8:
m._ o.o o.o Em- 9v
~82:qu a as 38838 “on
v.3 o.o od 3% 0.:
88883 a mm @8358 “on

ad od od od od

o.o od od od od
cd od o.o o.o od

mam
o6

g;-

od
od
od

2.3, 325 E .E zw
2:? 3:5 3 ._8 zo
£58 3 .E 20
02.; 3:5 E .E 20
€38 mm ._a 29
02? 325 E .E 20
$58 E 43. 20
8-53

02? 3%: E .E 20
”as, 325 R ._3 zw
$58 3 ._8 20
.9253 3 23

41

 

od od od od od o.o od o.o od od od o.o. od o.o 33> 328 “um .30 7E.
od od 06 ed od o.o o.o ad ad ad od o.o od od 33> 338 mm .30 7C.
od od od od o.o o.o o.o ed od od o.o od o.o o.o 33> 33:: NR .30 7:-
o.o 06 ed od od od o.o od od od od o.o od o.o 3302 ~&._0m 7E-

meg-£3
o.o od 06 ed od ed od od od od od o.o ed o.o 33> 328 vi .30 7E-
od 06 ed od o.o od 06 od 06 9o o.o od o.o od 33> 3E8 mm .30 2,—-
o.o od od od o.o od 0.0 od 06 ed od 06 ed od 33> 338 ~& .30 E
fifi od m.m- o.o Wm- o.o 0.0 od od od ONT o.o QC- od 3859 ~Q._0m Z-H

3-3m?
fivm o.o fio- od fio- o.o o.o ed 06 ed NA o.o Nd- od 33> 3E8 23—00 E
hem Baum fio- fio- fio- fio- od od we 9v NA NA md- Nd- 33> 338 mm .30 E

3m 00 :0- 00 3- 00 od 00 3 a S 00 No- 00 03? 808 E .8 E
od 3. od Qo 0.0 o.o od 3 ed 3 od 0.0 0.0 3. $58 E .8 E
8-8,;

06 ed 06 ed o.o od o.o ed 06 06 od 06 fio- fio 33> 338 wax-.30 7E-
fio od fio- od fio- od o.o od fio o.o Nd- o.o v.0- fio 33> 388 Mm .30 7E.
od o.o od 06 ed od 06 ed od od 06 ed od ed 33> 338 3:00 7C-

 

06 0.0 06 Cd od o.o od od od od 06 0.0 0.0 Cd $0559 NQAom E
:33;
- + - + - + - + - + - + - + 0330 00 82880000
UMEO<H 038 mmmm mmmm mom mmmm REE-ED mmmEomm m 20 m 7:

 

002830 aou33>m

 

.3: 090800 9 =8 9 2308 05 30— 0.0330 080m $3 320 30000 23 Ammo 86m 8300?:

88000 ANQV 320 8a .CS 83a 83038 8m 08 0003 830883 80303 8380—00- 05- .3 308083 83 A+V 3003088 0803
80~0833Q £38030 80:3 . Son 8 80.53 0034 mo 8033 3300-5308 0mm“ 05 .8.“ 2308 808000003 03on £308? 003— 05 89a
O§U<H 3080.88 83 .830 Mmmm “0:60:00 33:08 03080.38 3 mmmm 823183 3058.8: Ho.“ mmmm £00805 823393 .E
20v 328.88 “0.8-Em 3080—00 33m Am Zb 33808 “3.938 3080—00 33% Mom 033> 3:003 80¢ 00:0u0t3 8080A .o fl ._ 073,—.

42

 

E- 8 . 3 8 3 8 0.0 8 38- 8 m8 8 mm- 8 .88288E
3 3 od 3 08 3 08 3 08 3 08 3 3 3 83858833582

 

3.33
8.8 $- 2 38 no to 3 3 3.8 3.2- 3- 03 3- 3.2 .85 888 .88 20
gm :8 no- :0- 2- 3- od 3 3 3 3- 2 3 no- .85 8880 20
8.2 3- 3 S- 3 3- 08 08 3 3. 3m 3- 2 N8 .85 588 20
3m in S- :- 3- 2- 88 8.0 8.0 S a: S- 3. 3- .85 9:8 own E
8.2 8.2 8. _- 2 8. _- 2 3 8.0 2- 8. : 3 @8- od 3. .85 8880 E
mom «.3. no- 3 no- 3 0.0 8.0 8 3.0 3. 3- no 2- .85 588 E
38 $8 3- S- 3- 3- 88 3 3 3 S 2 N8- N8. 85 85888 8:82
8.33
0.: 8.2- 3 3- o.” 3. ed 3 3 0.2- 31 3 88- I: .85 888 088 20
3 2- mm- M: 2- M: 3 3 E 3- 3- 8.8 2- 8 .85 8880 20
no- no 38- no 38- no 3 3 no N8- 2 3 Z.- S .85 588 20
38 8 to 8 to 8 3 8 3 8 2: 8 38 8 .85 888 088 E
0:- 8.: 2- 3 E- 3 3 3 E _- I: G 3- 8.8 3- .85 8880 E
$- 8 38- 8 to- 8 0.8 8 3- 8 M: 8 N8 8 .85 888 E
3 3 3 3 3 88 3 3 08 8. 08 3 88 8.0 .85 05888 8882
5-883
- + - + - + - + - + - + - + 0383 .80 aozmtom0m

Gym-<92- 038 mmmm mmmm 80m mmmm @3385 @3805 m 26 m 72-
. 805330 :03:3>m
.35 030300 8 =8 08 20308 05 30— 83380 080m .3 30308003 83 fl; 3030008 0803 80803-8
83380 3005302 80:? . Son 8 808m 0034 .80 8033 3300-3308 cm“: 05 80.8 £0308 3083083 0303 :30an 003— 05 80¢
OMEU-cd- 30800.88 33 .038 Mmmm .3638 $3808 03080.38 3 mmmm 8033808 3088-83 80.8 mmmm 33805 803183 Am
29 $3808 809:3 30303.8 b3 Am 2.5 33808 808.38. 303038 DE 80.8 033> 08330 80¢ 030808-83 3080a A : 23H

 

 

43

 

_.m~
as
3‘
v.2-
n.3-
m.a-
dd

om
dd-
wd-

ad
«a-

c.2-
0.x.-
vd-
m.m~
a.a
dd
dd

m._ T
van
ad
md-
a.m~

9N-
d.m
v4
NS-
d.m
NA
dd

om
Nd
ad
om
v.m

NM
d.m-
ad-
v.»
WN-
ad-
dd

m.m
vd-
od-

fin
d.N-

9N-
d.m
v;
m8-
d.m
N;
dd

om
Nd
ad
ow
Ym

Nd
d.m-
ad-
vd
m.N-
ad-
dd

m.m
vd-
dd-

fin
dd-

dd
dd
dd
dd
dd
dd
dd

om
dd
dd
om
dd

dd
dd
dd
dd
dd
dd
dd

dd
dd
dd
dd
dd

ms
dd
v6
92-
Yd-
NS-
dd

om
N.—
_d-

om
a._-

ad-
ad-
5N-
a.a_
Nd
v6
dd

5m-
vd
md
Nd
mN

dd-
de
ms
v.3
md
dd
dd

om
9N-
ad

om
iv

NM“
3:
a.m-
fiN-
WN-
dd-
dd

fimd

_d
ad-
Wm-
ad-

mNT
9N-
m.-
5w

a.w~
:u
dd

om
N;-

“A

on
0.:

dad
Nam
04
d;
va-
WN-
dd

52
ad-
ad-
m.m-
dd-

.85 9:8 own 20
.85 tote 20
.85 :88 20

.8:— ano own 7E-
.BE tote 2,—-
92 n38 7:.

.85 €358 522

3.3:

.8:— ano own 20
.3: note 20

.8:— anu 20

.85 @800 um“ 72-
.BE tome E
.APEoov :4 2an

44

 

2 3 od 30 no no- 3 88888 88:82 8880
:.: 3 ca 3 2: 3.- S- 882: 8:8 E
08 a8 3 od 3: 2- 3:- 388 888850 8888:

um um ow um um um um teammbmcoog EoEflEooM

um um um um um om um um“ mo zofloca 30G: hwy-6550””

:8 3 od 3 ed 30 9o 8:8 :88 8 8388 8:82: 8882

88:3
8.8 3- am- o.o to- 2 S 88888 88:82 8880
on um um om om ow om HGoGOQEOQ @0258:— «EU

2 E E Q: 3 8.0m i 882: 8:8 E
8.8 to to 3 3 S o. 7 888 888880 8888:

8o 9: ed 3 2 :.o- :.:- 888888 0888.88
”.8 3 0.0 3 3 o: no- 08 8 888: 88: 8888:
SA 3- 3. ed 3 N: do- 8:8 :88 8 8888 8888 882

8-53
$7 0.: :.: 3 3:. SN 3”. 088808 88:8: 828:0

no 3 :.o :.: no 3. :.o 88888 88:3: 8:80

0d 0d 0d 0d 0d 0d 0d 000803 Egan HE

S no 2 od 2 no I. 888 Hem-8850 8888:

um om om om um um um cocmmbmcoog “cogmnEooM

M3 :0 :0 :8 3 9o 3 88 8 888: 88: 8:88,:

3 :.: 9: 3 q: :8 od 8:8 8:8 8 8888 8:808 882

_ d-EmB

 

DMEQE- 0:8 Mmmm mmmmfim Mmmm Baa-ED 000805 mZO m7:-

 

0028020 :oumfigm

0w§n0 do 2053.589

 

45

.35 030280 8 dam 8 20908 05 “02 m0w§20 080m 62-5008

003 0:30:50 3008 5:3 . Bow E 0852 00:5 do £0003 330-5005 dmfi 05 Sm 2088 80800030 0.0on £8.33?» 003 05 80¢
05003—- d0000_.o.a 2:0 .038 mmmm 6:60:00 3:00.58 00:20.30 00 Mmmm demand—Eon “00:0? .8.“ mmmm £00805 55033:: .E
20v 5:358 «Bu-Em “00000—00 bud Am 2.3 bag—08 “02-93 “00000—8 21a How 0029, 05—0009 Sod 00:0:0b6 0:020; .m_._ 030,—-

 

0.x 7 v. T v. 7 od 0.: T fim ad 80808800 802500:— 0020th
Yo :.o fio ed to 5o fio- 8088800 8008200.: 0885
w.m_ Nd Gm od 03 NS. QN- 80803 805$ 0E
wdm m6 mé o.o 5.3 5m- m.m- $008 :oz-cot0>0m 8088803:
wdm new we od 9mm m4- w.m- 808808008: 8088803:
:.m- 9o- 9o- od Nm- to o._ 0mm m0 808083 808— Q€§00m
06 ed o.o o.o o.o od ed 02? 3:8 8 88800 3:888 8:502
Sim?
m6 od 06 ed NN :1 _.o. 80800800 8005—00:— «02085

 

.3588 SA 050H

46

 

vv._wo.m 3m.~oo 93 :.o vmd ©3623” mood Ewe 6% 0:1? RES Va ._0m 20
8.3%,?

3.39m NwVJVoN Ed Ed mvd SWNQQN omd mmd 6% mwgonbfl E25383”
3-2m?

mmwovd m©~.mfim omd 2d wmd owmfimoam wmd Ed .oov €55an Eofizeuuom
5-3m?

 

83:35 Omguflw 2:8 mmmm mmmmgom mmmm cunmmcb @385 m 20 m E
333530 cozmflgm

owfifio mo :oumtomom

 

.86 333059“ 83 338%

EH 9.3 3355a? .83 Dragon—om duo—m 988m 05 mm? um. $8883 2285 $3628 2:. .E .258 uo>oEE_ mowcano

5:3 . Son 5 55$ 834 no @683 330-535 omfi 2t .8.“ 20on 308883 x08...“ :EBEB 8:: Qt 80a 83.? voonzoxn-wo_
gamma: 05 28 AB: Dagny»; @8832“ 6:8 Mmmm 6:628 3:888 85me E mmmm .coumamoa gamma: 8m Mmmw A35
mmmEoB scum—anon Am 20v >538:— §E=w 3828 31c .C 2.5 3:888 8993 vogue—om 3: 8m 829, 38:55 :2.— 2an

47

\

LA KE SUPER/0

 

N
° W
S

 

 

 

 

   
   

 

WFS-04
WFS-OS WFS-06 FS-07 ~. Canada
' JM ‘
Michigan “3'03
WFM-04 . “+03
WFM-Ol wmm . I . O
V . an ' ..
M-02 ' '0 d
‘ WFH-Ol wmm
LAKE fiflpH/GAN

 

 

(I? WFM-OS WFH-OS
LAKE

Michigan
HURON

 

 

\ ?

Figure 1.1. 1836 treaty-ceded waters and lake Whitefish management units in lakes
Huron, Michigan and Superior.

48

CHAPTER 2
EVALUATING METHODS FOR ESTIMATING PROCESS AND OBSERVATION

ERRORS IN STATISTICAL CATCH-AT—AGE ANALYSIS

Introduction

Modern statistically-based stock assessment models allow a stock assessment
analyst to explicitly account for process and observation errors. Observation errors
within statistical catch-at-age analysis (SCAA) commonly take the form of differences
between observed and true fishery catch or survey indices of abundance. Process errors
within SCAA generally take the form of annual deviations in recruitment, catchability, or
fishery selectivity. Errors within SCAA also can be combinations of observation and
process error. For instance, fishing effort can be predicted within SCAA using estimates
of annual fishing mortality rates on fully selected fish and fishery catchability. Some
analysts implicitly treat the deviations between observed and predicted effort as
observation error. In reality these deviations are due to interannual variation in
catchability (process error, which will often dominate) as well as errors in observing the
nominal amount of fishing effort. Similarly, deviations between model and observed
values of fishery catch per unit effort (CPUE) arise from a combination of observation
error and interannual variation in catchability. The variances associated with all of these
error sources or the ratios of those variances are used in SCAA to weight the different
data sources during the model fitting process (F oumier and Archibald 1982; Deriso et al.

1985).

49

It is important to understand how values for process and observation error
variances are obtained because those values can affect SCAA results. Deriso et al. (1985)
demonstrated that altering the assumed known ratio of catch variance to effort variance,
which they used to weight fishing effort and catch data, affected estimates of fully-
selected fishing mortality, surplus production and year-class strength of halibut. Chen
and Paloheimo (1998) found that misspecifying the ratio of catch variance to effort
variance could lead to increased estimation bias in catchability and natural mortality. The
National Research Council (1998) recognized the importance of correctly weighting
different data sources within stock assessments, and recommended that more research is
needed to determine how those weights should be set.

Process and observation error variance values can be derived either separately
from SCAA or estimated within SCAA. Derivation of error variance values separate
from SCAA is the more common approach, with these estimates or their ratio then
assumed known during the subsequent SCAA. A plausible estimate of observation error
variance for data subsets such as observed annual catch, effort, or abundance indices
often can be obtained through analysis of the raw data used to derive these quantities,
taking into account the sampling designs (Law and Kelton 1982; Sitar et al. 1999).
Process error variances cannot be estimated in the same way, by analysis of assessment
data subsets external to the model, because by themselves these data are not informative
about how model parameters such as catchability are varying. As a result, assessment
scientists often rely on expert opinion to obtain estimates of this component of variance.
Merritt and Quinn (2000) applied this expert opinion approach and other empirical data

weighting approaches to the assessment of a recreational fishery, and judged that the

50

expert opinion approach produced the best model based on an analytic hierarchy process,
a decision making technique. Since they were working with actual fishery data, Merritt
and Quinn (2000) could not evaluate the accuracy of the variance estimates produced by
expert opinion. With estimates (or educated guesses) of observation and process error
variances in hand, the assessment often then proceeds assuming these values or their ratio
is known. There are several potential disadvantages to such a two-step procedure. First,
uncertainty surrounding the error variance estimates is ignored in the subsequent SCAA
(Maunder 2001). Second, the reliability of process error variances based on expert
judgment can be questioned. Francis et al. (2003) discovered that the standard values
used in New Zealand stock assessments for the coefficients of variation (CV) for
commercial CPUE (effectively the CV for process errors in catchability), which primarily
are derived from expert opinion, typically were too low to be consistent with the resulting
interannual variation in assessment model estimates of fishing mortality. In contrast, they
found the prespecified trawl survey CVs were too large to be consistent with the resulting
deviations between assessment model estimates of catchable stock abundance and
observed survey indices.

Process and observation error variances generally are not estimated within SCAA
due to difficulties in estimating the variances as parameters. This task is particularly
difficult when multiple variances are being estimated. The potential advantages of
estimating variances in SCAA are that 1) all of the data in the analysis can be synthesized
to obtain the variance estimates, and 2) for some methods, uncertainty surrounding the
variance estimates can be quantified and accounted for in the analysis. Two main

statistical methods exist for estimating process and observation error variances in SCAA

51

(Schnute 1994). First, a Bayesian approach could be taken, in which prior information
about the process and observation error variances is incorporated into the analysis to
derive marginal posterior densities of the variance estimates as well as other parameters
and quantities of interest. Second, a mixed model approach could be taken, in which the
process errors are treated as random effects, rather than parameters, which can sometimes
allow for the estimation of both the process and observation error variances as model
parameters. Richards et al. (1997) suggested a third approach to estimating process and
observation error variances. This method requires a prior point estimate of observation
error variance. In essence the method is to repeatedly fit the assessment model, each time
using a different assumed known variance ratio, and choose the ratio that produces
deviations between observed data and model predictions that are most consistent with the
prior estimate of observation error variance. Unlike the other two methods, this approach
does not account for uncertainty in variance estimates in the analysis, and I refer to it as
the ad hoc method, because the approach to estimating the ratio of the variances was not
based on a formal statistical justification. To my knowledge, no previous attempt has
been made to compare these different approaches within SCAA.

My objective was to determine whether or not process and observation error
variances could be reliably estimated within SCAA. To answer this question, I evaluated
two different methods for estimating the variances associated with annual variations in
catchability (i.e., process error) and total catch (i.e., observation error) in SCAA. The
two methods I examined were the ad hoc approach described by Richards et al. (1997)
and a Bayesian approach. I looked at using both strongly and weakly informative priors

on catchability variance for the Bayesian approach. In addition, I initially attempted to

52

implement a mixed model approach using the random effects module for AD Model
Builder (Otter Research Limited 2005), but that estimation model failed to converge to a
solution for any of the simulated data sets with which I tested it. The mixed model’s
failure to converge likely was due to the highly complex nature of the model which
prevented the estimation of the random effects and associated variance parameters. The
random effects module of AD Model Builder has not been used to estimate variance
parameters in SCAA before to my knowledge. Additional experimentation with this
software may produce statistical catch-at-age mixed models with better convergence
properties. Monte Carlo simulations were used to investigate the performance of the

different methods.

Methods

I used a simulation study to evaluate different methods of estimating process and
observation error variances in SCAA. A data generating model was used to simulate data
sets from a hypothetical fish population. The estimation models, each using a different
error estimation method, were fit to the simulated data sets. The data generating model,
ad hoc and Bayesian estimation models were all built using AD Model Builder software
(Otter Research Limited 2004). For the following discussion, descriptions of all the
symbols are given in Table 2.1, while many of the equations describing my models are
given in Tables 2.2 and 2.3. I reference these equations as Equation x. y, where equation y

is found in Table x.

53

Data Generating Model

I developed a data generating model to simulate the dynamics of a hypothetical
fish population based on lake Whitefish stocks in the upper Great Lakes. The p0pulation
dynamics were described using abundance-at-age and age-specific mortality rates created
by the model. A gill net fishery operating on the population produced observed total
annual catch, age composition and fishing effort data.

I generated abundance-at-age using an exponential population function (Equation
2.2.1). To produce abundance-at-age in the first year (Equation 2.2.2), mortality was
applied to randomly generated numbers of age-1 fish, which were drawn from a
lognormal distribution (Table 2.4). The mean of the distribution was chosen by assuming
the population experienced equilibrium recruitment prior to the model time series.
Recruitment to the first age in subsequent years was calculated with a Ricker stock-
recruitment function (Equation 2.2.3; Table 2.4). The number of female spawners was
calculated as one-half of the number of fish age-3 and older, thereby assuming knife-edge
maturity and a 1:1 sex ratio.

Total mortality was partitioned into natural mortality and fishing mortality
sources (Equation 2.2.4). Natural mortality was a constant value for all years and ages
(Table 2.4). Fishing mortality was generated using a fully separable fishing mortality
model (Equation 2.2.5), where the age effect consisted of age-specific selectivity and the
year effect consisted of year-specific catchability and observed fishing effort. Age-
specific selectivity values were specified to create a dome-shaped selectivity curve,
which is typical of gill net fisheries (Table 2.4). Catchability varied from year to year

according to a lognormal white noise model (Table 2.4):

54

_ E

56M N N(O’°'§)°

The value for the standard deviation of log catchability 021 was randomly generated from

a lognormal distribution with three different means representing low, medium and high
levels of catchability variation (Table 2.4). An “observed” point estimate of the log
catchability standard deviation was generated, which simulated information that a stock
assessment analyst might possess. The observed point estimate was drawn from a

lognormal distribution:

(2) 0'

2
~ N 0,0" .
{q ( aq J
The log-scale standard deviation of the log catchability standard deviation 0' a q was

given the same value as was used to cause the true standard deviation of log catchability
to depart from its underlying median. Therefore, I effectively assumed that the observed
point estimate of the standard deviation of log catchability came from a lognormal
distribution with the same median as the true standard deviation of log catchability, but
with double the log-scale standard deviation as did the true standard deviation’s
generating distribution. Doubling the standard deviation represents the addition of
observation error to the measurement of the log catchability standard deviation. Fishing
effort was specified so that effort increases to a maximum in the middle of the time series
and then decreases to the end of the time series (Table 2.4). This fishing effort pattern

simulated a growing fishery that was regulated by effort limitations during the second

55

half of the time series. Total mortality Z0 used to produce abundance-at-age in the first

year (Equation 2.2.2) was generated with Equations 2.2.4 and 2.2.5 with the assumption
that fishing effort in years prior to the first year of the analysis was equal to fishing effort
in the first year of the analysis.

1 generated observed data from a gill net fishery from simulated abundance-at-age

and mortality rates. Catch-at-age was calculated using Baranov’s catch equation

(Equation 2.2.6). Observed total annual catch Cy was calculated by summing catch-at-

age Cyfi across ages for each year and incorporating observation error say:
_ :y
(3) Cy ’ Z Cy.a ‘3 ’
a=l

5C,y ~ N(0,ag).
I chose to use multiplicative lognormal errors because this is a standard assumption in

SCAA (Fournier and Archibald 1982; Deriso et al. 1985). The value for the standard

deviation of log total catch 0'C was randomly generated from a lognormal distribution

with two different means representing low and high levels of observation error (Table
2.4). An “observed” point estimate of the log total catch standard deviation was
generated, which simulated information that a stock assessment analyst likely would

possess. The observed estimate a'C was drawn from a lognormal distribution:

(4) 02: = ace“ ,
2
4C ~ N(0,ao_c )

56

The log-scale standard deviation of the log total catch standard deviation 0", C was given

the same value as was used to cause the true standard deviation of log total catch to
depart from its underlying median. Observed fishery age composition data was generated
by drawing a random sample from a multinomial distribution with a sample size of 100,
and proportions calculated from catch-at—age in the fishery (Equation 2.2.7). Any errors
in measuring fishing effort were lumped with interannual variation in catchability as

process error. Natural mortality was known without error.

Estimation Models

The estimation models used the same equations as the data generating model
except when estimating abundance-at-age in the first year, recruitment, and selectivity.
Annual recruitment was estimated as a mean recruitment parameter and a vector of
annual recruitment deviation parameters (i.e., a vector of deviations that must sum to
zero). Abundance-at-age in the first year was estimated as a mean abundance parameter
and a vector of abundance deviation parameters (i.e., a vector of deviations that must sum
to zero). Selectivity was estimated as a double logistic function of age (Equation 2.2.8).
Abundance-at-age (Equation 2.2.1), total mortality (Equation 2.2.4), fishing mortality
(Equation 2.2.5), catchability (Equation 2.1), catch-at-age (Equation 2.2.6), total catch
(Equation 3), and proportion of catch-at-age (Equation 2.2.7) were calculated as in the
data generating model. True parameter values produced by the data generating model
were used as starting values for parameters in the estimation models, to expedite
numerical searches during the simulations.

The estimation models differed from each other in the method used to estimate

variances for process error in catchability and observation error in total catch. First, an

57

ad hoc approach was used in which the proportion of process error variance was set so
that predicted observation error variance was consistent with an observed point estimate
of observation error variance (Richards et al. 1997). This approach has Bayesian aspects
because conditional on the value of the proportion of process error variance, point
estimates are obtained by maximizing the posterior density (Schnute 1994). Second, a
Bayesian approach with explicit priors on the variances was used in which the marginal
posterior densities of the variances were estimated. I considered two variants of the
Bayesian approach, one with an informative lognormal prior for log catchability variation

and the second with only a weakly informative lognormal prior for this variation.

Ad Hoc Estimation Model

In the ad hoc approach, I estimated the variances using a technique developed by
Richards et al. (1997). This approach requires repeated fits of the model with the

proportion of total variance due to log catchability variance p:

2
0'
(5) p=—%,
K
(6) K2=ag+ag,

being varied among fits. During each fit of the model, total variance was estimated as a

model parameter, and from this parameter the variances of log catchability a"; and log

total catch 03; were calculated as follows:

(7) 0.3 = pxz.

(& 03=0-pk2-

58

I varied the proportion of log catchability variance from 0.1 to 0.95 in increments of 0.05,

and refit the estimation model to a given data set for each value of p . I chose as best

among these model fits the one where the predicted standard deviation of log total catch
was closest to the “observed” point estimate of the log total catch standard deviation
created by the data generating model.

For a given model fit (specific p) using the ad hoc approach, highest posterior

density estimates of the parameters (a widely used approach, see Schnute 1994) were
obtained by maximizing the posterior density of the parameters conditional on the
observed data (Equations 2.3.1, 2.3.2a, and 2.3.3). I chose to minimize the negative log
posterior density (Equation 2.3.4) for ease of computation.

The probability density of the data conditional on the parameters was separated
into two components for total annual catch and proportion of catch-at-age (Equation
2.3.5). Total annual catch was assumed to follow a lognormal distribution, with the log

density (ignoring some additive constants) given by Equation 2.3.6. Proportion of catch—

at-age was assumed to follow a distribution that would arise if N E fish were observed,

with numbers observed at each age following a multinomial distribution, with the log
density (ignoring some additive constants) given by Equation 2.3.7. Note that the
probability density of the data conditional on the parameters is equivalent to the classical
likelihood function. Therefore, the highest posterior density parameter estimates are
equivalent to the maximum likelihood estimates.

The prior probability density of the parameters was separated into three

components for the general model parameters ¢ , catchability deviations sq, and total

59

variance K (Equation 2.3.8a). Devratrons 1n catchability were assumed to follow a

lognormal distribution, with the log prior density (ignoring some additive constants)

given by Equation 2.3.9. The prior densities of the log of all model parameters in ¢ and
2 . . . . . .
K were assrgned proper unrform prior densrtres, which follows common practice with the

intent of being weakly informative. Therefore, prior density of the log of a5 and K were

constants for all parameter values.

Bayesian Estimation Models

In the Bayesian approaches statistical inference was made on the posterior density
of the parameters conditional on the observed data (Equation 2.3.1) which was derived
using a Markov Chain Monte Carlo (MCMC) method. I chose to work with the negative
log posterior density for ease of computation (Equation 2.3.4). The standard deviations
of log-scale catchability and total catch were included as parameters to be estimated in
the model (Equation 2.3.2b). The probability density of the data conditional on the
parameters was separated into two components for total annual catch and proportion of
catch-at-age (Equation 2.3.5). The log densities for each of the components were the
same as in the ad hoc estimation model (Equations 2.3.6 and 2.3.7).

The prior probability density of the parameters was separated into four

components for the prior probability densities of the general model parameters ¢ ,

catchability deviations eq, log catchability standard deviation oq, and log total catch

standard deviation O'C (Equation 2.3.8b). Deviations in catchability were assumed to

follow a lognormal distribution as in the ad hoc estimation model (Equation 2.3.9). In

60

the first version of the full Bayesian approach, hereafier referred to as the informative
Bayesian approach, the standard deviations for log total catch and log catchability were
assumed to follow a lognormal distribution, with log prior density (ignoring some

additive constants) expressed as:

 

(9) ln[p(0'i)] = — 12 (Inc;- —ln0',-)2 —ln 0'01. ,
20'

0i
where i indexes the two error sources (i.e., total catch and catchability). The values for
the prior standard deviations for the standard deviations of log total catch and log
catchability were the same values used to create the true standard deviations in the data
generating model. The prior densities of the log of all general model parameters ¢were
assigned weakly informative proper uniform prior densities. Therefore, prior density of

the log of ¢ was a constant for all parameter values.

Marginal posterior densities for the standard deviations of total catch and
catchability were estimated using a MCMC method. The highest posterior density
parameter estimates served as starting values for the MCMC chain. A Metropolis-
Hastings algorithm with a scaled multivariate normal candidate generating distribution
was used to determine the marginal posterior densities (Gelman et al. 2004). The MCMC
chain was run for 500,000 cycles with values being saved every 25th cycle. The first
2,000 saved cycles of the MCMC chain were dropped as a bum-in period, in order to
remove the effect of the starting values (Gelman et al. 2004).

In reality, stock assessment analysts rarely have the data necessary to set such an
informative prior on the standard deviation of log catchability as I did in the informative

Bayesian estimation approach. Therefore, I also evaluated performance of the full

61

Bayesian method using a less informative prior. I refer to this as the objective Bayesian
approach. This approach was identical to the informative Bayesian approach except that
the prior density for the standard deviation of log catchability was assumed to follow a
lognormal distribution (Equation 9) with mean (0.35) and variance (0.49) specified so

that the prior density spanned all three levels of catchability variation.

Estimation Model Evaluation

My Monte Carlo simulation included six scenarios based on the three levels of
catchability variation and two levels of total catch variation. Five hundred data sets were
generated for each scenario for a total of 3,000 simulated data sets. Each estimation
model was fit to each of the simulated data sets. Estimation model runs were dropped
from the analysis if they exhibited poor convergence characteristics. After examining

preliminary results, ad hoc estimation model convergence was judged to be poor if the

. . -4 . . . .
maxrmum gradient component was greater than 1x10 . After examrmng preliminary

results, informative and objective Bayesian estimation model convergence was judged to
be poor if the effective sample size for log catchability standard deviation, log total catch
standard deviation, total abundance in the last year of analysis or highest posterior density
value was less than 200. Effective sample sizes were calculated from MCMC chains
using the method described by Thiebauz and Zwiers (1984) with lags out to 100 for
autocorrelation calculations.

The three approaches for estimating process and observation errors were
evaluated using the relative error (RE) of the standard deviations of log catchability,
standard deviation of log total catch and total abundance in the last year of the analysis.

The RE of the standard deviations of log catchability and log total catch indicated how

62

well the variances were estimated, while the RE of total abundance indicated how well
the approaches estimated a key management quantity. Relative error was calculated as

follows:

(10) RE: ,

 

where J? is a point estimate of the quantity of interest from the estimation model, and X
is the true value of the quantity of interest from the data generating model. For the
Bayesian methods I used the median of the marginal posterior distribution as a point
estimate, whereas for the ad hoc method the highest posterior density estimates were
used. The median of the relative errors (MRE) was used to examine systematic bias in
estimates from the estimation models. Median absolute relative error (MARE) , which
captures elements of bias and precision, was used to compare the range of relative errors

estimated by the estimation models.

Results

The following results are based on sample sizes of 500 model runs per scenario
for the ad hoc approach, 380 to 431 model runs per scenario for the informative Bayesian
approach, and 345 to 396 model runs per scenario for the objective Bayesian approach.
The number of poorly converged model runs for the informative and objective Bayesian
approaches is likely an artifact of my simulation study design. I had to limit the length of
the MCMC chains to reduce computational times and make the study feasible. Under
normal circumstances, an analyst would probably run longer MCMC chains or run

multiple chains from different starting points to improve convergence properties.

63

The informative Bayesian approach outperformed the ad hoc and objective
Bayesian approaches in the estimation of log total catch standard deviation. The
informative Bayesian approach was less biased than ad hoc and objective Bayesian
approaches in estimating standard deviation of log total catch (Figure 2.1). Informative
Bayesian approach MRE values for all six scenarios were close to zero and ranged from -
0.023 to 0.018. Objective Bayesian approach MRE values for all six scenarios exhibited
positive bias and ranged from 0.021 to 0.153. Ad hoc approach MRE values exhibited
negative bias and ranged from -0.338 to -0.019, except for the high catchability-low total
catch variance scenario (0.004).

Informative and objective Bayesian approaches demonstrated higher levels of
precision than the ad hoc approach in the estimation of log total catch standard deviation
(Figure 2.1). The differences in MARE values between informative Bayesian and ad hoc
approaches were small (-0.030 to 0.001) for medium catchability-low total catch, high
catchability-low total catch, and high catchability-high total catch variance scenarios
(Figure 2.2). The differences in MARE values between informative Bayesian and ad hoc
approaches were larger (-0.249 to -0.114) for low catchability-low total catch, low
catchability-high total catch, and medium catchability-high total catch variance scenarios
(Figure 2.2). The differences in MARE values between objective Bayesian and ad hoc
approaches were small (-0.074 to 0.050), except for the low catchability-high total catch
variance scenario (-0.246) (Figure 2.3).

The informative Bayesian approach also out performed the ad hoc and objective
Bayesian approaches in the estimation of the log catchability standard deviation. The

informative Bayesian approach was less biased than the ad hoc and objective Bayesian

64

approaches in estimating the standard deviation of log catchability (Figure 2.4).
Informative Bayesian approach MRE values for all six scenarios were close to zero, with
a small positive bias, ranging from 0.002 to 0.023. Objective Bayesian approach MRE
values generally were close to -l .0. The two exceptions were the medium catchability-
low total catch variance and high catchability-low total catch variance scenarios for
which the objective Bayesian approach MRE values were -0.048 and -0.064 respectively.
Ad hoc MRE values were negatively biased and ranged from -0.758 to -0.3 77.

The informative Bayesian approach was more precise than the ad hoc and
objective Bayesian approaches in estimating the standard deviation of log catchability
(Figure 2.4), although all methods had much lower precision for estimating the standard
deviation of catchability than for estimating the standard deviation of catch (note
difference in scale between Figure 2.1 and Figure 2.4). The differences in MARE values
between informative Bayesian and ad hoc approaches were substantial and ranged from -
0.637 to -0.292 (Figure 2.2), where the informative Bayesian approach was more precise.
The differences in MARE values between objective Bayesian and ad hoc approaches
generally were large and ranged from 0.235 to 0.548 (Figure 2.3), where the ad hoc
approach was more precise than the objective Bayesian approach. The two exceptions
were the medium catchability-low total catch variance and high catchability-low total
catch variance scenarios, -0.321 and -0.262 respectively, where the objective Bayesian
approach was more precise than the ad hoc approach.

Differences in performance between ad hoc, informative and objective Bayesian
approaches in the estimation of the total abundance in the last year of the analysis were

less marked than for variance estimates. For all three methods, the bias of the estimates

65

of total abundance in the last year increased at high catchability and total catch variance
levels (Figure 2.5). Ad hoc approach MRE values were negatively biased, and ranged
from -0.271 to —0.033. Informative Bayesian approach MRE values generally were close
to zero, positively biased, and ranged from 0.001 to 0.054. The one exception was the
high catchability-high total catch variance scenario which was 0.107. Objective Bayesian
approach MRE values generally were close to zero and ranged from -0.067 to 0.017. The
one exception was the high catchability-high total catch variance scenario which was -
0.124.

Precision of ad hoc, informative and objective Bayesian approach estimates of
total abundance in the last year decreased as catchability and total catch variance levels
increased (Figure 2.5). The Ad hoc approach was slightly less precise than the
informative Bayesian approach, with differences in MARE values ranging fiom -0.024 to
0.011 (Figure 2.2). Differences between objective Bayesian and ad hoc approach MARE

values were small and ranged from -0.015 to 0.059 (Figure 2.3).

Discussion

My results show that observation error variance will be more reliably estimated
than process error variance in SCAA. Observation error variance is better estimated due
to the availability of better prior information about observation errors. Estimates of
observation error variance obtained separately from SCAA, through analysis of the raw
data used to derive such quantities as observed total catch, provide a good source of prior
information for estimating observation error variance in SCAA. Such prior information
does not exist for process error variance because separate from SCAA the raw observed

data are not informative about how model parameters such as catchability vary. This was

66

demonstrated in my study when the ad hoc and objective Bayesian approaches produced
more accurate and precise estimates of log total catch standard deviation then of log
catchability standard deviation. These two approaches used more weakly informative
prior information or no prior information for log catchability standard deviation than for
log total catch standard deviation.

Use of the Bayesian approach allows for reliable estimation of both observation
and process error variances using a realistic, weakly informative prior for the process
error variance, when process error variability is greater than observation error variability.
Under this condition, the relatively strong informative prior for the observation error
variance and the strong signal for the process errors in the observed data allow SCAA to
reliably estimate the amount of total variance and successfully partition that variance
between observation and process error variances. In my study, this was evident when the
objective Bayesian approach was able to accurately estimate the log total catch and log
catchability standard deviations in scenarios where annual variability in catchability was
the dominant error source. Schnute and Richards (1995) found that, in general, their
catch-at-age estimation models performed better in a Monte Carlo simulation when
process error in recruitment was greater than observation error in an index of abundance.
Their estimation models estimated the process and observation error variances by
specifying the proportion of total variance due to process error variance, similar to the ad
[we approach, and obtaining maximum likelihood estimates of the variances analytically.
Unfortunately, Schnute and Richards (1995) did not look specifically at how their
estimation models performed at estimating the error variances. I hypothesize that process

error variability likely will be greater than observation error variability, and hence the

67

associated variances can be estimated, in any well monitored commercial fishery, as well
as most well monitored recreational fisheries. Chen and Paloheimo (1998) also have
suggested that errors due to environmental variation (i.e., process errors) may be greater
than observation errors for many fisheries. This finding emphasizes that another means
of improving the estimation of process and observation error variances, as well as stock
assessments in general, is to improve the quality of fishery monitoring data.

The ad hoc approach failed to reliably estimate the process and observation error
variances in my study. I was not surprised by this finding since the ad hoc approach
utilized the least amount of prior information (i.e., a single point estimate of log total
catch standard deviation) to estimate both of the standard deviations. More interesting
was the consistent underestimation of total variance in the system when using the ad hoc
approach. This negative bias might in part be explained by the statistical properties of the
estimator for the variances. Unlike the Bayesian approach which derived variance
estimates from the median of the posterior probability density, the ad hoc approach
simply used highest posterior density estimates to obtain variance estimates. Highest
posterior density parameter estimates share many similar properties with likelihood-based

parameter estimates, since highest posterior density estimates are obtained by
maximizing the probability density of the data given the parameters p(x|6) (Equation
2.3.1), which is identical to the likelihood fimction. Under this paradigm, the prior
probability densities p(6) could be thought of as penalty terms added to the likelihood

function. The maximum likelihood estimator of variance is known to be negatively
biased, thus the highest posterior density estimate of variance probably would possess the

same negative bias.

68

The ad hoc approach produced unbiased estimates of the log total catch standard
deviation in scenarios where catchability variation was greater than total catch variation,
but this apparent success is deceptive and potentially dangerous for analysts. Even in the
scenarios where catchability variation was dominant, the estimation model still
underestimated the total variance as evidenced by the associated negative bias in
estimates of the log catchability standard deviation. The estimation model was able to
match predicted and observed log total catch standard deviation values by adjusting the
proportion of total variance due to catchability variance, but the selected catchability
variance proportion did not match the true proportion from the data generating model. In
a real stock assessment where the true variances are unknown, such a result would lead
the analyst to believe that the total variance had been well estimated when, in fact, it had
been underestimated. This problem might be solved by correcting the predicted total
variance by the number of parameters estimated in the model, thus producing an unbiased
estimate of the total variance. Further study is needed to determine how well this total
variance correction would work, but it has the potential of making the ad hoc approach a
viable variance estimation technique.

I should point out that my study examined the ability of the ad hoc approach to
estimate one form of process error (i.e., catchability variation). The only other published
use of the ad hoc approach was to produce estimates of recruitment variability in a state-
space age-structured model, but the approach was applied to actual fishery data and its
performance was not quantified nor evaluated (Richards et al. 1997). The ad hoc
approach can be classified with other methods that use residual model error to estimate

associated variances, because the ad hoc approach employs the measured interannual

69

variation in the observed data as prior information for the stock assessment model
estimates of the process and observation error variances. As another example of this
class of methods, Francis et al. (2003) compared standard specified values of commercial
CPUE and trawl survey CVs to resulting residual variation between observed values and
stock assessment predictions of CPUE. This approach could be applied in an iterative
method to obtain variance estimates. An initial variance value would be specified and the
stock assessment model fit to the observed data. The resulting residual variation in
model results would be used to specify a new variance value for the next model run. This
process would be repeated until the specified variance value matched the resulting
residual variation in model results. The assessment models used for lake Whitefish in
1836 treaty waters have used a such an iterative approach to setting the variance
associated with variability about an assumed stock-recruitment relationship (Ebener et al.
2005). The Francis et al. (2003) study examined actual data from New Zealand fisheries,
and the Whitefish assessments use actual data also, so it is unknown how accurately
residual variation in stock assessment model results measures the true underlying
variance. Further study of the ability of these residual-based variance estimation
approaches to estimate other forms of process error variability, such as time-varying
selectivity and annual recruitment variations, would be useful and informative.

The ad hoc and Bayesian approaches performed equally well at estimating
numbers of fish in the last year of the analysis, even though the ad hoc approach
consistently underestimated the process and observation error variances. In theory, the
poor performance of the ad hoc approach in estimating the variances should have resulted

in poorer estimates of the final number of fish. To address this issue, it is necessary to

70

consider when it is important to properly estimate the error variances. Methot (1990)
suggested that when the different data sources used in SCAA do not trend over time, or
when trends are consistent between data sources, assessment model results will be less
sensitive to changes in the variances which are used to weight the different data sources.
It is when trends in the different data sources are inconsistent with each other that
assessment model results will be sensitive to changes in the variance values (Methot
1990). Therefore, it is most important to properly estimate the error variances when the
data sources are sending mixed signals about the population dynamics to the assessment
model. In my study, total catch did trend over time, but catchability did not since the
catchability deviations were generated using a white noise model. IfI had generated
catchability so that it trended over time, then it is likely that I would have seen
differences in the estimation model performances when it came to estimating the final
number of fish. Such an analysis was beyond the scope of this study, since I wanted to
determine if it were possible to estimate the error variances under the simplest conditions
I could imagine. Actual stock assessments are generally more complex, incorporating
multiple sources of observation and process error. As a result, I feel it would be
informative to evaluate the Bayesian and ad hoc approaches when estimating more than
two sources of variation.

I recommend that stock assessment analysts use the Bayesian approach when
attempting to estimate process and observation error variances in SCAA. The Bayesian
approach is fairly robust when existing data allow for the designation of strongly
informative priors for the error variances, particularly process error variance. The

Bayesian approach still can produce reliable estimates of the error variances with a

71

weakly informative prior for the process error variance, as long as high quality
monitoring data are available. I do not recommend the use of the ad hoc approach based
on my findings. The ad hoc approach consistently underestimates the error variances,
which could lead to biased estimates of important management quantities when the
different data sources send inconsistent signals concerning the dynamics of the

population.

72

References

Chen, Y., and J. E. Paloheimo. 1998. Can a more realistic model error structure improve
the parameter estimation in modeling the dynamics of fish populations? Fisheries
Research 38: 9-17.

Deriso, R. B., T. J. Quinn II, and P. R. Neal. 1985. Catch-at-age analysis with auxiliary
information. Canadian Journal of Fisheries and Aquatic Sciences 42: 815-824.

Ebener, M. P., J. R. Bence, K. Newman, and P. Schneeberger. 2005. Application of
statistical catch-at-age models to assess lake Whitefish stocks in the 1836 treaty-
ceded waters of the upper Great Lakes. Pages 271-309 in L. C. Mohr and T. F.
Nalepa, editors. Proceedings of a workshop on the dynamics of lake Whitefish
(Coregonus clupeaformis) and the amphipod Diporeia spp. in the Great Lakes.
Great Lakes Fishery Commission, Technical Report 66, Ann Arbor.

Fournier, D., and C. P. Archibald. 1982. A general theory for analyzing catch at age
data. Canadian Journal of Fisheries and Aquatic Sciences 39: 1195-1207.

Francis, R. I. C. C., R. J. Hurst, and J. A. Renwick. 2003. Quantifying annual variation
in catchability for commercial and research fishing. Fishery Bulletin 101: 293-
304.

Gelman, A., J. B. Carlin, H. S. Stern, and D. B. Rubin. 2004. Bayesian Data Analysis,
2nd edition. Chapman and Hall/CRC, Boca Raton, Florida.

Law, A. M., and W. D. Kelton. 1982. Simulation Modeling and Analysis. McGraw-
Hill, New York.

Maunder, M. N. 2001. A general framework for integrating the standardization of catch
per unit of effort into stock assessment models. Canadian Journal of Fisheries and
Aquatic Sciences 58: 795-803.

Merritt, M. F., T. J. Quinn, 11. 2000. Using perception of data accuracy and empirical
weighting of information: assessment of a recreational fish population. Canadian
Journal of Fisheries and Aquatic Sciences 57: 1459-1469.

Methot, R. D. 1990. Synthesis model: an adaptable framework for analysis of diverse
stock assessment data. International North Pacific Fisheries Commission Bulletin
50: 259-277.

National Research Council. 1998. Improving Fish Stock Assessments. National
Academy Press, Washington, DC.

73

Otter Research Limited. 2004. An introduction to AD Model Builder version 7.0.1 for
use in nonlinear modeling and statistics. Otter Research Limited, Sidney, British
Columbia.

Otter Research Limited. 2005. Random effects in AD Model Builder: ADMB-RE user
guide. Otter Research Limited, Sidney, British Columbia.

Quinn 11, T. J ., and R. B. Deriso. 1999. Quantitative Fish Dynamics. Oxford University
Press, New York.

Richards, L. J ., J. T. Schnute, and N. Olsen. 1997. Visualizing catch-age analysis: a case
study. Canadian Journal of Fisheries and Aquatic Sciences 54: 1646-1658.

Ricker, W. E. 1975. Computation and interpretation of biological statistics of fish
populations. Fisheries Research Board of Canada Bulletin 191.

Schnute, J .T. 1994. A general framework for developing sequential fisheries models.
Canadian Journal of Fisheries and Aquatic Sciences 51: 1676-1688.

Schnute, J. T., and L. J. Richards. 1995. The influence of error on population estimates
from catch-age models. Canadian Journal of Fisheries and Aquatic Sciences 52:
2063-2077. ‘

Sitar, S. P., J. R. Bence, J. E. Johnson, M. P. Ebener, and W. W. Taylor. 1999. Lake
trout mortality and abundance in southern Lake Huron. North American Journal
of Fisheries Management 19: 881-900.

Thiebauz, H. J ., and F. W. Zwiers. 1984. The interpretation and estimation of effective
sample size. Journal of Climate and Applied Meteorology 23: 800-811.

74

Table 2.1. Symbols and descriptions of variables used in data generating and estimation

 

 

models.
Symbols Description Application
Cyfl Number of fish caught by year and age Both
5 y Observed number of fish caught by year Both
Ey Fishery effort by year Both .
F J40 Instantaneous fishing mortality by year and age Both
M Instantaneous natural mortality Both
Nyfl Abundance by year and age Both
N0 Mean abundance for abundance in first year Estimation
N E Number of fish used to calculate age composition each year Both
Pyfl Proportion of catch by year and age Both
FLO Observed proportion of catch by year and age Both
R0 Mean recruitment Estimation
Sy Number of female spawners by year Generation
Zyfl Instantaneous total mortality by year and age Both
20,0 Instantaneous total mortality for abundance in first year by age Generation
b1 First inflection point of double logistic selectivity function EStimation
b2 First slope of double logistic selectivity funcion Estimation
b3 Second inflection point of double logistic selectivity function Estimation
b4 Second slope of double logistic selectivity funcion Estimation
m Total number of ages Both
n Total number of years Both
10(9 Ix ) Posterior probability density of parameters conditional on data Estimation
p(x|6) Probability density of data conditional on parameters Estimation
p(6) Prior probability density of parameters Estimation
qy Fishery catchability by year Both
(7 Median catchability Both
sa Fishery selectivity by age Both
0: Productivity parameter of Ricker recruitment function Generation

 

75

Table 2.1 (cont’d).

 

3

8R};

5w
3C.y

¢

2
K

#N
9

p

V’a

CC

Density dependent parameter of Ricker recruitment function
Process error in recruitment by year
Process error in catchability by year

Observation error in number of fish caught by year
Subset of model parameters common to both estimation
models

Total variance
Mean age-1 abundance for abundance in first year
Set of all model parameters

Proportion of total variance due to catchability variance
Standard deviation of age-l abundance for abundance in first
year

Standard deviation of log-scale recruitment

Standard deviation of log-scale catchability

Observed point estimate of log catchability standard deviation
Log-scale standard deviation of log catchability standard
deviation

Standard deviation of log-scale total catch

Observed point estimate of log total catch standard deviation
Log-scale standard deviation of log total catch standard
deviation

Process error in recruitment by year

Process error for abundance in first year by age
Observation error in log catchability standard deviation

Observation error in log total catch standard deviation

Generation
Generation
Both
Both

Estimation
Estimation
Generation
Estimation

Estimation

Generation

Generation
Both

Both

Generation
Both

Both
Generation

Estimation

Estimation
Generation

Generation

 

76

Table 2.2. Data generating and estimation model equations.

 

 

 

Equations Application
2.2.1 —Z Both
Ny+l,a+l = Ny,ae y,a
2.2.2 a—l Generation
- 220.}
N“, = N2_a,le 1:1 ;fora >1
2.2.3 Ny,l =aSy_1e“flSy_le£R’y §5R,y ~ N(0,o,2¢) Generation
2.2.4 Zy,a = M + Fy,a Both
2.2.5 Fy,a = Saquy Both
2.2.6 F -Z Both
_ y’a _ ysa
Cy” — Z (I e )Ny a
y,a
2.2.7 Cy 0 Both
Pya = ’
Cy
2.2.8 Estimation

 

 

l I
= 1-
S“ 1+e—b2(a-b1){ 1+e-b4(a-b3)]

 

77

Table 2.3. Posteriormrobability density equations for estimation models.

 

 

Equations Application
2-3-1 p(9|JC)°<= p(x|9)p(9) Both
2.3.2a e = if [5w L1 , K} Ad hoc
2.3.2b t9 = klqukzl , 0g 0C} Bayesian
2'3'3 ¢= {No [Vim ,-___ 1 ”KO [01]” ,-_ .16] b1 b2 b3 b4} BOth
2-3-4 -1n[p(6|x)loc-ln[p(x|6)l-lnlp(6)l Both
2-3-5 maxim]:nnmcwimwlai B...
2.3.6 ~ Both
ln[p(C|6)]=— 271C y: :l[(ln C y — In C y )2] — n In 0C
y= a=
2.3.83 ln[p(6)] = ln[p ¢)]+ lnlpieq )J+ ln[p(K)] Ad hOC
2.3.8b 1n[p(6)]=In[p(¢)]+lnlp(eq)J+Inlpla, )J+1n[p(oc)] Bayesian
2.3.9 ln[p(£q)] $2:qu y] ’11an Both

 

78

Table 2.4. Values of variables used in data generating model to create simulated data
sets.

 

 

Variable Level Value
n 20
m 8
#N 355000
0N 0.4
a 10.1
fl 5.10E-06
0R 0.4
M 0.24
E 0.1, 2.0, 3.0, 3.1, 3.3, 3.7, 4.4, 5.3, 6.5, 8.0, 8.0, 6.5, 5.3, 4.4, 3.7,
y 3.3, 3.1, 3.0, 2.0, 0.1
Sa 0.04, 0.15, 0.43, 0.85, 1.00, 0.82, 0.57, 0.37
(7 0.15
Eq Low 0.2
Medium 0.5
High 0.8
oq 0.2
a-‘C Low 0.25
High 0.75
o a C . 0.2
NE 100

 

79

 

 

 

 

 

 

 

1-0 ' 1:: Ad hoc
Informative Bayesian

0 5 — Objective Bayesian o
E .
if:
g o.o~
5
6'2

-O.5 -

'1.0 ' fl I I I in
Total catch: low high low high low high
Catchability: low medium high

Scenafio

Figure 2.1. Box plots representing relative error distributions for estimates of log total
catch standard deviation across different levels of catchability and total catch variance.
The bars represent median relative errors. The boxes, whiskers, and circles represent
25th and 75th, 10th and 90th, and 5th and 95th percentiles of the distributions,

respectively.

80

.0
is

0 Log total catch SD
0 Log catchability SD
0 Numbers last year

.0
N

 

 

 

 

 

(D
0
E, 0.0 ----~o---°---9---<>---e---°----
(D o O
a:
5 -02 - o
< '0.4 ‘
2 D
-0.6 - a U
'0.8 I I I I I I I
Total catch: low high low high low high
Catchability: low medium high
Scenano

Figure 2.2. Differences in median absolute relative errors (MARE) between informative
Bayesian approach and ad hoc approach across different levels of catchability and total
catch variance. Symbols represent informative Bayesian approach MARE values minus
ad hoc approach MARE values.

81

 

 

 

 

 

1-0 ' 0 Log total catch SD
0 8 . 0 Log catchability SD
' 0 Numbers last year

8 05 ..
C ' n
9
g 0.4 " D
O n
Lu (12‘- U
SE _,, o
E 0.0 '_——'8_-_°'-—'8—— o ——e ———————

412‘- o D

D

‘0.4 I I I I I r u
Total catch: low high low high low hig
Catchability: low medium high

Scenano

Figure 2.3. Differences in median absolute relative errors (MARE) between the objective
Bayesian approach and the ad hoc approach across different levels of catchability and
total catch variance. Symbols represent objective Bayesian approach MARE values
minus ad hoc approach MARE values.

82

1-0 ' E Ad hoc
Informative Bayesian
— Objective Bayesian

0.5 -

Relative error

 

 

-10 - '0- o-

 

l T l I I

Total catch: low high low high low high
Catchability: low medium high
Scenano
Figure 2.4. Box plots representing relative error distributions for estimates of log
catchability standard deviation across different levels of catchability and total catch

variance. The bars represent median relative errors. The boxes, whiskers, and circles
represent 25th and 75th, 10th and 90th, and 5th and 95th percentiles of the distributions,

respectively.

83

 

 

 

 

 

 

3 i :1 Ad hoc
Informative Bayesian °
2 - Objective Bayesian
L " O
8
5
a)
.2
E
a)
o:
-1 , a , , , '
Total catch: low high low high low high
Catchability: low medium high

Scenano

Figure 2.5. Box plots representing relative error distributions for estimates of total
abundance in the last year of analysis across different levels of catchability and total
catch variance. The bars represent median relative errors. The boxes, whiskers, and
circles represent 25th and 75th, 10th and 90th, and 5th and 95th percentiles of the
distributions, respectively.

84

CHAPTER 3
EVALUATING AND SELECTING METHODS FOR ESTIMATING TIME-VARYING

SELECTIVITY IN STATISTICAL CATCH-AT—AGE ANALYSIS

Introduction

Statistical catch-at-age analysis (SCAA) is a common method of fisheries stock
assessment. Age-structured catch data from a fishery are used to estimate quantities of
interest, such as population abundance and mortality rates, using likelihood methods
(Fournier and Archibald 1982). Auxiliary data that provide an index of abundance either
directly or indirectly, such as survey catch-per-unit-effort (CPE) or fishery effort, are
essential for reliable estimation (Deriso et al. 1985; Methot 1990). Estimated population
quantities from the last year of the analysis are typically used as a starting point for short-
term projections that are the basis for recommending harvest limits or targets.

In many SCAA models fishing mortality is assumed to be separable into year and
age effects, with their product being the fishing mortality rate for a given year and age
(Doubleday 1976). Here I refer to the year effect as fishing intensity, and to the age
effect as fishery selectivity. Selectivity refers to the relative vulnerability of specific ages
of fish to a fishery, so that age classes that are highly selected tend to be overrepresented
in the catch in comparison to their relative abundance in the population. Selectivity is
influenced by fishing gear characteristics, as well as fishing and fish behavior.

Selectivity often is modeled either as a function of age or it is allowed to vary
freely among ages. The parameters of the selectivity function or the selectivity values for

each age are estimated within SCAA along with other model parameters. Logistic

85

(Millar 1995; Punt et al. 2001), double logistic (Methot 1990; Ebener et al. 2005),
exponential-logistic (Thompson 1994), normal (Millar 1995), lognormal (Millar 1995),
gamma (Deriso et al. 1985; Millar 1995), and polynomials (F oumier 1983) are some of
the functions used to model selectivity. Regardless of how selectivity is modeled, a
restriction often must be applied to ensure a unique parameterization of the age and year
effects (Doubleday 1976). Selectivity functions generally are constrained by normalizing
the function to a reference age or to the age of maximum estimated selectivity. When
selectivity is allowed to vary freely with age, selectivity commonly is constrained by
setting selectivity at some reference age(s) equal to one.

The separability assumption can be relaxed, allowing selectivity to change over
time, when there is evidence to suggest that selectivity is not constant (i.e., gear
characteristics or fish behavior have changed). Separate selectivity values can be
estimated for different blocks of time within SCAA (Radomski et al. 2005). Some of the
selectivity function’s parameters can vary over time independently from year to year
(Bence and Rogers 1993), according to a polynomial in time (Ebener et al. 2005) or
random walk process (Gudmundsson 1994; Ianelli 1996). Nonadditive models have been
used to allow selectivity to vary with changes in fishing mortality (Myers and Quinn
2002; Radomski et a1. 2005).

Statistical catch-at-age analysis has been shown to be sensitive to the choice of
how selectivity is modeled. Incorrect assumptions about selectivity have been shown to
generate errors in SCAA estimates of biomass (Kimura 1990), spawning biomass (Punt et

al. 2002; Radomski et al. 2005), exploitation rate (Radomski et al. 2005), and the ratio of

86

stock biomass in the first year to the stock biomass in the final year of analysis (Yin and
Sampson 2004).

Radomski et al. (2005) looked specifically at how specification of time-varying
selectivity affected SCAA. They compared three methods for estimating selectivity:
constant, time-blocked and nonadditive and found no one method for estimating time-
selectivity performed best in all situations, but they did discover that time-varying
selectivity SCAA models performed as well as constant selectivity SCAA models when
selectivity was constant, and outperformed constant selectivity SCAA models when
selectivity varied with time. They speculated that allowing selectivity to vary according
to a random walk might improve the estimation of time-varying selectivity (Radomski et
al. 2005). Radomski et al. (2005) also recommended that research was needed to
determine the extent to which correct or adequate selectivity models could be identified.

The objective of my study was to compare the performance of different time-
varying selectivity functions within SCAA. In addition, I strove to identify a model
selection method that could allow analysts to select among alternative time-varying
selectivity functions within a specific SCAA. This contrasts with an objective of
determining a single “best” time-varying selectivity estimation method, which works well
in most situations. Of course, one possible outcome of my work could have been that an
omnibus procedure for modeling selectivity works better than selecting among
alternatives. I addressed my objectives through Monte Carlo simulations, in which 1
evaluated different methods of both modeling time-varying selectivity within a stock

assessment and selecting among the methods.

87

Methods

I used Monte Carlo simulations to compare four time-varying selectivity
estimation methods and evaluate three model selection techniques. I used a data
generating model to simulate data sets from a hypothetical fish population. The data
generating model used two different approaches to simulate time-varying selectivity: 1) a
double logistic function in which the first inflection point varied according to a first order
autoregressive process, and 2) selectivity for each age varied independently according to
a first order autoregressive process. I chose these two approaches to provide contrast in
how freely selectivity varies over time. The double logistic function is constrained so
that only selectivity of younger age fish changes over time. The age-specific selectivity
parameters allow selectivity to vary more freely, with age-specific selectivity values
changing independently of each other. I fit four estimation models, each using a different
time-varying selectivity estimation method, to the simulated data sets. The selectivity
estimation methods consisted of 1) a double logistic function in which the first inflection
point varied according to a random walk, 2) a double logistic function in which the first
and second inflection points varied according to random walks, 3) a double logistic
function in which all four parameters varied according to random walks, and 4)
selectivity for each age varied according to a random walk with a smoothing function
across ages. I chose these estimation approaches because they represent the two general
approaches for estimating selectivity in SCAA, namely modeling selectivity as a function
of age and estimating age-specific selectivity parameters. In addition, these four
estimation approaches form a continuum of increasing flexibility in how selectivity is

allowed to vary over time. The three model selection techniques included 1) root mean

88

square error (RMSE), 2) degree of retrospectivity, and 3) the Deviance Information
Criterion (DIC). The data generating model and four estimation models were all built
using AD Model Builder software (Otter Research Limited 2004). For the following
discussion, descriptions of all the symbols are given in Table 3.1, while most of the
equations describing my models are given in Tables 3.2 and 3.3. I reference equations as
Equation x. y, where equation y is found in Table x.

My Monte Carlo simulation included two scenarios based on two different
methods for generating time-varying selectivity. Five hundred data sets were generated
for each scenario for a total of 1,000 simulated data sets. Each of the four estimation
models was fit to each of the simulated data sets. I applied the three model selection

techniques to each estimation model fit to a simulated data set.

Data Generating Model

I developed a data generating model to simulate the dynamics of a hypothetical
fish population based on lake Whitefish stocks in the upper Great Lakes. The population
dynamics were described using abundance-at-age and age-specific mortality rates created
by the model. A gill net fishery operating on the population produced observed total
annual catch, age composition and fishing effort data. Each simulated data set included
20 years of data for fish ages 1 to 8+, where 8+ is a plus group containing all fish age-8
and older.

1 generated abundance-at-age using an exponential population equation (Equation
3.2.1). In order to produce abundance-at-age in the first year, mortality was applied to
randomly generated numbers of age-1 fish (Equation 3.2.2). The number of age-1 fish

was randomly drawn from a lognormal distribution (Table 3.4). I selected the mean of

89

the distribution by assuming the population experienced equilibrium recruitment prior to
the simulated time series. I calculated recruitment to the first age in each year with a
Ricker stock-recruitment function (Equation 3.2.3; Table 3.4). I calculated the number of
female spawners as one-half of the number of fish age-3 and older, thereby assuming
knife-edge maturity and a 1:1 sex ratio.

I partitioned total mortality into natural and fishing mortality components
(Equation 3.2.4). Natural mortality was a constant value for all years and ages (Table
3.4). I modeled fishing mortality by relaxing the assumption of full separability
(Equation 3.2.5). I generated year and age-specific selectivity using two different
methods to create a dome-shaped selectivity curve, which is typical of gill net fisheries. I
defined fishing intensity as a function of fishing effort (Equation 3.2.6). The errors
associated with fishing intensity were a combination of process error due to annual
variation in catchability and observation error in nominal fishing effort. I assumed that
variation in catchability would outweigh observation error in fishing effort and, therefore,

treat the fishing intensity errors as process error. The value for the standard deviation of
log fishing intensity 0,1 was randomly generated from a lognormal distribution for each

simulated data set (Table 3.4). I specified fishing effort so that effort increased to a
maximum in the middle of the time series and then decreased to the end of the time series
(Table 3.4). This fishing effort pattern simulated a growing fishery during the first half
of the time series that was regulated by effort limits during the second half of the time

series. 1 generated the total mortality used to produce abundance-at-age in the first year

Z0 (Equation 3.2.2) using Equations 3.2.4, 3.2.5 and 3.2.6 with the assumption that

fishing effort in years prior to the first year of the analysis was equal to fishing effort in

90

the first year of the analysis, and selectivity in years prior to the first year of the analysis
was constant at the initial values.

The two methods I chose to generate time-varying selectivity provide contrast in
how selectivity changes over time. For the first method, I generated selectivity using a

double logistic function of age (Methot 1990):

1 l
s — L ‘ 1— .
(1) a 1+e—b2la‘bl,yll: 1+e—b4(a-b3)]

 

 

I varied the first inflection point over time from an initial value according to a first order

autoregressive process (Table 3.4): 6y ~ N (0, 0;)
(2) loge bl,y+1 = loge bl + pl (loge bl,y '— loge bl l+ 6y ,
2
5}, ~ N (0, 0'5 ).
I randomly drew the initial value of the first inflection point from a lognormal

distribution with mean 51 and standard deviation 05 The value for the standard

deviation of log first inflection point 05 was randomly generated from a lognormal

distribution for each simulated data set (Table 3.4). I normalized age-specific selectivity
in a given year using the maximum generated age-specific selectivity value for that year.
By allowing the first inflection point to vary over time, I was simulating a scenario in
which the vulnerability of young fish to the fishery was changing over time. For the
second method, I chose a more flexible approach to generating time-varying selectivity
based on a method used by Butterworth et al. (2003). In this approach, age-specific
selectivity varied over time from initial values according to a first order autoregressive

process (Table 3.4):

91

(3) loge Sy+l,a = loge s2, + p2 (loge 3),", —loge s;)+ y)“, ,

2
7y,a ~ N(O"7r l
I used the same correlation and standard deviation parameters for all ages. I randomly

drew the initial values for selectivity at each age from lognormal distributions with means

Ea and standard deviation 0,. The value for the standard deviation of log selectivity 0'),

was randomly generated from a lognormal distribution for each simulated data set (Table
3.4). I normalized age-specific selectivity in a given year using the maximum generated
age-specific selectivity value for that year. By allowing age-specific selectivity values to
vary over time, I simulated a scenario in which the vulnerability of all age classes of fish
to the fishery changed independently over time.

I generated observed data from a gill net fishery from simulated abundance-at-age
and mortality rates. I calculated catch-at-age using Baranov’s catch equation (Equation
3.2.7). I calculated observed total annual catch by summing catch-at-age across ages for

each year and incorporating observation errors (Equation 3.2.8; Table 3.4). The value
for the standard deviation of log total catch 0", was randomly generated from a lognormal

distribution for each simulated data set (Table 3.4). I generated observed fishery age
composition data by drawing a random sample from a multinomial distribution with a
sample size of 400, and proportions calculated from catch-at-age in the fishery (Equation
3.2.9). Natural mortality and observed fishing effort were known without error (Table

3.4).

92

Estimation Models

The estimation models used the same equations as the data generating model
except when estimating abundance-at-age in the first year, recruitment, and selectivity. I
estimated annual recruitment as a mean recruitment parameter and a vector of annual
recruitment deviation parameters (i.e., a vector of deviations that must sum to zero). I
estimated abundance-at-age in the first year as a mean abundance parameter and a vector
of abundance deviation parameters (i.e., a vector of deviations that must sum to zero). I
calculated abundance-at-age (Equation 3.2.1), total mortality (Equation 3.2.4), fishing
mortality (Equation 3.2.5), fishing intensity (Equation 3.2.6), catch-at-age (Equation
3.2.7), total catch (Equation 3.2.8), and proportion of catch-at-age (Equation 3.2.9) using
the equations described for the data generating model. I used true parameter values
produced by the data generating model as starting values for parameters in the estimation
models to expedite numerical convergence during simulations.

The estimation models differed from each other in the method used to estimate
time-varying selectivity for the fishery. The four methods I chose represent increasing
flexibility in the estimation of time-varying selectivity. The cost associated with
increased flexibility is an increase in the number of parameters that must be estimated. In
the first estimation approach, I allowed the first inflection point of the double logistic

function to vary over time according to a random walk:

 

 

1 [ l :l
sa = _ , _ a 1— __ _ ,
(4) 1+e bZl“ bl,yl l+e b4(a b3)

2
logewa =loge b1,y +ny, 77y ~ N(0,0',7 ).

93

This approach is the least flexible of those I examined since it only changes the lower
ages at which selectivity increases most rapidly over time. In the second estimation
approach, I allowed the first and second inflection points of the double logistic function

to vary over time according to random walks:

1 1
5 = , . 1_ . . ’
( ) 5“ 'bzla‘blyllr 1+e‘b4la‘53,yl]

1+e

 

 

loge bi,y+l = loge bi,y +77i,y , 771;)» ~ N(0,03),
where i indexes the inflection points of the double logistic function (i.e., b1.y and b3’y). I

made the simplifying assumption that the standard deviations of the two log-scale
inflection points were equal. This approach of varying the two inflection points allows
the ascending and descending limbs of the selectivity curve to expand and contract over
the course of time. In the third estimation approach, I allowed the two inflection points
and the two slopes of the double logistic function to vary over time according to random

walks:

 

 

1 I
(6) s = . . 1— . . ,
a 1+e—b2’y(a—bl’y) 1+e—b4’y(a—b3’y)

2
loge bi,y+l = loge bi,y + "Ly, 77i,y ~ N(0,0',7 )9
10 b- =10 b- +r- r- ~N(0 02)
3e J,y+1 ge by 1:)” by r r r
where i indexes the inflection points and j indexes the slopes (i.e., bz’y and b4'y) of the

double logistic function. Again, as I did for the infection points, I assumed that the
standard deviations of the two log-scale slopes were equal. This approach of allowing all

of the double logistic function parameters to vary over time provides maximum flexibility

94

in the estimation of time-varying selectivity for this functional form. In the fourth
estimation approach, I allowed age specific selectivity values to vary over time according

to random walks (Butterworth 2003):

(7) loge Sy+l,a =loge Sy,a +wy,av

to)” ~ N(0,a§, ).
I made the simplifying assumption that the year-specific standard deviations of log
selectivity were equal for all years. I constrained age-specific selectivity with a curvature
penalty using squared third-differences to ensure smoothness in selectivity across age

classes (Butterworth 2003):

 

‘ 1-31 I
(8) g(S 5y a;0_ 02): i1," Z] 3We y, a+3- Oges y, ,2a+:.0-+3 Oges y, a+l loges y, a)2 .
y]: a]: 9’

I made the simplifying assumption that the age-specific standard deviation of log-scale
selectivity was the same for all ages. I added this curvature penalty term to the negative
log posterior density. In all four time-varying selectivity estimation approaches, 1
normalized age-specific selectivity using the maximum estimated age-specific selectivity
value. I estimated the variances associated with log total catch, log fishing intensity and
log selectivity using a Bayesian approach in which the marginal posterior densities were
estimated with Markov Chain Monte Carlo simulations.

I made statistical inference on the posterior density of the parameters conditional
on the observed data (Equation 3.3.1) which was derived using a Markov Chain Monte
Carlo (MCMC) method. More specifically, I used MCMC with the Metropolis-Hastings
algorithm as it is implemented in AD Model Builder (Otter Research Limited 2004).

Maximum likelihood parameter estimates were used as starting values for each MCMC

95

chain. I ran the MCMC chain for each model for 1,000,000 cycles, saving parameter
values every 10th cycle. I dropped the first 40,000 cycles from the chain of saved
MCMC values as a burn in period, which reduced the effect of chain starting values on
final MCMC estimates (Gelman et al. 2004). I dropped model runs with poor
convergence properties from the analysis. I judged MCMC chain convergence to be poor
if the effective sample size for the highest posterior density value was less than 300. I
selected the highest posterior density value because it provides an overall measure of how
the MCMC chains are mixing. Effective sample sizes were calculated from MCMC
chains using the method described by Thiebauz and Zwiers (1984) with lags out to 100
for autocorrelation calculations. I chose to minimize the negative log posterior density
(Equation 3.3.2a) for ease of computation. For the fourth estimation approach in which
age-specific selectivity values varied over time, I added the curvature penalty term
(Equation 8) to my negative log posterior density (Equation 3.3.2b). The parameters
estimated in the model (Equation 3.3.3) included the subset of parameters common to all
of the estimation models and the subset of time-varying selectivity parameters at specific
to each estimation model.

The subset of parameters used to model time-varying selectivity depended upon
the method used to estimate selectivity. For the first estimation approach in which the
first inflection point of the double logistic function varied with time, the selectivity
parameters included the first inflection point in the first year, annual deviations in the
first inflection point, standard deviation of the log-scale first inflection point, and the
other three parameters of the double logistic function (Equation 3.3.4a). For the second

estimation approach in which both inflection points of the double logistic function varied

96

with time, the second inflection point was replaced by a second inflection point in the
first year, annual deviations in the second inflection point, and a standard deviation of the
log-scale second inflection point selectivity deviations (Equation 3.3.4b). For the third
estimation approach in which all four parameters of the double logistic function varied
with time, both slopes were also replaced by corresponding slopes in the first year, annual
deviations for each of these parameters (Equation 3.3.4c). For the fourth estimation
approach in which the age-specific selectivity values varied with time, the selectivity
parameters included the age-specific selectivity values in the first year, annual deviations
for each age-specific selectivity value, and standard deviations for the year and age-
specific log selectivity values (Equation 3.3.4d) ,

I separated the probability density of the data conditional on the parameters into
two components for total catch and proportion of catch-at-age (Equation 3.3.5). I
assumed total annual catch followed a lognormal distribution, with the log density

(ignoring some additive constants) given by Equation 3.3.6. I assumed proportion of

catch-at-age followed a distribution that would arise if N E fish were observed, with

numbers observed at each age following a multinomial distribution, with the log density
(ignoring some additive constants) expressed by Equation 3.3.7.

For all of the time-varying selectivity estimation approaches, I assumed the prior
probability densities of the random walk deviations for selectivity parameters followed
lognormal distributions, with the log prior densities (ignoring some additive constants)

expressed as:

<9) ln[p(z.-)]=--1—5§[z§yl-(n-l)ln0i,

97

where i indexes the time-varying selectivity parameters (e. g., first inflection point of
double logistic function). I assumed the prior probability densities of the log total catch,
log catchability, and log selectivity standard deviations followed lognormal distributions,

with the log prior densities (ignoring some additive constants) expressed as:

(10) 1n[p(a,.)]= —-—17(1na,'- 4110,)2 —1n.9,-,
2.9,-

where i indexes the error sources (e. g., observation errors in total catch). I assigned a
strong informative prior density (i.e., identical to the generating distribution from the data
generating model) to the log total catch standard deviation (Table 3.5). Thus, I assumed
the analyst had good prior information on how observation errors in total catch were
distributed, which is a reasonable assumption for a well monitored commercial fishery. I
assumed the analyst would not have such strong prior information for the other standard
deviations. Therefore, I assigned more weakly informative prior densities which allowed
the remaining standard deviations to vary over a realistic range of values (Table 3.5).

The time-varying age-specific selectivity parameter estimation model failed to converge

to a solution when weakly informative prior densities were assigned to the year and age-

specific log selectivity standard deviations, am and 0‘0 respectively. As a result, I fixed

the values for the year and age-specific log selectivity standard deviations at 0.15 and
0.08 respectively for all simulations. This solution followed the common practice of
assuming variances to be known when they cannot be estimated in the estimation model.
I assigned weakly informative uniform prior densities to the logs of all other model

parameters. Therefore, prior densities for each log-scale model parameter, excluding the

98

selectivity random walk deviations and their associated variances, were constants for all
parameter values.

I compared the performance of the four estimation models by calculating the
relative error (RE) of population biomass and exploitation rate in the last year of the

analysis, for each simulated data set:
(11) RE = — ,

where If is the point estimate of the quantity of interest from the estimation model, and
X is the true value of the quantity of interest from the data generating model. I used the
median of the marginal posterior distribution as a point estimate. Estimated biomass and
exploitation rate in the last year often play an important role when stock assessment
results are used to inform management actions. In addition, I used the median of the
relative errors (MRE) to examine whether there was systematic bias in estimates from the
estimation models. I used the median absolute relative error (MARE), which captures
elements of bias and precision, to compare the range of relative errors made when using

the estimation models.

Model Selection Methods

I evaluated the performance of three model selection techniques to determine
which technique(s), if any, could identify consistently the “best” time-varying selectivity
estimation approach. The three model selection techniques I used to identify the best
time-varying selectivity estimation approach were RMSE, degree of retrospectivity, and
DIC. By best selectivity estimation approach, I mean the estimation approach which

most closely predicts the true fish population as produced by the data generating model.

99

More specifically, for each simulated data set I measured relative performance of the
different estimation models based on the RE of population biomass and exploitation rate
in the last year of the analysis. I used three definitions of the best or nearly best
estimation model(s) for a given simulation run: 1) the estimation model producing the
lowest final population biomass or exploitation rate RE, 2) estimation models producing
RES within 0.05 of the lowest RE, and 3) estimation model producing RES within 0.1 of
the lowest RE. I allowed for this relaxation in the definition of best or nearly best
estimation model because in a real stock assessment, where the true population
characteristics are unknown, alternative estimation models which produce similar results
often would be treated as equally viable. In particular, I chose the values 0.05 and 0.1
because they represented a difference in model results that most analysts would consider
negligible. In addition, I used the MRE and MARE to examine bias and precision in
estimates from the estimation models chosen by each selection method. Comparison of
the model selection methods was made using the subset of simulation runs in which all
four estimation models converged on adequate solutions to avoid problems with different
convergence rates between the estimation models.

My first model selection procedure focused on proportion of catch-at-age
residuals, with the selected model minimizing the RMSE for these residuals. I chose this
as one possible method because I thought generally large pr0portion of catch-at-age
residuals might occur for estimation models that incorrectly modeled selectivity patterns.
These residuals were calculated from the posterior medians of predicted proportions of

catch—at-age.

100

My second model selection method is based on retrospective analysis, which
involves the comparison of successive estimates of model output quantities as additional
years of data are added to the stock assessment (Parma 1993; Mohn 1999). For this
selection method I selected the model that minimized the absolute value of Mohn’s

(1999) degree of retrospectivity statistic:

X(1:y),y "thany

 

(12) DR =

9

y=n-10 X(1:y),y

where X (1 ,y), y is the predicted value of quantity X in year y estimated from the data set

spanning year 1 to year y and X (1 _. n y y is the predicted value of quantity X in year y

estimated from the data set spanning year 1 to the last year of the full data set n. Here I
conducted a retrospective analysis for each estimation model-simulated data set fit by
dropping a year of data from the simulated data set and refitting the estimation model,
repeating this process until the last 10 years of data had been sequentially removed from
the analysis. Systematic retrospective patterns in model quantities can occur when time-
varying processes are modeled as being constant over time (Mohn 1999). Though all of
my estimation models allowed selectivity to change over time, I expected to see
retrospective patterns in cases where an estimation model had difficulty tracking changes
in selectivity. To make this approach practical, I used highest posterior density estimates
of the parameters with the variance parameters fixed at their point estimates from the
analysis of the full data set.

My final selection method was to select the model that minimized the Deviance
Information Criterion (Spiegelhalter et al. 2002). Information-theoretic model selection

criteria generally work by balancing model goodness of fit against model complexity

101

(i.e., the number of parameters in the model). The effective number of parameters in
complex models, such as my SCAA models, is often less than the actual number of
parameters due to various constraints placed upon those parameters. I chose to use DIC,
as opposed to the more commonly used Akaike’s Information Criterion (AIC; Akaike
1973) and Bayesian Information Criterion (BIC; Schwartz 1978), because DIC provides a
means of estimating the effective number of parameters. Wilberg (2005) demonstrated in
a different SCAA application that selection by DIC could result in estimates with lower
mean square errors, than always using any particular single model.

Deviance Information Criterion is composed of two components (Spiegelhalter et

al. 2002):

(13) DIC=I§+pD,
where D is the average deviance and p D is the effective number of parameters. I

estimated the average deviance as (Spiegelhalter et al. 2002):

_ C
(14) D =é—Z—2log. p(xl6c),

0:]

where C is the number of MCMC cycles saved minus the burn in and p(x|t9c) is the

probability of data x conditional on parameters 6 from MCMC cycle c. I estimated the

effective number of parameters as (Wilberg 2005):
(15) PD =D—‘D(9ML)a
where D( 9M1.) is the deviance evaluated at the maximum likelihood parameter estimates

and the other variables are described above.

102

For each model selection method, I calculated in what percentage of the
simulation runs the method selected the best or nearly best estimation models. In
addition, for each selection method I examined the distribution of RES for final
population biomass and exploitation rate estimates.

1 compared the performance of using estimation models selected by degree of
retrospectivity to the performance of always using the same estimation model, for each of
the estimation models. The objective here was to determine if this model selection
technique outperformed the omnibus approach of always using the same estimation
model. I used model selection by degree of retrospectivity in this evaluation because of
the good performance of this model selection method (see Results). To properly make
this comparison, I used degree of retrospectivity to select the best estimation based on
final p0pulation biomass and exploitation rate for each Simulation run, rather than for the
subset of Simulation runs where all four estimation models converged on solutions.
Comparisons were made using MRE and MARE values for final population biomass and
exploitation rate selected by degree of retrospectivity and estimated by each of the

estimation models.

Results

Model runs exhibiting poor convergence characteristics were dropped from the
analysis. The following results are based on sample sizes of 333 to 425 model runs per
scenario (Table 3.6). All of the dropped model runs failed to converge to highest
posterior density solutions, thus MCMC simulations could not be run. I suspect that with
sufficient effort, which would be warranted for a real assessment, an analyst could have

made adjustments in many of these cases to achieve convergence. This was not practical

103

in the context of this simulation study. It should be noted that the subsets of simulation
runs demonstrating poor convergence characteristics generally were different for each of
the estimation models (i.e., incidents of poor convergence were not due to characteristics

of particular simulated data sets).

Estimation Models

There was little difference in the biases of the four estimation models’ estimates
of population biomass in the last year of analysis within each data generating scenario
(Figure 3.1). The four estimation models produced less biased estimates of the final
population biomass in the double logistic generating scenario compared to the age-
Specific selectivity parameters generating scenario. Median relative error values for final
population biomass ranged from 0.01 to 0.13 for the double logistic function generating
scenario (Table 3.6). In contrast, MRE values for final population biomass ranged from -
0.23 to 0.55 for the age-specific selectivity parameters generating scenario. The
estimation model using the double logistic function with four time-varying parameters
produced the most biased estimates of population biomass in both data generating
scenarios.

The four estimation models produced more precise estimates of population
biomass in the last year of analysis when the estimation models more accurately
represented the true underlying population (i.e., when selectivity estimation and data
generating models were similar) (Figure 3.1). Median absolute relative error values for
final population biomass varied from 0.20 to 0.26 for the three double logistic function
estimation models in the double logistic function generating scenario (Table 3.6). On the

other hand, the age-specific selectivity parameters estimation model had a final

104

population biomass MARE value of 0.54 for the double logistic function generating
scenario. The age-specific selectivity parameters estimation model had final population
biomass MARE value of 0.35 for the age-specific selectivity parameters generating
function. In contrast, the three double logistic function estimation models had final
population biomass MARE values ranging from 0.50 to 0.61 for the age-specific
selectivity parameters generating scenario. The estimation model using the double
logistic function with four time-varying parameters produced estimates of final
population biomass that were less precise than the estimation models using double
logistic functions with one and two time-varying parameters (Figure 3.1).

There was little difference in the biases of the four estimation models’ estimates
of exploitation rate in the last year of analysis within each data generating scenario
(Figure 3.2). The four estimation models produced less biased estimates of the final
exploitation rate in the double logistic generating scenario compared to the age-specific
selectivity parameters generating scenario. Median relative error values for final
exploitation rate ranged from -0.10 to -0.02 for all four of the estimation models in the
double logistic function generating scenario (Table 3.6). In contrast, the MRE values for
exploitation rate ranged from -0.36 to -0.18 for all four estimation models in the age-
specific selectivity parameters generating scenario. The estimation model using the
double logistic function with four time-varying parameters produced the most biased
estimates of population biomass in both data generating scenarios.

The four estimation models produced more precise estimates of exploitation rate
in the last year of analysis when the estimation models more accurately represented the

true underlying population, though the difference was not as pronounced in the age-

105

specific selectivity parameters generating scenario (Figure 3.2). Median absolute relative
error values for the final exploitation rate varied from 0.20 to 0.25 for the three double
logistic function estimation models in the double logistic generating scenario (Table 3.6).
On the other hand, the age-specific selectivity parameters estimation model had an
exploitation rate MARE value of 0.58 for the double logistic generating scenario. The
age-Specific selectivity parameters estimation model had a final exploitation rate MARE
value of 0.38 for the age-specific selectivity parameters generating function. In contrast,
the three double logistic function estimation models had final exploitation rate MARE
values ranging from 0.48 to 0.57 for the age-specific selectivity parameters generating

scenario.

Model Selection

I compared the performance of the model selection methods by examining the
subset of simulation runs where all four of the estimation models exhibited good
convergence properties. All of the estimation models converged on good solutions for
438 of the 1,000 Simulation runs.

Degree of retrospectivity selected the best or nearly best estimation model, based
on final population biomass and exploitation rate RES, as often as or more often than DIC
and RMSE (Figure 3.3). Degree of retrospectivity selected the best or nearly best model
in 34—57% of the simulation runs when the best or nearly best model was chosen based on
final population biomass RE, and in 33-52% of the simulation runs based on final
exploitation rate RE. Deviance information criterion selected the best or nearly best
model in 27-48% of the simulation runs when the best or nearly best model was chosen

based on final population biomass RE, and in 27-50% of the Simulation runs based on

106

final exploitation rate RE. Root mean square error selected the best or nearly best model
in 29-49% of the simulation runs when the best or nearly best model was chosen based on
final population biomass RE, and in 33-49% of the simulation runs based on final
exploitation rate RE.

Selecting estimation models using degree of retrospectivity produced estimates of
population biomass and exploitation rate in the last year of analysis that were as biased
and precise as or less biased and more precise than estimation models selected using DIC
and RMSE (Figures 3.5 and 3.6). In particular, degree of retrospectivity selected
estimation models that produced final population biomass and exploitation rate estimates
that were less biased and more precise than estimates selected by DIC and RMSE in the
age-specific selectivity parameters generating scenario (Table 3.7).

Degree of retrospectivity performed as well as or better than the individual
estimation models at estimating final population biomass and exploitation rate. Degree
of retrospecitivity produced a final population biomass MRE of 0.05 and MARE of 0.24
in the double logistic generating scenario, which is comparable to the estimation
performances of the three time-varying double logistic functions in that same scenario
(Table 3.6). Degree of retrospecitivity produced a final population biomass MRE of -
0.01 and MARE of 0.40 in the age-specific selectivity parameters generating scenario,
which is less biased than any of the individual estimation models and of intermediate
precision between the age-specific selectivity parameters and double logistic function
estimation models in that same scenario (Table 3.6). Degree of retrospecitivity produced
a final exploitation rate MRE of -0.05 and MARE of 0.24 in the double logistic

generating scenario, which is comparable to the estimation performances of the three

107

time-varying double logistic functions in that same scenario (Table 3.6). Degree of
retrospecitivity produced a final exploitation rate MRE of 0.03 and MARE of 0.42 in the
age-specific selectivity parameters generating scenario, which is less biased than any of
the individual estimation models and of intermediate precision between the age-specific
selectivity parameters and double logistic function estimation models in that same

scenario (Table 3.6).

Discussion

There was no single time-varying selectivity estimation model that outperformed
the others in all situations that I examined. Rather, the estimation model(s) that produced
the estimates most tightly distributed about true population biomass and exploitation rate
in the last year of analysis was the one that most closely represented the true underlying
population. The three estimation models that used variants of the double logistic function
to model time-varying selectivity produced better estimates of final population biomass
and exploitation rate than the age-specific selectivity parameters estimation model when
the selectivity of the true population was generated with a double logistic function.
Likewise, the age-specific selectivity parameters estimation model produced better
estimates of final population biomass and exploitation rate than the three double logistic
function estimation models when the selectivity of the true population was generated with
age-specific selectivity parameters. This sort of result is common to simulation studies
where there are similarities between data generating and estimation models (e. g.,
Radomski et al. 2005; Wilberg and Bence 2006).

My study suggests that if an analyst knows the underlying form that selectivity

takes in a fish population, then he or she can model time-varying-selectivity reasonably

108

well. This begs the question, how well can an analyst know true selectivity patterns and
how they vary over time? Many fishing gears such as gill nets and trap nets are size
selective. If one makes the assumption that the size and age of fish are correlated and by
extension selectivity and age of fish are correlated, then modeling selectivity as some
function of ages, which produces a smooth selectivity curve, is a reasonable approach. It
is more difficult to think of situations where age-specific selectivity values vary relatively
independently of each other over time. Kimura (1990) demonstrated that estimating age-
specific selectivity parameters outperformed the use of a selectivity function when the
function was incorrectly Specified. The approachiof Butterworth et al. (2003) that I used
in this study is an extension of Kimura’s (1990) approach, which allows the age-Specific
selectivity parameters to vary over time. Further study is needed to determine whether
estimation models that assume time-varying age-specific selectivity parameters
outperform a time-varying selectivity function when the function is misspecified.

Model complexity is another issue that must be addressed when evaluating
different time-varying selectivity models. Increased model complexity means an
increased number of parameters that must be estimated, which can lead to over-
pararneterization of the model. An over-pararneterized model can produce poor
parameter estimates with high variances (Bumham and Anderson 2002). In my study,
the issue of model complexity was most clearly demonstrated in the performance of the
double logistic function with four time-varying parameters and two associated variances.
I expected the four time-varying parameter selectivity estimation method to outperform
the other double logistic function approaches when the observed data were generated

using age-specific selectivity parameters, due to the increased flexibility granted by

109

allowing all four double logistic function parameters to vary over time. Instead, I found
that the four time-varying parameter double logistic function produced more biased
estimates of final population biomass than the other double logistic function estimation
approaches when the observed data were generated using age-Specific selectivity
parameters. One of the two variances associated with the slopes and inflection points of
the four time-varying parameter double logistic function was estimated as nearly zero
(i.e., making it effectively the same as the two time-varying parameter double logistic
function) in many of the simulation runs, which suggests that the observed data were not
informative enough to estimate all of the selectivity parameters.

The performance of the four time-varying parameter double logistic function in
my study could be due to my data generating model design. The observed data were
generated by allowing selectivity parameters to vary over time according to a first order
autoregressive process, which did not follow any trend over time, and for which the
deviations among ages were not correlated. The performance of the time-varying double
logistic methods, may have improved had the generating selectivity function produced
correlated changes in selectivity for adjacent ages, like those that would be generated by
variations in one or more parameters of a firnction.

I was surprised to see how well degree of retrospectivity performed as a time-
varying selectivity model selection method compared to DIC and RMSE. Estimation
models selected using degree of retrospectivity produced final population biomass and
exploitation rate estimates that were more or equally accurate and precise compared to
estimates selected by DIC and RMSE for the data generating scenarios I examined. In

particular, degree of retrospectivity selected final p0pulation biomass estimates that were

110

much more accurate and precise than those estimates selected by DIC and RMSE when
observed data sets were generated using age-specific selectivity parameters. The
robustness of degree of retrospectivity as a time-varying selectivity estimation model
selection method probably is due to the fact that it can detect consistent patterns in model
estimates over time (i.e., as new years of data are added to the model). As I expected,
these consistent or retrospective patterns do appear to be indicative of an estimation
model that has difficulty estimating time-varying selectivity. Deviance information
criterion and RMSE lack this ability to detect retrospective patterns since they merely
evaluate the model fit to the complete time series of observed data. Parma (1993)
developed an alternative metric for identifying retrospective patterns using the square
root of the mean square error between the retrospective estimate of a model quantity and
a corresponding reference estimate on the log scale. Mohn (1999) points out that this
mean square error metric is unable to differentiate between retrospective and random
patterns since it uses a mean square, rather than signed sum, in its calculation. Though I
did not test Parma’s (1993) metric, I suspect that it would perform similarly to my DIC
and RMSE methods. I recommend that degree of retrospectivity be used to select
between estimation models using different methods of estimating time-varying
selectivity, based on its performance in my study.

Selecting from multiple estimation models using degree of retrospectivity worked
better than choosing a single estimation model in my study. Nothing is lost in estimation
performance by using degree of retrospectivity, even if an analyst is able to correctly
specify time-varying selectivity. In addition, degree of retrospectivity outperforms

estimation models which misspecifiy time-varying selectivity (i.e., assuming a double

1]]

logisitic function when true age-specific selectivity values vary over time). Therefore, I
recommend that degree of retrospectivity be used to select between time-varying
selectivity models.

I should note that my study only looked at the performance of model selection
methods on an individual basis. The ability to select the best estimation model may be
improved by using combinations of different selection techniques. For example,
estimation models could be ranked based on their degree of retrospectivity. If multiple
estimation models have equal or nearly equal degree of retrospecitvity values, then DIC
or RMSE values could be used to select between those models with degree of
retrospecitvity values close to zero. Further study of using such multiple model selection

methods would be informative.

112

References

Akaike, H. 1973. lnforrnation theory as an extension of the maximum likelihood
principle. Pages 267-281 in B.N. Petrov and F. Csaki, editors. Second
International Symposium on Information Theory. Budapest: Akademiai Kiado.

Bence, J .R. and J .E. Hightower. 1993. Status of Bocaccio in the
Conception/Monterey/Eureka INPFC areas in 1992 and recommendations for
management in 1993. Appendix B in: Status of the Pacific coast groundfish
fishery through 1992 and recommended acceptable biological catches for 1993.
Pacific Fishery Management Council, Portland, OR.

Bumharn, K.P., and DR. Anderson. 2002. Model Selection and Multimodel Inference: a
Practical Information-Theoretic Approach. Springer-Verlag, New York.

Butterworth, D.S., J .N. Ianelli, and R. Hilbom. 2003. A statistical model for stock
assessment of southern bluefin tuna with temporal changes in selectivity. African
Journal of Marine Science 25: 331-361.

Deriso, R.B., T.J. Quinn II, and RR. Neal. 1985. Catch-at-age analysis with auxiliary
information. Canadian Journal of Fisheries and Aquatic Sciences 42:815-824.

Doubleday, W.G. 1976. A least squares approach to analyzing catch at age data.
Research Bulletin of the International Commission of Northwest Atlantic
Fisheries 12:69-81.

Ebener, M. P., J. R. Bence, K. Newman, and P. Schneeberger. 2005. Application of
statistical catch-at-age models to assess lake Whitefish stocks in the 1836 treaty-
ceded waters of the upper Great Lakes. Pages 271-309 in LG Mohr and TR
Nalepa, editors. Proceedings of a workshop on the dynamics of lake Whitefish
(Coregonus clupeaformis) and the amphipod Diporeia spp. in the Great Lakes.
Great Lakes Fishery Commission Technical Report 66.

Fournier, D. 1983. An analysis of the Hecate Strait Pacific cod fishery using an age-
structured model incorporating density-dependent effects. Canadian Journal of
Fisheries and Aquatic Sciences 40: 1233-1243.

Fournier, D., and GP. Archibald. 1982. A general theory for analyzing catch at age data.
Canadian Journal of Fisheries and Aquatic Sciences 39:1195-1207.

Gelman, A., J .B. Carlin, H.S. Stern, and DB. Rubin. 2004. Bayesian Data Analysis, 2nd
edition. Chapman and Hall, Boca Raton.

Gudmundsson, G. 1994. Time series analysis of catch-at-age observations. Applied
Statistics 43:117-126.

113

Ianelli, J. 1996. An alternative stock assessment model of the eastern Bering Sea pollock
fishery. Pages 81-104 in Stock assessment and fishery evaluation report for the
groundfish resources of the Bering Sea/Aleutian Islands regions, November 1996.
North Pacific Fisheries Management Council, Anchorage, Alaska.

Kimura, OK. 1990. Approaches to age-structured separable sequential population
analysis. Canadian Journal of Fisheries and Aquatic Sciences 47:2364-2374.

Methot, RD. 1990. Synthesis model: an adaptable framework for analysis of diverse
stock assessment data. Pages 259-277 in L. Low, editor. Proceedings of the
symposium on applications of stock assessment techniques to gadids.
International North Pacific Fisheries Commission Bulletin 50.

Millar, RB. 1995. The functional form of hook and gillnet selection curves cannot be
determined from comparative catch data alone. Canadian Journal of Fisheries and
Aquatic Sciences 52:883-891.

Mohn, R. 1999. The retrospective problem in sequential population analysis: an
investigation using cod fishery and simulated data. ICES Journal of Marine
Science 56:473-488.

Myers, R.A., and T.J. Quinn II. 2002. Estimating and testing non-additivity in fishing
mortality: implications for detecting a fisheries collapse. Canadian Journal of
Fisheries and Aquatic Sciences 59:597-601.

Parma, A.N. 1993. Retrospective catch-at-age analysis of pacific halibut: implications
on assessment of harvesting policies. Pages 247-265 in G. Kruse, D.M. Eggers,
C. Pautzke, R.J. Marasco, and T.J. Quinn H, editors. Proceedings of the
International Symposium on Management Strategies for Exploited Fish
Populations. Alaska Sea Grant College Program.

Punt, A. E., D. C. Smith, R. B. Thomson, M. Haddon, X. He, and J. M. Lyle. 2001.
Stock assessment of the blue grenadier Macruronus novaezelandiae resource off
south-eastern Australia. Marine and Freshwater Resources 52:701-717.

Punt, A.E., A.D.M. Smith, and G. Cui. 2002. Evaluation of management tools for
Australia’s South East Fishery: 2. how well can management quantities be
estimated? Marine and Freshwater Research 53:631-644.

Radomski, P., J .R. Bence, and T.J. Quinn 11. 2005. Comparison of virtual population
analysis and statistical kill-at-age analysis for a recreational kill-dominated

fishery. Canadian Journal of Fisheries and Aquatic Sciences 62:436-452.

Schwartz, G. 1978. Estimating the dimension of a model. Annals of Statistics 6: 461-
464.

114

Spiegelhalter, D.J., N.G. Best, B.P. Carlin, and A. van der Linde. 2002. Bayesian
measures of model complexity and fit. Journal of the Royal Statistical Society
Series B 63: 583-639.

Thiebauz, H. J ., and F. W. Zwiers. 1984. The interpretation and estimation of effective
sample size. Journal of Climate and Applied Meteorology 23: 800-811.

Thompson, G.G. 1994. Confounding of gear selectivity and the natural mortality rate in
cases where the former is a nonmonotone function of age. Canadian Journal of
Fisheries and Aquatic Sciences 51:2654—2664.

Wilberg, M.J. 2005. Improving statistical catch-at-age stock assessments. Doctoral
Dissertation. Michigan State University, East Lansing, Michigan.

Wilberg, M.J., and J .R. Bence. 2006. Performance of time-varying catchability
estimators in statistical catch-at-age analysis. Canadian Journal of Fisheries and
Aquatic Sciences 63: 2275-2285.

Yin, Y., and DB. Sampson. 2004. Bias and precision of estimates from an age-

structured stock assessment program in relation to stock and data characteristics.
North American Journal of Fisheries Management 24:865-879.

115

Table 3.1. Symbols and descriptions of variables used in data generating and estimation
models.

 

Symbol Description Application
Cyfl Number of fish caught by year and age Both

6:y Observed number of fish caught by year Both

Ey Fishery effort by year Both

F yfl Instantaneous fishing mortality by year and age Both

M Instantaneous natural mortality Both

Nyfl Abundance by year and age Both

N0 Mean abundance for abundance in first year Estimation
N E Number of fish used to calculate age composition each year Both

Pyfl Proportion of catch by year and age Both

PM, Observed proportion of catch by year and age Both

R0 Mean recruitment Estimation
S), Number of female spawners by year Generation
Zyfl Instantaneous total mortality by year and age Both

ZQa Instantaneous total mortality for abundance in first year by age Generation
b1, y First inflection pt. of double logistic selectivity function by year Both

b2 First slope of double logistic selectivity funcion Both

b3 Second inflection pt. of double logistic selectivity function Both

b4 Second slope of double logistic selectivity funcion Both

bi Mean of first inflection pt. of double logistic selectivity function Estimation
13, Fishing intensity by year Both

m Total number of ages Both

n Total number of years Both

116

Table 3.1 (cont’d).

 

19891")
phW)
19(9)
q

Posterior probability density of parameters conditional on data
Probability density of data conditional on parameters

Prior probability density of parameters

Fishery catchability

Fishery selectivity by year and age

Mean fishery selectivity by age

Mean fish weight by age

Productivity parameter of Ricker recruitment function

Density dependent parameter of Ricker recruitment function
Process error in selectivity parameter i by year

Process error in first inflection point of double logistic function
by year
Process error in recruitment by year

Subset of time-varying selectivity parameters
Process error in selectivity by year and age

Process error in inflection points of double logistic function by
year
Error in fishing intensity by year

Mean number of age-1 fish for abundance in first year
Set of all model parameters

Prior standard deviation of log-scale fishing intensity standard
deviation

Prior standard deviation of log-scale inflection points standard
deviation

Prior standard deviation of log-scale total catch standard
deviation

Prior standard deviation of log-scale slopes standard deviation

First correlation parameter for first order autoregressive process

Estimation
Estimation
Estimation
Both

Both

Generation
Both

Generation
Generation
Estimation

Generation

Generation
Estimation
Generation

Estimation

Both
Generation
Estimation

Estimation
Estimation
Estimation

Estimation

Generation

 

117

Table 3.1 (cont’d).

 

P2

0N

Second correlation parameter for first order autoregressive
process

Standard deviation of number of age-1 fish for abundance in
first year

Standard deviation of log-scale first inflection point

Generating mean of log-scale first inflection point standard
deviation
Standard deviation of log-scale recruitment

Standard deviation of log-scale selectivity

Generating mean of log-scale selectivity standard deviation
Standard deviation of log-scale inflection points

Prior mean of log-scale inflection points standard deviation
Age-specific standard deviation of log-scale selectivity
Standard deviation of log-scale fishing intensity

Generating and prior mean of log-scale fishing intensity
standard deviation
Standard deviation of log-scale slopes

Prior mean of log-scale slopes standard deviation
Standard deviation of log-scale total catch

Generating and prior mean of log-scale total catch standard
deviation
Year-specific standard deviation of log-scale selectivity

Process error in slopes of double logistic function by year
Observation error in number of fish caught by year
Process error in selectivity by year and age

Process error in recruitment by year

Process error for abundance in first year by age

Generating standard deviation of log-scale first inflection point
standard deviation

Generation

Generation

Generation

Generation

Generation
Generation
Generation
Estimation
Estimation
Estimation
Both

Both

Estimation
Estimation
Both
Both

Estimation
Estimation
Both

Estimation
Estimation
Estimation

Generation

 

118

Table 3.1 (cont’d).

 

C Generating standard deviation of log-scale selectivity standard Generation
7 deviation
C). Generating standard deviation of log-scale fishing intensity Generation
standard deviation
4V Generating standard deviation of log-scale total catch standard Generation

deviation

119

Table 3.2. Data generating and estimation model equations.

 

 

 

Equation Application
3.2.1 —2 Both
Ny+l,a+1 = Ny,ae y,a
3.2.2 a—l Generation
_ 220,].
Nl,a = N2—a,le 1:]
3.2.3 -,BS _ e 2) Generation
Ny’l = y_]e y ‘e as, ~N(0,a,,
3.2.4 Zyfl = M + Fyfl Both
3.2.5 Fyfl = Sy,afy Both
3.2.6 fy : que/ly ,4); ~ N(0,oi‘) Both
3.2.7 F _ Both
Cy,a = _y,a Ny’a(1 —e Zy’a)
Zyfl
3.2.8 ~ U Both
Cy = ZCJ’fl e y,uy ~ N(0,03)
a
3.2.9 C Both
P, a = dd
y

 

120

Table 3.3. Posterior probability density equations for estimation models.

Equation

 

3.3.1 p(6|x) oc p(x|t9)p(0)
3.3.22: —1n[p(6lx)]oc«Labial—mime]
3.3.2b —1n[p(t9|x) + 8(3y,a :0; )0: -1n[P(x|9)]— ln[p(6)]+ g(sy,a;a%)

3.3.3 6 = holy/013:1,R0,[a)y};=29qr¢}
3.3.4a ¢ = {bl’l’[,7y}"y:1,o-n,b2,b3,b4l
3.3.4b ¢_ {bu [my]:l 1,,b2 193,, {my}: ' b4 0 a}

334C ¢= $11 [Ulyk= 11’ [)2], [T2y};:_ll 9,9173] [7,3, y;ll}ny ”b4,1 [14,)»:3’ 11,0",01}

3.3.4d ¢ ={Phallrflwy’a};]:=1’0”"08}
3-3-5 maxim]: 1nLv(6I6)1+1nLa(PI6)1

3.3.6
1n[p(C|6)]=—n[(ln Cy - In C y )2 ] — n In 0C
2;C y: _1

3.3.7

maple): ZN $103.. 1. P...)

 

121

Table 3.4. Values of quantities used in data generating model to create simulation data
sets.

 

Quantity Value

 

n 20

m 8

.UN 355,000

O'N 0.4

a 10.1

fl 5.10E-06

as 0.4

[Wa 123:1 0.20, 0.48, 0.73, 0.91, 1.32, 1.52, 1.76, 2.15

M 0.24

[E }n 0.8, 1.6, 2.4, 3.2, 4.0, 4.8, 5.6, 6.4, 7.2, 8.0, 8.0, 7.2, 6.4, 5.6, 4.8, 4.0, 3.2,
y y=l 2.4, 1.6, 0.8

q 0.15

0:1 0.4

C}. 0.1

bf 4.01

[bi]?=2 1.40, 3.49, 0.50

[3; 31:] 0.04, 0.15, 0.43, 0.85, 1.00, 0.82, 0.57, 0.37

p] 0.9
02; 0.2
45 0.1
p2 0.9
0’7 0.15
4, 0.1
0;, 0.05
(V 0.1
NE 400

 

122

Table 3.5. Values used to define prior probability densities in estimation models.

 

Quantity Value

 

0', 0.4
.9 1. 0.25
0;, 0.2
3,, 0.42
0;, 0.05
.9, 0.1
a; 0.2
.9, 0.42

 

123

Table 3.6. Median relative errors (MRE), median absolute relative errors (MARE), and
number of replicates (N) for estimates of final population biomass and exploitation rate
produced by the time-varying selectivity estimation models: double logistic functions
with one (DLl), two (DL2), and four (DL4) time-varying parameters, and time-varying
age-specific selectivity parameters (ASP).

Population Biomass

 

 

Estimation DLl ASP

Model MRE MARE N MRE MARE N
DLl 0.05 0.20 414 0.23 0.50 411
DL2 0.06 0.22 361 0.33 0.57 425
DL4 0.13 0.26 333 0.55 0.61 430
ASP 0.01 0.54 382 -0.23 0.35 409

Exploitation Rate

DLl -0.04 0.20 414 -O.l8 0.50 411
DL2 -0.06 0.21 361 -0.24 0.51 425
DL4 -0.10 0.25 333 -0.36 0.48 430
ASP -0.02 0.58 382 0.30 0.38 409

 

124

Table 3.7. Median relative errors (MRE), median absolute relative errors (MARE), and
number of replicates (N) for estimates of final population biomass and exploitation rate
chosen by the model selection methods: root mean square error (RMSE), deviance
information criterion (DIC), and degree of retrospectivity (DR).

Population Biomass

 

 

 

 

Model DLl ASP

Selection MRE MARE N MRE MARE N
RMSE 0.10 0.21 179 0.43 0.53 259
DIC 0.10 0.22 179 0.40 0.51 259
DR 0.07 0.22 179 -0.05 0.35 259

Exploitation Rate

RMSE -0.08 0.21 179 -0.32 0.54 259
DIC -0.09 0.22 179 -0.29 0.50 259
DR -0.06 0.24 179 0.10 0.37 259

 

 

125

3000 - [:1 DL1 o o
1500 - m DL2
/— DL4

Relative Error

 

 

 

.1}.
1+»
11

 

 

Data generating model

Figure 3.1. Box plots representing relative error distributions for estimates of population
biomass in the last year of analysis across different data generating models. The data
generating and estimation models include double logistic functions with one (DLl), two
(DL2), and four (DL4) time-varying parameters, and time-varying age-specific selectivity
parameters (ASP). The bars represent median relative errors. The boxes, whiskers, and
circles represent 25th and 75th, 10th and 90th, and 5th and 95th percentiles of the
distributions, respectively.

126

 

Relative error
N

 

   

 

 

DL1 ASP
Data generating model

Figure 3.2. Box plots representing relative error distributions for estimates of
exploitation rate in the last year of analysis across different data generating models. The
data generating and estimation models include double logistic fimctions with one (DL1),
two (DL2), and four (DL4) time-varying parameters, and time-varying age-specific
selectivity parameters (ASP). The bars represent median relative errors. The boxes,
whiskers, and circles represent 25th and 75th, 10th and 90th, and 5th and 95th percentiles
of the distributions, respectively.

127

Figure 3.3. The percentage of model runs when the model selection methods chose the
best or nearly best estimation model based on estimates of final population biomass. The
model selection methods include root mean square error (RMSE), deviance information
criterion (DIC), and degree of retrospectivity (DR). The best or nearly best estimation
model(s) is defined as the model(s) producing A) the lowest final population biomass
relative error, B) within 5% of the lowest final population biomass relative error, and C)
within 10% of the lowest final population biomass relative error.

128

Figure 3.4. The percentage of model runs when the model selection methods chose the
best or nearly best estimation model based on estimates of final exploitation rate. The
model selection methods include root mean square error (RMSE), deviance information
criterion (DIC), and degree of retrospectivity (DR). The best or nearly best estimation
model(s) is defined as the model(s) producing A) the lowest final exploitation rate
relative error, B) within 5% of the lowest final exploitation rate relative error, and C)
within 10% of the lowest final exploitation rate relative error.

130

 

 

 

 

 

 

131

 

3000 ' : RMSE O O
1500 - m DIC

10/ l

 

 

 

 

 

T
§ 3
e 6- \
a 4. \
m u .
2- s
o ---e~atae----e~§-m--

 

DL1 ASP
Data generating model

Figure 3.5. Box plots representing relative error distributions for estimates of population
biomass in the last year of analysis chosen by model selection methods across different
data generating models. The data generating models include double logistic functions
with one time-varying parameter (DL1) and time-varying age-specific selectivity
parameters (ASP). The model selection methods include root mean square error
(RMSE), deviance information criterion (DIC), and degree of retrospectivity (DR). The
bars represent median relative errors. The boxes, whiskers, and circles represent 25th and
75th, 10th and 90th, and 5th and 95th percentiles of the distributions, respectively.

132

6' IZIIRMSE

mDIC
—DR
L4.
8
5 .
.f‘>_’2-
g 0 O o o o
o:

 

 

 

DL1 ASP
Data generating model

Figure 3.6. Box plots representing relative error distributions for estimates of
exploitation rate in the last year of analysis chosen by model selection methods across
different data generating models. The data generating models include double logistic
functions with one time-varying parameter (DL1) and time-varying age-specific
selectivity parameters (ASP). The model selection methods include root mean square
error (RMSE), deviance information criterion (DIC), and degree of retrospectivity (DR).
The bars represent median relative errors. The boxes, whiskers, and circles represent
25th and 75th, 10th and 90th, and 5th and 95th percentiles of the distributions,
respectively.

133

    
 

V

1111111”

 

1111111111

 

956