**Beyond Theory to Application and Evaluation: Diffusion Approximations
for Population Viability Analysis**

Elizabeth E. Holmes

National Marine Fisheries Service

Northwest Fisheries Science Center

2725 Montlake Blvd. E.

Seattle, WA 98112

206(860-3369)

eli.holmes@noaa.gov

# Abstract

Census data on endangered species are often plagued by problems that
make quantitative population viability analysis (PVA) a challenge. This
paper addresses four such problems: sampling error, density-dependence,
non-stable age structure and population supplementation that masks the
true population status. Estimating trends and extinction risks using
such corrupted data presents serious parameter estimation difficulties.
Here I review diffusion approximation (DA) methods for estimating
population status and risks from time series data. A variety of
parameterization methods are available for DA models; some correct for
data corruption and others do not. I illustrate how stochastic Leslie
matrix models can be used to evaluate the performance of a proposed DA
model and to select among different DA parameterization methods for a
given application. Presenting the uncertainty in estimated risks is
critical, especially when the data are highly corrupted and estimated
parameters are more uncertain. Using a Bayesian approach, I demonstrate
how the level of data support for different risk levels can be
calculated using DA parameter likelihood functions.

Key words: population viability analysis, diffusion approximation,
process error, stochasticity, extinction, Dennis model, population
models, matrix models, salmon

Limited data are a common stumbling block for quantitative population
viability analysis (PVA). Sufficient data for developing detailed life
history models are often unavailable; indeed often the only available or
planned data are a simple time series of counts (Morris et al. 2002). In
the last decade, diffusion approximation (DA) methods have been
developed which use count data alone for the estimation of PVA risk
metrics, such as the probability of crossing extinction thresholds, mean
passage times, and average long-term rates of population growth or
decline (Lande and Orzack 1988, Dennis et al. 1991). These methods have
since been used to estimate extinction risks for numerous species of
conservation concern (Dennis et al. 1991, Nicholls et al. 1996, Gerber
et al. 1999, Morris et al. 1999, McClure et al. 2003). The appeal of DA
methods from an applied standpoint is their simplicity and their
reliance on simple census data alone, and they have become one of the
basic quantitative tools presented in recent books on PVA methods
(Morris and Doak 2002, Lande et al. 2003).

Diffusion approximation methods stem from theory concerning the behavior
of stochastic age-structured population models with no
density-dependence,

\[\]

where **A***~t~* is a stochastic population transition matrix, e.g. a
Leslie matrix, for time *t*. The **A***~t~*'s are assumed to be drawn
from some unspecified stationary statistical distribution, although they
need not be drawn randomly. For example, **A***~t+1~* might be drawn
conditioned on **A***~t~* if there is temporal correlation in
year-to-year environmental conditions (i.e. good years follow good
years). For such models, the asymptotic behavior of any weighted census
of the population, , is a stochastic exponential process (Tuljapurkar
and Orzack 1980, Tuljapurkar 1989):

\[\]

where e*~t~* is distributed normally with mean 0 and variance s^2^*t*
for *t* big. The statistical distribution of future population sizes,
here log *N~t~*/*N*~0~, is distributed normally with mean m*t* and
variance s^2^*t*, for *t* big. The weights *w~i~* can be equal to 1 so
that *N~t~* is simply the total number of individuals alive at time *t*,
but other weightings, such *w~i~* equals the fraction of age *i*
individuals that eventually become reproductive, may be more appropriate
when calculating extinction metrics. The parameter m in Eqn. determines
the rate at which the median log population size, log *N~t~*, increases
through time, while s^2^ determines the rate at which the distribution
spreads, or in other words, the variability of potential population
sizes at time *t*. Throughout this paper, s^2^ will be referred to as
process error. See Caswell (2001, section 14.3.3) for a review of the
theory of products of random matrices in the context of stochastic
age-structured models, and Morris and Doak (2002, chapter 3) for a
discussion of stochastic exponential models.

Diffusion approximation methods assume that Eqn. holds for all *t* \> 0
including small *t* and that the e*~t~* are independently and
identically distributed. This allows one to model the population as a
diffusion process (Lande and Orzack 1988):

\[\]

P(*y*\|*x*) means the probability of *y* given *x* The diffusion model
has the property that log *N~t~*/*N*~0~ is distributed normally with
mean m*t* and variance s^2^*t*, like the stochastic exponential process
it is used to approximate. See Dennis et al. (1991) for a fuller
discussion of the diffusion approximation.

This approximation opens a toolbox of parameterization methods for
linear models with normal error. It also provides analytical estimates
of passage probabilities, i.e. the probability of crossing a particular
threshold within a given time frame, for example extinction estimates
(to 1 individual) or quasi-extinction estimates to some critical
population size greater than 1. Strictly speaking, however, an
age-structured population process is not a diffusion process. For a time
series from an age-structured population, log *N~t~*~+t~/*N~t~* where t
= 1 is not normally distributed, only log *N~t~*~+t~/*N~t~* where t is
big, is normally distributed. Also in an age-structured population time
series, the temporal independence assumptions for t = 1 are violated
since the ratios, log *N~t~*~+1~/*N~t~* and log *N~t~*~+2~/*N~t~*~+1~,
will be correlated, even if the environment is uncorrelated.

Despite these assumption violations, the diffusion model approximates
many types of stochastic age-structured population processes, as seen
both from simulated and real data (Lande and Orzack 1988, Dennis et al.
1991, Holmes and Fagan 2002, Fagan et al. unpublished manuscript). There
are, however, well-known cases where the diffusion approximation
performs poorly, such as when year-to-year variability is high and
populations are small (Ludwig 1996b, 1999), when demographic
stochasticity introduces strong non-linearities (Wilcox and Possingham
2002), and when density-dependence is extreme (Sabo et al. 2003). Thus,
it is critical to carefully evaluate the appropriateness of the
diffusion model given the life history of the species at hand and given
the types of extinction or quasi-extinction metrics one would like to
estimate. In addition, careful selection of parameterization methods for
estimating the DA parameters is essential since many real population
time series will contain extraneous variability whether due to
age-structure cycles, density-dependent feedback, or sampling error, and
this will confound proper parameterization. In this paper, I discuss and
illustrate the application of DA methods with the biologist faced with
assessing risks for a real population in mind. I discuss the available
parameterization methods and present a quantitative approach for
assessing the performance the diffusion approximation and its
parameterization.

Once a parameterization method is chosen, evaluation of which risk
metrics to use is needed. Some metrics will be robust and others may be
useful only over certain time frames. In particular, estimates of the
probability of crossing population thresholds (e.g. extinction metrics)
can be extremely variable for certain parameter ranges (Fieberg and
Ellner 2000, Ellner et al. 2002). This is not a problem specific to the
DA approach; it is a feature of stochastic population processes in
certain parameter ranges. Nonetheless, conservation biologists are often
asked to make an assessment of extinction or quasi-extinction risk;
indeed, such an assessment may be mandated by law when working with
species being considered under the U.S. Endangered Species Act. One of
the positive features of DA methods is that uncertainty is well
characterized and can be estimated from the data. Thus, while we may be
uncertain about the point estimates of the risk metrics, we can be
precise about the level of uncertainty. In this paper, I discuss two
different methods for representing uncertainty: confidence intervals and
posterior probability distributions.

These methods for evaluating and presenting DA risk metrics are
illustrated using three simulated populations of different salmonid
species. These concrete examples are intended to illustrate analyses and
models that can be used to assess the performance of a proposed DA
approach given the life histories and data corruption one expects in a
particular PVA application. The paper is divided into four steps,
reflecting the steps of evaluating a proposed DA application: 1)
evaluating whether a diffusion approximation is reasonable for the
population at hand, 2) evaluating and selecting a parameterization
method, 3) evaluating the sensitivity of proposed risk metrics to
parameterization errors, and 4) estimating parameters and risk metrics
given a census time series from the concerned population and presenting
the uncertainty in these estimates.

**Evaluating the diffusion approximation for the population**

*Simulated data*

The first step in evaluating a proposed DA application is to develop a
model that captures the basic life history and that has realistic
environmental variability, including the temporal correlation in the
environment. Although a precise model for the population will likely not
be available, given that diffusion approximations are typically used
where data are limited, often enough natural history information will be
available to roughly parameterize a model for the purpose of evaluating
the performance of proposed methods. Such an evaluation is illustrated
here using stochastic matrix models for Snake River spring/summer
chinook, Snake River fall chinook, and upper Columbia River steelhead,
which are currently listed under the U.S. Endangered Species Act.
Although this paper uses matrix models, the simulated data could be
produced by any variety of models or by using a bootstrapping method to
simulate from actual time series data, e.g. Hinrichsen (2002).

The basic model used here tracks individuals from age 1 to spawners:

\[\]

where **A***~t~* is the Leslie matrix model describing the transition
probabilites from age *i* to *i*+1. The stochastic matrices used for the
three salmonid populations are listed below with the parameter values
defined and given in Table 1.

[Snake River spring/summer chinook]{.underline} The following matrix
model was adapted from Kareiva et al. (2000) and modified to separate
out the spawner class. The parameters estimates were based on estimates
for Marsh Creek 1980-1998 from M. McClure (National Marine Fisheries
Service):

[Snake River fall chinook]{.underline} This model is also adapted from
Kareiva et al. (2000) with parameter estimates based on 1980-1998 from
M. McClure (National Marine Fisheries Service):

[Upper Columbia steelhead]{.underline} The following Leslie matrix was
estimated for Methow River steelhead for pre-1998 conditions (T. Cooney,
National Marine Fisheries Service).

Stochasticity was added to survivorships and fecundity by multiplying
the survivorship or fecundity terms by a random variable e*~i~* with the
maximum possible survivorship capped at 1.0. The e~1~ term represents
variability in spawner to age 1 ratios. For all matrices, e~1~ variables
were drawn from a lognormal distribution specified by exp(normal(mean=0,
variance=0.13)). The variance was estimated from spawner to smolt (age
1.5) data for spring/summer chinook in the Snake River (P. Levin,
National Marine Fisheries Service). The e~2\ ~variables associated with
age 1 to age 2 survivorship for spring/summer chinook and steelhead were
drawn from a lognormal distribution specified by exp(normal(mean = 0,
variance = 0.08)). The variance was estimated from parr to smolt (age
0.5 to age 1.5) survivorship data for spring/summer chinook in the Snake
River (Achord et al. 2003). The e~1~ and e~2\ ~variables were set to
co-vary slightly (correlation = 0.2) since both are associated with
basin conditions, although one is associated more with spawning ground
conditions and the other is associated more with downstream migration.
The e~o\ ~variables associated with ocean survivorship were drawn from a
lognormal distribution specified by exp(normal(mean = 0, variance =
0.02)). The variance for e~o~ was chosen such that the resulting s^2^
estimates at the population level reflected approximately that estimated
from actual time series. The e~o~ variables associated with different
ages were set to co-vary strongly (correlation coefficient = 0.8) since
they are all associated with ocean conditions.

In addition, studies of ocean condition indices and studies of cohort
(i.e. all young of year *t*) survivorship data for Columbia River
salmonids suggest that there is also strong year-to-year correlation in
ocean survivorship, i.e. good ocean survivorship one year tends to be
followed by good ocean survivorship the next year. To model this, the
mean ocean survivorship across ages was chosen such that it co-varied
strongly (correlation coefficient = 0.7) with the mean ocean
survivorship in the previous year. The level of correlation was chosen
from the levels seen in time series of cohort survivorship data.
Temporally correlated juvenile survivorship (in-stream survivorship) was
not included since evidence of this was not see in the study by Achord
et al. (2003).

Although current densities of salmon in these three populations are well
below historic levels, density-dependence in parr to smolt survivorship
has been shown from studies on Snake River spring/summer chinook (Achord
et al. 2003). Achord et al. postulate that juveniles continue to show
density-dependence at low density because of nutrient limitation from
the associated low densities of spawner carcasses. Density-dependence
was incorporated into the egg to age-1 survivorship term (*s*~1~) by
fitting a linear model to the juvenile survivorship data in Achord et
al. (their figure 3), giving

\[\]

where *s*~1~ is the mean egg to age-1 survivorship in Table 1 and
*N~s,t~* is the number of spawners in year *t*. The minimum was set at
0.39*s*~1~ and the maximum at 1.17*s*~1~. The effect of Eqn. was an
increase in of up to 17% over *s*~1~ when spawner density was lower than
the previous 10-year mean spawner density and up to 61% lower than
*s*~1~ survivorship when spawner density was higher than the 10-year
average.

Using the stochastic matrices and Eqn. , one thousand population time
series were generated for each species. The time series were started
from a stable age structure drawn at random from the stochastic
distribution of stable age structures. The distribution of log
*N~t~*/*N*~0~ from the one thousand time series describes the
statistical distribution of population size at time *t*.

### Performance of the diffusion approximation for simulated salmonid time series

The goal was to determine whether a diffusion model exists that
describes the future distribution and quasi-extinction risks for the
number of current and future spawners within a salmon population. A
weighted total population at time *t* was defined as where *w~i~* equals
the average fraction of age *i* fish that will eventually become
spawners. This was used instead of simply since the age structure is
heavily dominated by age-1 individuals due to low age-1 to age-2
survivorship and using would mean effectively tracking the dynamics of
age-1 individuals alone rather than the individuals with reproductive
potential across all ages.

First the predicted versus observed distribution of log *N~t~*/*N*~0~
was tested. The predicted distribution is normal with a mean and
variance that increase linearly with time, such that (1/*t*) times the
mean and variance of log *N~t~*/*N*~0~ equal constants. Figure (a-c)
shows the mean and variance of log *N~t~*/*N*~0~, divided by *t* for the
simulated salmon populations. This figure illustrates that the mean and
variance increased linearly with *t* as predicted by theory. The mean
and variance of log *N~t~*/*N*~0~ divided by *t* reached a constant
after approximately 5 years (Fig. a-c). The distribution of log
*N~t~*/*N*~0~ approached normality on the same time frame, as determined
by the difference between the cumulative distribution of log
*N~t~*/*N*~0~ versus the cumulative distribution of a normal with mean
m*t* and variance s^2^*t* (not shown). This evaluation indicates that
for time periods past 5 years in the future, there exists a diffusion
model that will properly describe the statistical distribution of
population sizes at time *t* \> 5. The parameters, m and s^2^, for this
diffusion model are given respectively by the two constants: (1/*t*)
times the mean and variance of log *N~t~*/*N*~0~ as *t* increases. The
next step is to ask whether and over what time frames, a diffusion model
gives reasonable approximations for the extinction or quasi-extinction
metric of interest.

Even though the distribution of log *N~t~*/*N*~0~ for *t* not too small
may be properly characterized by a diffusion model, the probability of
crossing thresholds may not be since the variance of log *N~t~*/*N*~0~
is higher than s^2^*t* for *t* \< 5 and since there is temporal
correlation in the population process due to age structure in the
population model. Thus the extinction or quasi-extinction metric itself
needs to be tested. For this paper, the extinction metric of interest
was the probability of 90% decline within different time frames.
Although estimating extinction to 1 individual is a popular risk metric,
and unfortunately sometimes mandated, the reader is cautioned about
using the diffusion approximation is to estimate extinction to very low
numbers since factors that drive dynamics at very low population sizes
(such as demographic stochasticity) and the catastrophes often
associated with ultimate extinction will likely be poorly represented in
a time series of a relatively larger population declining to low numbers
(e.g. Fagan et al. unpublished manuscript).

The diffusion approximation was tested by comparing the actual
probabilities from the one thousand simulations to those predicted by
the diffusion model with parameters m and s^2^ from Fig. (a-c). For the
diffusion model, the probability of decline to a threshold population
size *N*~crit~ within *t* years is (Lande and Orzack 1988, Dennis et al.
1991):

. \[\]

For 90% decline, *N*~0~/*N*~crit~ = *N*~0~/(0.1 *N*~0~) = 10. F(*y*) is
the value of the cumulative distribution function at *y* of a standard
normal distribution with mean 0 and variance 1. Equation gives the
probability of 90% decline at any time within 0 to *t*, however in the
simulations, population size is only observed at discrete yearly
intervals. In this application, the difference between Eqn. and the
technically correct comparison, the probability of observing
quasi-extinction at a discrete set of yearly intervals, is minor, but
this is not always the case. If needed, the probability of observing
quasi-extinction at discrete intervals can be calculated numerically by
simulating the diffusion process.

Figure (d-f) shows the actual (gray lines) versus diffusion
approximation estimates (black lines) of the probability of 90% decline.
The plots illustrate that the diffusion model correctly described the
probability of 90% decline within different time frames for the
simulated time series of current and future spawners. This was true
despite the fact that these time series were not strictly speaking a
diffusion process. These results are for a particular passage metric of
interest, i.e. the probability of 90% decline. The results should not be
over-generalized to say that a simple diffusion model would correctly
describe all extinction and quasi-extinction metrics for salmonids. It
is important to test particular metrics of interest such as done here
for the probability of 90% decline. Also it is important to consider
what segment of the population to track. 'Future and current spawners'
integrates over multiple age-classes without over-representing any one
age class and smoothes out the year-to-year boom-bust cycles that salmon
are prone to. This effectively limits the non-process error in the time
series and enables Eqn. , which is based on a diffusion model with no
non-process error, to work well.

**Parameter estimation for diffusion approximations**

This first evaluation indicated that there exists a diffusion model that
correctly approximates the behavior of these simulated salmon
populations for the purposes of estimating the distribution of future
reproductive population sizes and for estimating the probability of
hitting a critical level of 10% of current size. The next question is
how well the parameters of this diffusion model can be estimated given
realistic data constraints. Here I review currently available methods
for parameterizing a diffusion model: maximum likelihood methods
assuming low data corruption (Dennis et al. 1991), the same methods but
using a running sum transformation of the data (Holmes 2001), maximum
likelihood methods using a Kalman filter (Harvey 1989, Lindley 2003),
slope methods (Holmes 2001, Holmes and Fagan 2002), and asymptotically
unbiased estimators for random matrix products (Heyde and Cohen 1985).
Code for these methods is given in the electronic Supplement. The
different methods have their pluses and minuses depending on the quality
and length of data available. Below these methods are described and in
the next section, I evaluate the performance of these different methods
given the data constraints faced by a PVA of Columbia River endangered
and threatened salmonid populations (McClure et al. 2003).

All these methods start with a time series of censuses, *O*~1~ *O*~2~
*O*~3~ ... *O~k~*, of the population. For the purposes of this paper, I
assume *O~t~* represents yearly censuses. The census need not enumerate
the entire population; index counts, such as a segment of the population
or a specific age or stage class, can also be used.

### Maximum likelihood methods

I use the following state-space model to model the relationship between
the diffusion model we wish to estimate from the data, the actual
population counts, and the observed counts:

*X~t~~\ ~*= *X~t-1~* exp(*m + e~p~*) where *e~p\ ~*\~
normal(mean=0,variance=*s*^2^)

*N~t~~\ ~*= *X~t~* exp(*e~n~*) where *e~n\ ~*\~ g(mean=0,variance=) \[\]

*O~t~* = *N~t~* exp(*e~se~*)^\ ^ where *e~se\ ~*\~ f(mean=b,variance=)

where '\~' means 'is distributed as' and g(^.^) and f(^.^) are some
unspecified statistical distributions. *X~t~* represents the diffusion
process, which correctly approximates the age-structured population
process for *t* not too small. *N~t~* are the true population counts,
which include process error (the *X~t~* part) and non-process error, ,
which is variability due for example to density-dependent feedback or
age-structure cycles. Non-process error is apparent in the time series,
but does not multiply through time as e~p~ does in Eqn. . Observations
of *N~t~* are made and these observations are corrupted by some sampling
error, e~se~. Equation reduces to:

*X~t~~\ ~*= *X~t-1~* exp(*m + e~p~*)

*O~t~* = *X~t~* exp(*e~n~* + *e~se~*)^\ ^ where *e~n~* + *e~se\ ~*\~
h(mean=b,variance=) \[\]

where h(^.^) is an unspecified statistical distribution. In this paper,
I refer to this state-space model as a corrupted diffusion process. The
total non-process error, , equals non-process error in the
age-structured population process, , plus non-process error from
sampling error in the censuses of the population, .

If we assume that there is no non-process error in our observations, = 0
and *O~t~* = *X~t~*. Then the maximum likelihood estimates of m* *and*
*s*^2^* are (Dennis et al. 1991):

\[\]

Equation assumes that there are no missing data; cf. Dennis et al.
(1991) for methods when there are missing data. Although the assumption
of no non-process error in the data will never strictly hold for real
data, analysis of a variety of non-salmon time series suggests that it
can be a reasonable approximation (Fagan et al. manuscript).

The estimates and are the maximum likelihood estimates for the diffusion
approximation (Eqn. ) not of the age-structured process itself. This
creates problems if the life history of the species of concern is such
that there is high non-process error in the time series due to
age-structure cycling. If there is also sampling error in the data, this
adds additional non-process error variability. Thus Eqn. will produce
biased estimates of s^2^ because it attributes the non-process error in
the data to process error, and the expected value of , E(), is *s*^2\ ^+
2. In this case, it is more appropriate to use the full state-space
model (Eqn. ) with not assumed to be zero. Kalman filters are widely
used in time series analysis for maximum likelihood estimation for
state-space models of the form of Eqn. with e~np~ assumed to be normally
distributed (Harvey 1989, section 3.4). Lindley (2003) gives a specific
case of such an algorithm for population processes of the form of Eqn. .

This approach for maximum likelihood estimation for Eqn. can be
summarized as follows. Assume h(**^.^**) in Eqn. is normal with b = 0
(if b is non-zero it can be factored out) and define , then the
distribution of *y~t~* given all the data prior to *y~t~* is:

\[\]

where E(**^.^**) is 'expected value' and *F~t~* denotes the variance of
. The likelihood of a particular set of m, and s^2^ given the data, , is

and the log likelihood is

+constant \[\]

where . The Kalman filter calculates *v~t~* and *F~t~* given a
particular set of m, and s^2^. The Kalman filter for this application is
given in Appendix A or see Harvey (1989, chapter 3). The maximum
likelihood estimates are the set of m, and s^2^ that maximize Eqn. and
are denoted here as and .

A variant of the maximum likelihood methods uses a running sum
transformation of the data in place of *O~t~* in Eqn. . A running sum is
defined in this case as

\[\]

where the *a~i~* are weightings. The running sum estimate for m is

\[\]

Often the purpose of the running sum is to transform the data into a
population-level count (see examples in Dennis et al. 1991 and Holmes
2001). But the transformation has an independent benefit when data are
corrupted by extraneous variability. In this case, can have lower
variability and bias than (Holmes 2001, Holmes and Fagan 2002). The
process error, s^2^, cannot be estimated using *R~t~* in Eqn. . The
running sum transformation filters out not only the extraneous
variability but also the process error that we are trying to estimate.

### Slope method

The slope method was developed to deal with situations where *O~t~* is
corrupted with high non-process variability, , that is up to one to two
orders of magnitude greater than s^2^ and has an unknown statistical
distribution (Holmes 2001, Holmes and Fagan 2002). It also uses the
corrupted diffusion model, but estimates the parameters by examining how
the distribution of changes with t. From Eqn. :

\[\]

where the t are chosen such that consecutive *O~t+~~t~~\ ~/O~t~* ratios
do not overlap. The slopes of the linear regressions of the mean and
variance of log(*O~t+~~t~~\ ~/O~t~)* versus *t* give estimates of m and
s^2^, respectively. An unbiased estimator for s^2^ can devised using
Eqn. with a long enough time series, however when the available time
series is short (such as \< 20 years) such an estimator is prone to
negative variance estimates. To stabilize the estimates, regression can
be done on the running sum transformed data, *R~t+~~t~~\ ~/R~t~* instead
(Holmes and Fagan 2002):

\[\]

This slope estimate for s^2^ is biased; from numerical simulations
E(0.5*s*^2^+0.15. However, if non-process error is less than an order of
magnitude of the process error, the result is mean estimates that are
within 50-200% of the correct value.

*Asymptotically unbiased estimator for s*^2^

The previous methods are all based on the diffusion approximation for
log *N~t~*~+t~/*N~t~* which has variance s^2^t for very small t. The
s^2^ estimates from these methods will tend to be biased to some degree
even if the time series is infinitely long. This problem is worst for ,
which relies on the variance at t = 1 for estimation of s^2^. Heyde and
Cohen (1985) provide an estimator for s^2^ that is not based on the
diffusion approximation but rather on products of random matrices:

\[\]

This estimator is asymptotically unbiased, i.e. unbiased as the length
of the census time series goes to infinity, although not necessarily
unbiased for short time series as shall be seen. The asymptotically
unbiased estimator for m is the same as .

*Correcting for inputs*

All these methods assume that there are no external inputs into (or
subtractions from) the population. If inputs are occurring, they will
mask the true m and thus must be corrected for to prevent bias in the
estimation of m. In Appendix B, a general method to correct for the
input problem is derived. The formulation of the correction depends on
the life history of the species at hand and the ages at which
individuals are externally input into the population. Thus a specific
correction will have to be derived following the methods in Appendix B
for the species of concern. Here the calculation is illustrated for
salmon populations experiencing regular hatchery fish introductions.

Hatchery-reared juvenile fish are regularly released into Pacific salmon
stocks as part of remediation for impacts or to provide fishing
opportunities. If these hatchery-born fish return to reproduce in the
wild, their offspring are indistinguishable from the offspring of
wild-born fish, for typical census purposes at least. The goal of the
input correction is to estimate what m would be if no supplementation
were occurring. As described in Appendix B, correcting for inputs
requires determining the relationship between the true mean number of
offspring per wild-born spawner at year *t*, denoted *R~0,t\ ~*, to the
apparent mean number, denoted , based on the population growth rate with
inputs occurring. The species in this study reproduce once and die,
which makes a relatively straightforward calculation (cf. Appendix B):

\[\]

where *f~w,t~* is an estimate of the fraction of returning spawners at
year *t* that were wild-born. For notational simplicity, it is assumed
here that hatchery fish reproduce at the same rate as wild-born fish. If
hatchery fish reproduce at lower effectiveness, *S~h,t~* can be
re-expressed as wild-equivalents by multiplying *S~h,t~* by the
reproductive effectiveness of hatchery fish relative to wild-born fish.

To solve for the no supplementation population growth rate at year *t*,
denoted l*~t~*, in terms of the observed with supplementation growth
rate, denoted , the relationship was used. Here *T* is mean generation
time (Caswell 2001, pgs. 126-130). Combining this relationship with Eqn.
gives:

\[\]

As discussed in Appendix B, could be estimated from the ratio of wild
spawners at year *t+1* versus year *t*, e.g. (), or by the ratio of
age-1 fish, e.g. (), where *F~t~* is the fecundity of spawners. However,
for Pacific salmon populations, *f~w,t~* is often approximate, and
yearly estimates of *F~t~* are unusual. Instead, I will use here the
total spawner ratios,, to estimate of . This is an approximation since
it tacitly assumes that either *f~w,t~* or *F~t~* remain relatively
constant from year to year, but reflects an approximation that is
required given real data constraints. The input-corrected m estimate is
then:

\[\]

*Overview of parameterization performance*

Table 2 gives an overview of the four different parameterization methods
in terms of whether they correct for extraneous variability in the time
series and some rough guidelines regarding their data needs. The table
emphasizes that there is a trade-off between data needs and the ability
to correct for error in the data. When the data contain high extraneous
variability, the reduction in parameter bias warrants the increase in
the variability of parameter estimates. The flip side of this is that
using a method that corrects for error on time series with low errors
will lead to an unnecessary loss of precision.

**Performance of parameterization methods given data constraints**

Although Table 2 can help select likely candidate methods for a
particular application, a more quantitative evaluation of the
performance of particular parameterization methods is needed since
actual census data have a level extraneous error somewhere between none
and severe and thus the trade-off between precision and bias is unclear.
Such an evaluation is illustrated here using the stochastic matrix
models for the three salmonid populations. This evaluation incorporates
the following data constraints faced by a PVA of these populations
(McClure et al. 2003): 1) counts of only the spawning segment of the
populations, 2) time series limited to 20 years, 3) severe age-structure
perturbations in the beginning of some time series due to reproductive
collapses during dam construction (Williams et al. 2001), and 4) high
sampling error.

To examine the robustness of the parameter estimation methods to these
constraints, 20-year time series of spawner only counts were generated
using the matrix models. The time series were generated using either an
initial stable age structure (drawn randomly from the stable set) or an
initial age structure with no age-1 individuals. The simulated spawner
time series, *S~t~*, were then corrupted with lognormal sampling error:

*O~t~* = *S~t~* exp(e~se~) with e~se~ \~ normal(mean=a, variance=). \[\]

The observed spawner count at year *t* is *O~t~*. The DA parameters, m
and s*^2^*, were estimated from each 20-year simulated time series using
the Dennis, Kalman, running sum, Heyde-Cohen or slope methods. For the
running sum and slope methods, spawner counts were transformed by adding
four consecutive counts: *R~t~* = . Parameter estimates were compared
with the correct DA parameters determined from Fig. .

Realistic levels of sampling error to add were estimated from studies on
sampling error in spawner and redd surveys (Jones et al. 1998, Dunham
and Rieman 2001) and from an examination of the average non-process
error in Columbia River redd count data. Redds are the egg nests made by
spawning salmon, and redd counts were the most common data type for the
Columbia River PVA in McClure et al. (2003). Using a lognormal model for
observation errors, Jones et al. (1998) found levels of total within
plus between observer sampling error variability of = 0.09 to 0.78.
These levels of lognormal sampling variability are consistent with those
found in another study of variability between different observers'
counts of redds (Dunham and Rieman 2001). To find the average total
non-process error (which includes sampling variability plus other
non-process error) in actual time series, I examined a collection of 44
twenty-year redd-count time series from different spawning areas in the
Snake River basin, Idaho (data in supplements of McClure et al. 2003). I
estimated an average total non-process error of = 0.56 in this
redd-count data. The empirical studies on sampling errors combined with
the average total non-process error within redd-count time series
suggest that year-to-year sampling error variability is likely to be in
the range of = 0.1 to 0.75, with 0.1 on the low end given minimum
observer variability alone and 0.75 on the high end given that this is
more than observed levels of total non-process error. These levels of
sampling error are approximately 10 to 100 times larger than the median
process error estimates from actual salmonid time series in the Columbia
River basin.

To model these sampling error levels, lognormal sampling error was added
to the spawner counts using Eqn. with low ( = 0.1), medium ( = 0.5) or
high ( = 0.75) sampling error. The bias in error, a, drops out since the
parameter estimates always use the difference between log counts (e.g.
log O*~t+1~* -- log O*~t~*). Jones et al. (1998) found that observers
have a non-linear tendency to underestimate as the number of objects to
be counted increases. However, they found this for counts of tens of
thousands of salmon; Dunham and Rieman (2001) did not find this pattern
for redd counts in the hundreds which are more typical of the counts for
the endangered and threatened salmon populations simulated in this
study. Thus, I assumed that the sampling error was independent of total
spawner numbers.

Figure shows the performance of the m estimators with either no sampling
error and a stable starting age structure, with only an age-structure
perturbation and no sampling error, or a stable initial age structure
and low, medium or high sampling error. In the box plots, the middle
line is the median and the box encloses 75% of the estimates from the
1000 simulations. All methods gave unbiased estimates of m for the
simulations started with a stable age structure. The Dennis and slope
estimates of m were most variable, and the runsum and Kalman estimates
were generally least variable. The age-structure perturbation led to an
underestimation of m that was most apparent for the Dennis estimate,
while the other three methods showed only slight underestimation. For
the runsum and Kalman estimates, 75% of the estimates were within +/-
0.05 of the correct value for low and medium sampling errors. For high
sampling error, the range of the estimates increased to 75% of estimates
between +/- 0.08 of the correct value. These errors can be translated
into errors in percent yearly growth rate by multiplying by 100.

Figure shows the performance of the s*^2^* estimators. The Dennis method
overestimated s^2^ for all cases, with and without sampling error. For
example, with no sampling error and a stable starting age structure, the
median estimates were 0.35, 0.12, and 0.08 versus the correct values of
0.033, 0.042, and 0.040 for spring/summer chinook, fall chinook, and
steelhead, respectively. The bias was due to the tendency of the salmon
simulations, like real salmon populations, to exhibit boom and bust
cycles. The result was high non-process error in the time series even
without added sampling error. With a perturbed age structure, the more
highly overestimated s^2^, and with sampling error added in addition,
the overestimation was severe. The median were 2.02, 1.74 and 1.71 for
the three species with the high added sampling error. The Heyde-Cohen
estimates were similarly highly biased with high non-process error.

In contrast to the overestimation by these estimators, the Kalman
estimator was prone to underestimating the correct s^2^; median Kalman
estimates were close to 10^-16^ for spring/summer chinook across all
corruption levels and for fall chinook and steelhead for the higher
sampling error levels. Even when the median estimates were close to the
true values (i.e. for fall chinook and steelhead with low sampling
error), s^2^ estimates close to zero were still common and comprised
20-30% of the estimates. Overall for the 20-year time series especially
the time series with high sampling error levels, the data were usually
most likely under a scenario with all the variability attributed to
non-process error and thus s^2\ ^= 0. In contrast, the mean estimate of
s^2^ (rather than the median) was uniformly close to the true s^2^
(similar to the results by Lindley (2003)). This difference between the
mean and median estimates occurred because the distribution of estimates
was highly skewed, with most estimates essentially zero except for a few
very high estimates, which brought up the mean. This difference between
the median Kalman estimates and the correct DA values diminished if a
longer time series was used, and largely disappeared for 100-year time
series.

The slope method gave estimates that were overall closest to the true
s^2^ for the 20-year time series, although with no sampling error added
to the time series, they tended to underestimate the true value as
predicted. With a stable initial age structure and no sampling error,
the median estimates were 0.025, 0.022, 0.018 relative to the correct
values of 0.033, 0.042, and 0.040 for the three species, respectively.
The estimates were not appreciatively changed by the age-structure
perturbation nor low sampling error. As sampling error increased to the
medium and high levels, however, the slope estimator began to increase
and for the highest sampling error the median estimates were 0.12, 0.12,
0.12, for the three species. These estimates were still substantially
closer to the true values than the Dennis and Heyde-Cohen estimators at
high sampling error levels.

*Testing the input corrector*

To study the robustness of the input correction (Eqn. ), simulations
were run as above but with a random number of age-1 hatchery fish added
each year. Between 1 and 100,000 fish were added each year; the number
drawn from a uniform random distribution. All simulations were started
with a stable age structure. Each simulation produced wild and hatchery
spawner counts of which observations, *O~w,t~ *and *O~h,t~*, were made
with medium sampling error. The parameter m was estimated using with
either no input correction (Eqn. with *O~t~* = *O~w,t~* + *O~h,t~*) or
with input correction (Eqn. ). The wild fraction was estimated from the
*O~w,t~* and *O~h,t~* counts:

\[\]

The slope method (Eqn. with ) was used to estimate s*^2^*.

The top panel of Fig. contrasts the input-corrected versus the
uncorrected m estimates. Without correction for inputs, the median
estimates of m were 0.00, 0.02 and 0.00 indicating stable or increasing
populations with the inputs relative to the correct values of --0.093,
-0.015 and --0.18, indicative of declining wild dynamics for the three
species, respectively. With correction, the median estimates were
--0.096, -0.023, and --0.18. Thus, the input correction successfully
extracted the correct value from population trajectories that otherwise
appeared to come from a process with much higher m due to hatchery
inputs. The lower panel illustrates that the estimate with no explicit
hatchery correction gave a relatively unbiased s^2^ estimate despite the
variable hatchery inputs.

# Sensitivity of risk metrics to parameter misestimation

The errors in parameter estimation seen in Figs. and translate into
errors in the estimated risk metrics. Two risk metrics were evaluated
for their sensitivity to parameterization errors: the long-term average
population growth rate, l = exp(m) and the probability of a 90% decline
within a given time frame (Eqn. ). Figure compares the estimates of l.
The figure shows the probability distribution, i.e. P(a \< l \< a +
da)/da. Since the variability of m estimates was not dramatically
different among the estimators, the variability in l using only is
shown. The figure illustrates that 50% of the estimates were within
+/-0.04 of the correct value even for high sampling error; although,
much lower and higher estimates were not uncommon. With low sampling
error 10% of estimates were more than +/- 0.07 greater than the correct
values for the three species. For high sampling error, this increased to
+/- 0.10. This level of variability indicates that unless the population
is declining very rapidly, a 20-year time series is unlikely to be
sufficient to reject a null hypothesis that the population is stable or
increasing (l \> 1) at the P=0.05 level using a standard t-test.

Figure compares the estimates of the probability of 90% decline from
simulated time series with a stable initial age-structure and no
sampling error added. Three methods for estimating the parameters were
contrasted. The first used the Dennis estimators, which assumed zero
non-process error (as in Dennis et al. 1991). The second method used the
Kalman estimates with normal non-process error. The third used the slope
method for s^2^ and the running sum method for m (as in Holmes and Fagan
2002). Due to high overestimation of s*^2^* by , the median probability
estimates using the Dennis estimators were highly biased relative to the
observed probabilities from the simulated spring/summer chinook and fall
chinook time series (Fig. , left panels). The probability of 90% decline
for steelhead was dominated by the very low m, and was less affected by
overestimation of s^2^. The majority of Kalman estimates for s^2^ were
severe underestimates for spring/summer chinook (Fig. ), and the median
probability estimates were correspondingly biased for this species (Fig.
, middle panels). For steelhead and fall chinook, the median Kalman
estimates for s^2^ were close to the true values when no sampling error
was added (Fig. 3), and in this case the median probability estimates
were correspondingly close to the true probabilities. The slope
estimates of s^2^ were generally close to the true values for all
species, and the median probability estimates using matched the observed
probabilities across species.

Estimation of the probability of 90% decline using the parameter
estimates, and , was also tested for sensitivity to sampling error
levels. Estimation was generally robust to low to medium sampling error
(except for fall chinook), but a mismatch between median predictions and
the true probabilities became apparent with high sampling error (Fig. )
when the difference between and s^2^ increased (Fig. ). This was
especially true for fall chinook with the predicted probabilities
bearing little relationship to the true low probabilities. Variability
in the estimated probabilities was highest for fall chinook and lowest
for steelhead. These differences in the variability of estimates between
species was likely due to differences between their rates of decline; l
was closest to 1.0 for fall chinook (with l = 0.98) and farthest for
steelhead (with l = 0.84). Figure illustrates that in some cases, e.g.
fall chinook with l close to 1.0, estimated probabilities of severe
declines are likely to be so variable as to provide little useful risk
information. In other cases, however, e.g. spring/summer chinook and
steelhead with ls farther from 1.0, the median estimated probabilities
correctly track the true probabilities and although broad, the 75% range
of estimates tracks the true probabilities. Nonetheless the probability
estimates were variable and accurately representing the variability is
critical to the use of probability metrics in applications, as discussed
in the next section.

# Application to actual data

The previous sections have focused on using simulated data to evaluate
the appropriateness of the diffusion approximation for a given
application and to select parameterization methods. Once this is
completed, the next step is to use the diffusion approximation to
estimate risk metrics using the actual time series data from the
population of concern. This final step is illustrated using a 38-year
times series from spring chinook in the Upper Columbia River basin (data
from T. Cooney, National Marine Fisheries Service). To analyze this
data, two approaches were used. The first focuses on point estimates and
confidence intervals. The second uses a Bayesian approach and estimated
the probability of the parameters and risk metrics given the data. This
approach gives a measure of the data support for different potential
true risk levels. Preliminary diagnostic tests for the data, namely
testing log *N~t+1~*/*N~t~* for normality, outliers, and serial
correlation, are not discussed here, but are reviewed in Dennis et al.
(1991).

*Point estimates and confidence intervals*

The approach here consists of four steps: 1) select estimators and make
the needed assumptions about their expected values, 2) calculate the
point estimates, 3) specify their expected distributions, and 4) use
parametric bootstrapping to estimate confidence intervals.

1\) Estimators and their expected values: The previous analyses suggest
that given the expected levels of process and non-process error in
salmonid data:

\[\]

2\) Point estimates: From the time series, the following estimated were
made: , and . Using the assumptions in Eqn. , the point estimates of the
risk metrics can be calculated. The point estimate of l for this
population is , indicating a 7% per year decline. The point estimate for
a 90% decline within 25 years uses Eqn. with and : P(*N*~25~ \< 0.1 x
*N*~0~) = 0.44.

3\) Their estimated distributions: The statistical distributions of the
parameter estimates are themselves estimated by using the known
distribution of the parameter estimates (Holmes and Fagan 2002) and in
place of the expected or true value of the parameters, which appear in
these distributions, using the point estimates for the parameters (from
step 2):

is the estimated distribution

and using as an estimate of

is the estimated distribution

given and using \[\]

is the estimated distribution

and using the estimated s^2^ and

The parameter *n* is the length of the time series, *L* is the number of
counts summed together to form a running sum, and *df* 0.333 + 0.212*n*
-- 0.387*L*. The distribution gamma(b,w) denotes a gamma distribution
with scale b and shape w, a c^2^(a) denotes a chi-squared distribution
with a degrees of freedom and a normal(*a*,*b*) denotes a normal
distribution with mean *a* and variance *b*.

4\) Confidence intervals via parametric bootstrapping: The confidence
intervals were estimated by generating thousands of parameter estimates
from their estimated distributions (Eqn. ). First and were drawn using
their estimated distributions, and then a was generated using the and
draws. For each parameter set, l and P(*N*~25~ \< 0.1 x *N*~0~) were
calculated. This was repeated 10,000 times. The 95% confidence intervals
contain 95% of the estimates: l 95% confidence intervals were (0.85,
1.04) and P(*N*~25~ \< 0.1 x *N*~0~) 95% confidence intervals were
(0.04, 0.94).

### Posterior probability distributions

The point estimates do not take into account the uncertainty in the
parameter estimates and the confidence intervals, which do, are
difficult to use in a decision-making framework. For example, the 95%
confidence intervals on l include 1.0. However, this population has
declined 93% over the last 38 years, and the data are most consistent
with a l less than 1.0, even though it is possible for such a decline to
have occurred by chance in a population with long-term non-declining
dynamics (i.e. true l \> 1). Statistical decision theory (e.g. Berger
1985) presents a framework for incorporating uncertainty in the true
risk levels and the costs of different management decisions given
different true risk levels. Wade (2000) and Dorazio and Johnson (2003)
provide recent discussions of this Bayesian decision framework in
conservation biology and resource management contexts. This framework
relies on calculating the probability of different risk levels (i.e. the
posterior probability distributions) given the uncertainty in the
underlying true parameters, in this case in m, s^2^, and .

Calculation of the posterior probability distribution for a given risk
metric,Y, involves calculating the probability that the risk metric is
some particular value, f, given the data. This is calculated by
integrating "the probability of the data given a particular set of
parameter values, q, times the probability of that q set" over all sets
of parameter sets, q, for which the risk metric equals f:

\[\]

where p(Y =f\| data) is the probability density at Y =f given the data
(i.e. the posterior probability density function for Y), L(q\|data) is
the probability of the data given the parameters q (i.e. the likelihood
function), *p*(q) is the initial assumption regarding the probability of
different true process parameters (i.e. the prior), and h(data) is a
normalizing constant. There are many texts on Bayesian statistics. The
sections in Hilborn and Mangel (1997) on Bayesian methods are
particularly accessible for ecological applications.

The following algorithm calculates using as the data. Similar algorithms
can be found in Hilborn and Mangel (1997), pages 256-260. The code to
run this algorithm and produce Figs. and is given in the digital
appendix.

[Step 1]{.underline}. Specify prior distributions for m, s^2^, and . For
this example, I used a uniform distribution of m on (-0.2, 0.2) and a
uniform distribution of on (0, 1). As will be seen, the data provide
much information on m and , and so the posterior distributions are quite
different than the priors. The same cannot be said for s^2^. My
subjective prior having examined hundreds of salmonid time series, is
that s^2^ is between 0.001 and 0.1. To express this, I used a prior for
s^2^ of a gamma distribution with scale of 0.25 and shape of 5. I
compared this with an alternate prior of a uniform distribution on (0,
0.5). This prior says the prior belief is that s^2^ \> 0.1 by a factor
of four to one, which seems unlikely. A standard reference prior for
normal variances is uniform on log s^2^ (cf. Lee 1989, section 2.7), but
for the corrupted diffusion model, *L*(s^2^ \| data) goes to a constant
as log s^2^ goes to negative infinity since s^2^ = 0 is perfectly
plausible from the model's standpoint. Thus this prior will lead to
infinite probability densities when Eqn. is integrated.

[Step 2]{.underline}. Randomly draw values of m, s^2^, and from their
priors.

[Step 3]{.underline}. Calculate the Y using the values of m, s^2^ from
step 2. If the risk of interest is the probability of 90% decline in a
given time frame, this means putting these parameters into Eqn. to
calculate Y. If the risk metric of interest is l, Y=exp(m).

[Step 4]{.underline}. Calculate the total likelihood, **L**, of the m,
s^2^, and parameters generated in step 2 given the actual estimates from
the data. Note that is independent of and . To simplify this step, I
also assumed that and are independent, given a particular set of
parameters. This is not entirely true. In simulations, I have found
correlation coefficients of approximately 0.2.

**L** =x x

The likelihoods are calculated from the probability of given the m,
s^2^, and parameters generated in step 2. The distributions of and are
given in Eqn. and from these the likelihoods are given. The equals the
probability density at of a normal with mean m and variance . The equals
the probability density at of a gamma with scale b and shape w where and
. Recall from Eqn. that E() . The equals the probability density at of a
gamma with scale b and shape w where and . Recall the introduction of
that E() . See any statistics text for the probability density functions
for normal and gamma distributions. The gamma distribution arises given
that , thus \~c^2^(*q*) E(s^2^)/*q*, which equals a gamma distribution
with scale = 2 E(s^2^)/*q* and shape = *q*/2.

[Step 5]{.underline}. Repeat steps 2-4 thousands of times, recording the
parameters (step 2), the risk values (step 3), and total likelihood
(step 4) for each iteration.

[Step 6]{.underline}. Divide the outputs of interest, e.g. risk values
and parameters, into discrete intervals and calculate the proportion of
the total likelihood (**L** summed over all the iterations) within each
interval. Divide the probabilities by the size of the interval. This
gives the posterior probability distribution for the output of interest.

Figure shows the posterior probability distributions for m, s^2^, and
relative to the priors for the 1960-1998 Methow River spring chinook
data. The plots (d-f) show the posterior distributions using the uniform
prior on s^2^. The posterior on s^2^ is not uniform and gives more
weight to s^2^ \< 0.2 than s^2^ \> 0.2. The plots (a-b) show the
posterior distributions using my quite informative prior on s^2^ based
on previous analyses. The posterior distribution for s^2\ ^is the same
as the prior in this case. Together, this indicates that the data are
not providing much information beyond that s^2^ is more likely to be
less than 0.2. For m and , the data provide much more information, and
the posteriors are quite different than the uniform priors. Putting a
strong prior on small s^2^ decreases the uncertainty in m and , and
their posterior distributions are more peaked, but even with a prior
that gives weight to higher s^2^, there is still considerable
information on m and .

Figure shows the posterior probability distributions for the following
risk metrics: l = exp(m), the probability that the population ever
experiences a 90% at some point in the future (this is p' in Eqn. ), and
the expected value of the probability of 90% decline within different
time frames. The expected value of the probability of extinction
multiplies the probability of extinction given a particular set of true
parameters by the probability of those parameters (specified by the
posterior probability distribution) and integrates this over all
possible parameter values. See Ludwig (1996a) for an extensive
discussion of the expected value of the probability of extinction, which
he terms the Bayesian probability. A simple algorithm to calculate the
expected probability is to divide 0 to 1 into regular intervals, take
the value of the posterior probability density of P(90% decline in *x*
years) from step 6 at the center of those intervals and multiply that by
the width of the interval, and then sum these together to calculate the
expected value of P(90% decline in *x* years).

Most of the area of posterior probability distribution of l is in the
region l \< 1 (Figs. a, d) indicating that the data give high support to
long-term declining dynamics (l \< 1). The support is considerably
higher when strong prior information on s^2^ is included (Fig. 9a), but
is still high when less specific prior information on s^2^ is used (Fig.
9d). Figures (b,e) show the posterior distributions of the probability
of eventual 90% decline. The prior (based on the parameter priors) gave
roughly equal weight to a 0 or 1 probability. The posterior distribution
was strongly skewed towards a probability of 1 for eventual 90% decline
indicating that the data strongly support a high probability of eventual
90% decline. This occurred because either m \< 0, which means 90%
decline is certain eventually, or if m \> 0, s^2^ was large since that's
the only way we could have observed a 7% yearly decline over the last 38
years by chance with m \> 0. Thus the probability of eventual 90%
decline was still high. Given that data strongly support eventual
probability of a 90% decline, when is it likely to occur? Figures (c,f)
show the expected value of the probability of 90% decline within
different time frames. Fifty percent of the population trajectories
generated by sampling from the posterior distributions of the parameters
would have declined 90% within 30 years. Over 65% percent would have
declined 90% within 50 years.

Overall, l appeared to be the most straightforward metric to use. It
behaved well and is not overly sensitive to estimation of s^2^, which is
imprecise. The l metric is also transparent. If l \< 1, the population
is declining; if l \> 0.9, its collapsing. If the population needs to be
recovered from low levels, l must be increased above 1.0. Risk metrics
involving probabilities are much more problematic. They are highly
variable and sensitive to estimation of s^2^. This section gave an
example of how to incorporate uncertainty into probability metrics and
how to express that uncertainty, however much work remains to be done in
this area. It is tempting to argue that because of their high
uncertainty, metrics involving extinction or quasi-extinction
probabilities are bad and should simply be abandoned; indeed others have
argued as much. However, a measure of extinction risk is integral to the
legal definition of "jeopardy" under the U.S. Endangered Species Act.
While this is true, research into how to quantify and express the
extinction risk in a way that properly characterizes what we do and do
not know about the probability of population declines to critical levels
is essential to science-based decision making for threatened and
endangered species.

# Concluding remarks

The diffusion approximation approach has a solid foundation on theory
concerning the behavior of stochastic age-structured models and opens a
large toolbox of quantitative methods for linear stochastic models. In
addition, the fact that many population processes with a multitude of
parameters can be approximated by a three parameter model means that
parameters can estimated with increased statistical power -- relative to
trying to estimate a large number of parameters with limited data.
However, the diffusion approximation approach makes simplifying
assumptions, which will always be violated to some degree for real
populations. This paper describes a method for evaluating the
appropriateness of the diffusion approximation for a particular PVA
application and for selecting among the currently available
parameterization methods. The salmonid populations modeled in this paper
exhibit substantial deviations from the diffusion model, such as
non-stable age structure, high temporal autocorrelation in ocean
survivorship and juvenile density dependence, yet the diffusion model
does approximate the behavior of the simulated population trajectories,
in terms of probabilities of crossing severe decline thresholds and of
the predicted distribution of future population sizes.

The analyses also illustrate that careful estimation of the diffusion
parameters is critical and that poor performance of estimators should
not be confused with poor performance of the diffusion approximation. In
particular, one should be cautious about assuming that non-process error
within the data is low even if sampling error is known to be low. In the
simulated salmonid populations, density-dependence in the juvenile stage
and perturbation of the age structure created feedbacks that led to
non-process error variability within the data. The result was high
non-process error variability even with no sampling error added. This
led to overestimation of the process error with parameterization methods
that are not designed to deal with high non-process error. Two
estimators designed to deal with high non-process error were tested: a
Kalman filter maximum-likelihood estimator and a slope estimator. For
the 20-year salmon time series with small process error, the majority of
Kalman filter maximum likelihood estimates of *s*^2^ were essentially
zero. In other words, a model with non-process error alone would have
the highest likelihood of producing the data. An alternative estimator,
the slope estimator for s^2^, performed better for the simulated salmon
populations and gave s^2\ ^estimates that were within 50-150% of the
true value. For low to medium sampling error levels relative to those
expected within actual salmon censuses (but still five to fifty times
higher than the expected process error), gave median estimates that did
not appear to substantially degrade quasi-extinction estimates, however
for high sampling error, s^2^ estimates were an order of magnitude too
high and this did degrade the estimates. The robustness of to
non-process error due to density dependence is encouraging; similar
robustness was also found in a cross-validation of
diffusion-approximation methods using a large collection of real salmon
time series (Holmes and Fagan 2002).

The purpose of this paper is to illustrate a method for investigating
the appropriateness of diffusion approximation for a specific population
and to evaluate its parameterization given the data constraints for that
population. New parameterization methods are being continually refined,
and no doubt better methods than those presented here will eventually be
available. The methods described here can be used to evaluate and
compare these new methods to existing ones. Experience suggests that the
best parameterization method is application specific; that is, it may
not be possible to find a method that is best for all situations. In
particular, the relative performances of the different parameterization
methods for the simulated salmonid time series in this paper should not
be over-generalized to their performance for other population time
series. The maximum likelihood estimates assuming zero non-process error
did especially poorly for the salmonid time series. However, a study of
the performance of these estimators on a large collection of real time
series found that they provided good estimates for time series of a wide
variety of species -- just not for salmon (Fagan et al. unpublished
manuscript). Presumably, the level of non-process error in many
non-salmonid time series is not so large as to significantly bias s^2^
estimates. Similarly, maximum likelihood estimation assuming non-zero
non-process error (the Kalman estimates) produced severe underestimates
of s^2^ for the twenty-year salmonid time series. However this was not a
problem with longer simulations (100+ years), and for non-salmonid
population processes, it may perform well for short time series. It is
also worth noting that the Kalman algorithm does estimate the parameters
with the highest likelihood of producing the data -- assuming all
parameters are equally likely. In the real world, all parameters are not
equally likely, specifically process error equaling essentially zero is
unlikely given that year-to-year variability in survivorship and
fecundity is occurring. Incorporating informative priors on s^2^ into
the Kalman algorithm may be a way to formally avoid unrealistic
parameter estimates.

Regardless of the PVA model used, a certain amount of variability in
estimated parameters and risk metrics is an unavoidable aspect of the
analysis of stochastic population processes, simply due to the nature of
these processes. One of the strengths of DA methods is that the
statistical distributions of the estimated parameters are known. As a
result, the uncertainty in the estimated risks can be calculated. This
is often not the case for other PVA approaches, such as Leslie matrix
models or individual-based simulations, where uncertainty in the
estimated model parameters is often poorly known, if known at all. Even
though the uncertainty in DA risk metrics can be calculated, this
uncertainty is definitely high. In this situation, examining either the
likelihood functions or the posterior probability distributions for the
risk metrics, rather than simply the point estimates and confidence
intervals, will help to clarify the level of data support for different
risks and to choose risk metrics that are most informative. Statistical
decision theory (e.g. Berger 1985) provides a framework for integrating
estimates of the data support for different risk levels with the
consequences of different true risk levels. Taking an approach that
assesses the degree of data support for questions of conservation
concern emphasizes that while uncertainty does exist, being uncertain
does not mean we know nothing.

# Acknowledgements

I thank R. Hinrichsen, C. Jordan, P. Kareiva, P. Levin, M. McClure, P.
McElhany, and B. Sanderson for helpful discussions throughout the
development of these methods. T. Cooney provided the impetus and initial
simulations that led to the development of the input corrections. T.
Cooney, P. Kareiva , and M. Marvier supplied the salmon matrices. I also
thank S. Heppell, B. Morris and an anonymous reviewer for reviews that
helped improve this paper, and J. Butzerin for her copy-editing.

**Literature Cited**

Achord, S., P. S. Levin, and R. W. Zabel. 2003. Density-dependent
mortality in Pacific salmon: the ghost of impacts past? Ecology Letters
6: 335-342.

Berger, J. O. 1985. Statistical decision theory and Bayesian analysis.
Springer Verlag, NY, NY, USA.

Caswell, H. 2001. Matrix Population Models. Sinauer, Sunderland, Mass.

Dennis, B., P. L. Munholland, and J. M. Scott. 1991. Estimation of
growth and extinction parameters for endangered species. Ecological
Monographs 61:115-143.

Dorazio, R. M., and F. A. Johnson. 2003. Bayesian inference and decision
theory ---a framework for decision making in natural resource
management. Ecological Applications 13:556--563.

Dunham, J. and B. Rieman. 2001. Sources and magnitude of sampling error
in redd counts for Bull Trout. North American Journal of Fisheries
Management 21:343-352.

Ellner, S.P., J. Fieberg, D. Ludwig and C. Wilcox. 2002. Precision of
population viability analysis. Conservation Biology 16: 258-261.

Fagan, W. F., J. Rango, A. Folarin, J. Sorensen, J. Lippe, and N. E.
McIntyre. Cross-validation of quasi-extinction risks from real time
series: an examination of diffusion approximation methods. Manuscript.

Fieberg, J. and S. Ellner. 2000. When is it meaningful to estimate
extinction probability? Ecology 81:2040-2047.

Gerber, L., D. DeMaster and P. Kareiva. 1999. Grey whales and the value
of monitoring data in implementing the U.S. Endangered Species Act.
Conservation Biology 13:1215-1219.

Harvey, A. C. 1989. Forecasting, structural time series models and the
Kalman filter. Cambridge University Press, Cambridge, UK.

Heyde, C.C. and J. E. Cohen. 1985. Confidence intervals for demographic
projections based on products of random matrices. Theoretical Population
Biology 27: 120-153.

Hinrichsen, R. A. 2002. The accuracy of alternative stochastic growth
rate estimates for salmon populations. Canadian Journal Fisheries and
Aquatic Sciences 59:1014--1023.

Hilborn, R. and M. Mangel. 1997. The ecological detective: confronting
models with data. Princeton University Press, Princeton, NJ.

Holmes, E. E. 2001. Estimating risks in declining populations with poor
data. Proceedings of the National Academy of Science USA 98:5072-77.

Holmes, E. E. and W. F. Fagan. 2002. Validating population viability
analysis for corrupted data sets. Ecology 83: 2379-2386.

Jones, E. L., T. J. Quinn, and B. W. Van Alen. 1998. Observer accuracy
and precision in aerial and foot survey counts of pink salmon in a
Southeast Alaska stream. North American Journal of Fisheries Management
18:832-846.

Kareiva, P., M. Marvier, and M. McClure. 2000. Recovery and management
options for spring/summer chinook salmon in the Columbia River Basin.
Science 290: 977-979.

Lande, R. and S. H. Orzack. 1988. Extinction dynamics of age-structured
populations in a fluctuating environment. Proceedings of the National
Academy of Science USA 85:7418-7421.

Lande, R., S. Engen, and B. Saether. 2003. Stochastic population models
in ecology and conservation: an introduction. Oxford University Press,
Oxford, UK.

Lee, P. M. 1989. Bayesian statistics: an introduction. Halsted Press,
N.Y., N.Y.

Lindley, S. T. 2003. Estimation of population growth and extinction
parameters from noisy data. Ecological Applications 13: 806-813.

Ludwig, D. 1996a. Uncertainty and the assessment of extinction
probabilities. Ecological Applications 6: 1067-1076.

Ludwig, D. 1996b. The distribution of population survival times.
American Naturalist 147: 506-526.

Ludwig, D. 1999. Is it meaningful to estimate extinction probabilities?
Ecology 80:298-310.

McClure, M. M., E. E. Holmes, B. L. Sanderson, and C. E. Jordan. 2003. A
large-scale, multi-species risk assessment: anadromous salmonids in the
Columbia River Basin. Ecological Applications 13:964-989.

Morris, W. F. and D. F. Doak. 2002. Quantitative conservation biology:
theory and practice of population viability analysis. Sinauer Press,
Sunderland, MA.

Morris, W., D. Doak, M. Groom, P. Kareiva, J. Fieberg, L. Gerber, P.
Murphy, and D. Thomson. 1999. A practical handbook for population
viability analysis. The Nature Conservancy.

Morris, W. F., P. L. Bloch, B. R. Hudgens, L. C. Moyle, and J. R.
Stinchcombe. 2002. The use of population viability analysis in
endangered species recovery planning: Past trends and recommendations
for future improvement. Ecological Applications 12: 708-712.

Nicholls, A. O., P. C. Viljoen, M. H. Knight and A. S. Van Jaarsveld.
1996. Evaluationg population persistence of censused and unmanaged
herbivore populations from the Kruger National Park, South Africa.
Biological Conservation 76:57-67.

Sabo, J. L., E. E. Holmes, and P. Kareiva. 2003. The efficacy of simple
viability models in ecological risk assessment: Does density dependence
matter? Ecology, *in press*.

Tuljapurkar, S. D. 1989. An uncertain life: demography in random
environments. Theoretical Population Biology 35:227-294.

Tuljapurkar, S. D. and S. H. Orzack. 1980. Population dynamics in
variable environments. I. Long-run growth rates and extinction.
Theoretical Population Biology 18:314-342.

Wade, P. R. 2000. Bayesian methods in conservation biology. Conservation
Biology 14: 1308-1316.

Wilcox, C. and H. Possingham. 2002. Do life history traits affect the
accuracy of diffusion approximations for mean time to extinction?
Ecological Applications 12: 1163-1179.

Williams, J. G., S. G. Smith, and W. D. Muir. 2001. Survival estimates
for downstream migrant yearling juvenile salmonids through the Snake and
Columbia River hydropower system, 1966-1980 and 1993-1999. North
American Journal of Fisheries Management 21: 310-317.

Table 1. Parameter estimates for the matrix models from (M. McClure,
National Marine Fisheries Service) and (T. Cooney, National Marine
Fisheries Service).

  ------------------------------------------------ ----------------------- -------------- -----------
  Parameter                                        spring/summer chinook   fall chinook   steelhead
                                                                                          
  *p*\*, mouth to spawning ground survival         0.4815                  0.36           0.58
                                                                                          
  *s*~1~, egg to age-1 survival                    0.018                   0.0044167      0.047
  *s*~2~, mean age-1 to age-2 survival in-stream   0.044                   NA             0.009
  *s*~o~, ocean survivorship                       0.8                     0.8            0.8
  *b*~2~, fraction of age 2 that spawn             0                       0              0.009
  *b*~3~, fraction of age 3 that spawn             0                       0.081          0.333
  *b*~4~, fraction of age 4 that spawn             0.216                   0.65           0.693
  *b*~5~, fraction of age 5 that spawn             1.0                     0.863          0.923
  *b*~6~, fraction of age 6 that spawn             NA                      1.0            1.0
  *h*~2~, in ocean harvest of age 2                NA                      0.0123         NA
  *h*~3~, in ocean harvest of age 3                NA                      0.0465         NA
  *h*~4~, in ocean harvest of age 4                NA                      0.1368         NA
  *h*~5~, in ocean harvest of age 5                NA                      0.1838         NA
  *h*~6~, in ocean harvest of age 6                NA                      0.1953         NA
  *m*, average ♀ eggs per ♀ spawner                2747                    1500           2500
                                                                                          
  ------------------------------------------------ ----------------------- -------------- -----------

\**p* = fraction that survive harvest during upstream migration ×
fraction that survive upstream migration to spawning ground × survival
on the spawning ground

Table 2. Overview of estimator performance and limitations.

  ------------------------------ ----------------------------- ----------------------------------------------------- -----------------------------------------
                                 Corrects for sampling error   Expected value                                        Number of years needed (rules of thumb)
                                                                                                                     
  m* estimators*                                                                                                     
                                                                                                                     
  Dennis (same as Heyde-Cohen)   no                            m, sensitive to age-structure perturbations in data   10+
  Kalman                         yes                           m                                                     10+
                                                                                                                     
  Runsum                         yes                           m                                                     15+
                                                                                                                     
  Slope                          yes                           m                                                     20+
                                                                                                                     
  s*^2^ estimators*                                                                                                  
  Dennis                         no                                                                                  10+
  Kalman                         yes                           (\*)                                                  50+ (\*)
  Slope                          yes                                                                                 20+
  Heyde-Cohen                    no                            for short series                                      100+
                                                                                                                     
  ------------------------------ ----------------------------- ----------------------------------------------------- -----------------------------------------

\* Although the mean estimates are equal to the true value, i.e. are
unbiased, the distribution of estimates can be highly skewed for short
time series with low s^2^ and high , such that the vast majority of
estimates are 0 and a few rare extremely high estimates bring the mean
up to s^2^.

**Figure Legends**

Figure . Behavior of the population trajectories from the age-structured
salmonid models versus the behavior of the diffusion approximation. For
a diffusion model, the mean and variance of log population size
(*N~t~*/*N~0~*) should both increase linearly with *t*. The
age-structured models had this property beginning at approximately *t*=5
(a-c). The diffusion model also correctly estimated the probability of
90% decline (*N~t~*/*N~0~ = *1/10) within different time frames (d-f).

Figure . Performance of the m estimates from the Dennis, runsum, Kalman,
and slope methods. One thousand 20-year spawner time series were
generated from the spring/summer chinook, fall chinook, and steelhead
stochastic matrix models. The column labels indicate the type of error
added: none = stable starting age structure and no sampling error, age =
no sampling error but perturbed age structure, low, med, high = stable
age structure plus lognormal sampling error with = 0.1, 0.5 or 0.75,
respectively. The box plots summarize the estimates among all 1,000
simulations. The line in the box shows the median estimates and the box
encloses the middle 75% of the estimates, the whiskers enclose 95% of
the estimates.

Figure . Performance of the *s*^2^ estimates from the Dennis,
Heyde-Cohen, Kalman and slope methods for the spring/summer chinook,
fall chinook, and steelhead simulations. See Fig. for details.

Figure . Error in *m* and *s*^2^ estimates due to yearly hatchery inputs
into the salmon simulations. The starting age structure was perturbed
and low sampling error ( was added. The parameter m was estimated using
with either no input correction (Eqn. ) or with input correction (Eqn.
). The s^2^ parameter was estimated using versus with no explicit
hatchery correction.

Figure . Variability of l estimates using . See Fig. for simulation
details.

Figure . Estimated versus actual probability of 90% decline within a
given time horizon. Probability estimates were calculated using Eqn.
with and ('Dennis'), Eqn. with and ('Kalman'), or Eqn. with and
('Runsum/Slope'). The parameters were estimated from simulated spawner
counts. The true probability was calculated from one thousand 100-year
time series generated from each stochastic matrix model. Results are
shown for simulations with no sampling error added to the spawner
counts. Simulations were started with a stable initial age-structure.
The solid line shows the median errors observed in the simulations and
the dashed and dotted lines encompass 75% and 95% of the estimates
respectively.

Figure . Effect of sampling error on estimated versus actual probability
of 90% decline within a given time horizon. Probabilities were
calculated using Eqn. with and . Results are shown for simulations with
low, medium or high sampling error added to the spawner counts.
Simulations were started with a stable initial age-structure. The solid
line shows the median errors observed in the simulations and the dashed
and dotted lines encompass 75% and 95% of the estimates respectively.

Figure . Estimated posterior probability distributions and prior for m,
s^2^, and . This analysis used a 1960-1998 spawner count time series
from Methow River spring chinook in the Upper Columbia River basin. The
right panels show the estimated distributions using an informative
prior, which put most probability on s^2^ between 0 and 0.1. The left
panels show the estimated distributions using a uniform prior which put
even probability on s^2^ between 0 and 0.5.

Figure . Estimated posterior probability distributions and prior for l,
the probability that the population ever experiences a 90% decline at
some point in the future, and the expected value of the probability of
90% decline within different time frames.