**[Viability analysis for endangered metapopulations: a diffusion
approximation approach]{.underline}**

[E. E. Holmes]{.underline}^[1]{.underline}^

[National Marine Fisheries Service]{.underline}

[Northwest Fisheries Science Center]{.underline}

[2725 Montlake Blvd E]{.underline}

[Seattle, WA 98112]{.underline}

[B. Semmens]{.underline}

[Zoology Department]{.underline}

[University of Washinton]{.underline}

[Seattle, WA 98195]{.underline}

[1 Corresponding author]{.underline}

[e-mail: ]{.underline}

[phn: 206-860-3369]{.underline}

[fax: 206-860-3467]{.underline}

# [23.1. Introduction]{.underline}

[Population viability analysis (PVA) assesses the rate of population
decline and the risks of extinction or quasi-extinction over a defined
time horizon for a population of concern (Gilpin and Soule 1986; Boyce
1992; Morris and Doak 2002). Although the techniques employed to conduct
PVA are varied, they typically involve building quantitative models that
are parameterized by demographic and environmental data. PVA was first
used in the early 1980s (Shaffer 1981), and in the past decade it has
gained broad acceptance in the conservation community as a useful tool
for assessing and managing 'at risk' species (Beissinger 2002; Morris
and Doak 2002; Reed et al. 2002). This is particularly true for
demogaphic PVAs, due in large part to the advancements in monte carlo
techniques and desktop computers (Beissinger 2002). The International
Union for the Conservation of Nature (IUCN)'s Red List Criteria,
probably the most widely applied set of decision rules for determining
the status of at risk species, is partially defined by metrics that
require some form of PVA (IUCN 1994). For instance under one of the Red
List criteria, a taxon may be classified as endangered if a "reduction
of at least 50%, projected or suspected to be met within the next ten
years or three generations" is predicted.]{.underline}

[Although many PVAs are focused on single populations in single sites,
there are often needs for spatially explicit PVAs: many populations of
conservation concern are distributed across multiple sites and
additionally, the primary anthropogenic threats facing at risk species
are habitat destruction and alteration, which are fundamentally spatial
processes (Wilcove et al. 1998). Several software packages have been
written for spatially explicit PVA, including
RAMAS]{.underline}^[®]{.underline}[\ ]{.underline}^[Metapop (Akçakaya
1997) and RAMAS]{.underline}^[®]{.underline}^[ GIS (Boyce 1996), ALEX
(Possingham and Davies 1995), and VORTEX (Lacy 1993). These models
incorporate a diversity of demographic and spatial attributes such as
distance-dependent migration, allee effects, social population
structure, habitat quality and spatial arrangement, and genetic
variability. The development of flexible sophisticated PVA software
packages such as these has made the construction and simulation of
spatially explicit PVA models feasible for those who are not highly
skilled programmers and has greatly increased the number of managers and
scientists capable of using spatially realistic PVA models.
]{.underline}

[As the use of PVA has grown in conservation science, so have concerns
that PVAs are often overextended given limited data sets (Reed et al.
2002). Beissinger and Westpahl (1998) suggested that PVA should be
limited to assessing relative risks over short time frames using the
simplest models that can reasonably be justified. For single species
with spatially simple structure, data needs can often be met when
Beissinger and Westpahl's call for model moderation and simplicity are
heeded. When one is faced with species with more complex spatial
structure, a much larger amount of data is needed to parameterize the
dynamics of individual local populations, the levels and patterns of
dispersal, and the spatial pattern of temporal correlations among local
populations (e.g. Ralls et al. 2002). Unfortunately, collection of the
data needed to parameterize a spatial model is rare for species of
conservation concern, at least in the United States (Morris et al.
2002), and there is a disconnect between the parameter requirements for
spatially explicit PVA models, and the willingness and/or ability of
management agencies to collect the types of data needed to appropriately
apply such tools. Since it is usually impossible to retroactively
fulfill data requirements for a spatial PVA and there will always be
cases where collection of spatial data is infeasible, managers require
PVA tools that that can help guide conservation of metapopulations in
the absence of large amounts of spatial data.]{.underline}

#### [Diffusion approximation for metapopulations]{.underline}

[One approach to the problem of limited population data is to find a
diffusion approximation that correctly models the long-run statistical
properties of a complex population process. This approach has been used
successfully for single population models (Karlin and Taylor 1981; Lande
and Orzack 1988; Lande 1993; Dennis et al. 1991; Hill et al. 2002; see
also Morris and Doak 2002, Ch. 3, and Lande et al. 2003) and reduces the
problem of parameterizing a large model with many parameters to the much
simpler task of parameterizing a two-parameter diffusion model. One of
the main practical implications of the diffusion approximation approach
is that it is not necessary to know the multitude of parameters
describing the local dynamics, dispersal levels, spatial patterns of
dispersal, and spatial synchrony between local populations in order to
make basic predictions about the statistical distribution of the
long-term metapopulation or local population trajectories. The relevant
two parameters for the diffusion approximation can be estimated from a
simple time series of counts from the population process. ]{.underline}

[In this chapter the diffusion approximation approach is used to model
the long-run behavior of spatially structured populations. Our focus is
on stochastic metapopulations characterized by structured population
size, density-independent local dynamics, and, in keeping with the
assumption of density independence, a metapopulation that is declining
as a whole. Local populations are assumed to have patch-specific
structured local dynamics and dispersal rates, with spatial structure
among local populations both in terms of their local dynamics and
dispersal patterns. Description of the long-run statistical behavior of
the metapopulation trajectories using a diffusion approximation allows
the estimation of PVA risk metrics such as the long-term rate of
metapopulation decline and the probability of reaching different
threshold declines over different time horizons (i.e. probabilities of
extinction or quasi-extinction). These methods for estimating
metapopulation PVA metrics are illustrated using data from two chinook
salmon metapopulations in the U.S. Pacific Northwest. ]{.underline}

# [23.2. A stochastic metapopulation model]{.underline}

Our focus is on declining metapopulations, and thus what has been termed
non-equilibrium metapopulations. We model a collection of local
populations connected by dispersal where local populations have
density-independent local dynamics, which may be "sources" or "sinks",
but the metapopulation as a whole is declining. Dispersal levels could
very low, resulting in basically independent local populations, or
extremely high, resulting in essentially one population. From a
practical standpoint, this approach is most appropriate when dispersal
is not insignificant (e.g. above 2% per year localized dispersal or 0.1%
global dispersal), otherwise parameterization of the model requires
inordinately long time series. Data from this type of metapopulation
would be characterized by fluctuating local population trajectories, but
actual extinctions would be unusual until the metapopulation has very
few individuals. Our model assumes no density-dependence nor carrying
capacities within the individual local populations. Such a model is only
appropriate in cases where the population is declining and all local
populations are well below their carrying capacities. Our example using
data on chinook salmon illustrates a situation that is likely to be well
modeled as this type of metapopulation.

[In the following section, a rather parameter intensive mathematical
description of a stochastic, declining metapopulation is given. However,
the reader should keep in mind that this model will not be
parameterized. Rather the asymptotic behavior of this model's
trajectories will be derived and that information will be used to
develop a diffusion approximation of the process. Time series data will
then used to parameterize the diffusion approximation.]{.underline}

*[The model]{.underline}*

Consider an individual local population *i* with stochastic yearly
growth and stochastic dispersal to and from other local populations.
Such a local population's numbers in year *t*, *N~i~(t)*, could be
described as follows:

[ \[\]]{.underline}

[where
]{.underline}*[z]{.underline}~[i]{.underline}~*[(]{.underline}*[t]{.underline}*[)
is the stochastic growth rate of local population
]{.underline}*[i]{.underline}*[ in year ]{.underline}*[t]{.underline}*[
and is a random variable with some unspecified statistical distribution
with mean ]{.underline}*[m]{.underline}~[i]{.underline}~*[ and variance
. The ]{.underline}*[m]{.underline}~[i]{.underline}~*[ term will be
referred to as the local population's intrinsic growth rate; it will not
be observed since the local population is subject to immigration and
emigration. Some fraction of individuals,
]{.underline}*[d]{.underline}~[i]{.underline}~*[(]{.underline}*[t]{.underline}*[),
leaves local population ]{.underline}*[i]{.underline}*[ at year
]{.underline}*[t]{.underline}*[ and disperses to other local
populations, and dispersal into local population
]{.underline}*[i]{.underline}*[ occurs from other local populations. The
fraction of dispersers from local population
]{.underline}*[j]{.underline}*[ that go to local population
]{.underline}*[i]{.underline}*[ in year ]{.underline}*[t]{.underline}*[
is
]{.underline}*[a]{.underline}~[ji]{.underline}~*[(]{.underline}*[t]{.underline}*[)]{.underline}~[\ ]{.underline}~[and
can vary depending on the destination, ]{.underline}*[i]{.underline}*[,
thus allowing for spatially-structured dispersal. The dispersal
parameters,
]{.underline}*[d]{.underline}~[i]{.underline}~*[(]{.underline}*[t]{.underline}*[),
and
]{.underline}*[a]{.underline}~[ji]{.underline}~*[(]{.underline}*[t]{.underline}*[),
are assumed to be temporally random variables from some unspecified
statistical distribution.]{.underline}

#### [The model in matrix form]{.underline}

The model for the entire metapopulation can be written using a random
transition matrix, **A**(*t*), which encapsulates both the dispersal and
local growth:

[ \[\]]{.underline}

[where]{.underline}

\[\]

[The '(]{.underline}*[t]{.underline}*[)' on the
]{.underline}*[d]{.underline}*['s, ]{.underline}*[a]{.underline}*['s,
and ]{.underline}*[z]{.underline}*['s have been left off to remove
clutter. There may be any level or spatial pattern of temporal
correlation between the intrinsic local growth rates,
]{.underline}*[z]{.underline}~[i]{.underline}~*['s, dispersal rates,
]{.underline}*[d]{.underline}~[i]{.underline}~*['s, and dispersal
patterns,
]{.underline}*[a]{.underline}~[ji]{.underline}~*['s.]{.underline}

[In the matrix model, each row represents 1 unit of habitat. Local
populations with multiple units of habitat appear as multiple rows with
very high dispersal between the units of habitat in that local
population. The habitat units within a local population could vary in
quality (i.e. habitat within a local population need not be uniform),
and different local populations certainly differ in the number of
habitat units they contain. The
]{.underline}*[d]{.underline}~[i]{.underline}~*['s and
]{.underline}*[a]{.underline}~[ji]{.underline}~*['s are assumed to be
drawn from some distribution that can be different for each local
population or local population pair. Although the
]{.underline}*[d]{.underline}~[i]{.underline}~*['s,
]{.underline}*[a]{.underline}~[ji]{.underline}~*~[,]{.underline}~[ and
]{.underline}*[z]{.underline}~[i]{.underline}~*['s are temporally random
variables, they are assumed to be stationary, i.e. that there is no
overall change in the mean values over time. For the purposes of this
paper, it will be assumed that the
]{.underline}*[d]{.underline}~[i]{.underline}~*['s,
]{.underline}*[a]{.underline}~[ji]{.underline}~*['s]{.underline}~[,]{.underline}~[
and ]{.underline}*[z]{.underline}~[i]{.underline}~*['s are all strictly
postitive, which this means that all local populations are connected to
each other to some (although possibly very low) degree and that mean
yearly geometric growth rates,
exp(]{.underline}*[m]{.underline}~[i]{.underline}~*[)'s, while possibly
very small are not zero. These assumptions imply that the
]{.underline}**[A]{.underline}**[(]{.underline}*[t]{.underline}*[)
describe an ergodic set of matrices (Caswell 2001, page 375). The
assumption of strict positivity is not strictly necessary. It is
possible for
]{.underline}**[A]{.underline}**[(]{.underline}*[t]{.underline}*[) to
describe an ergodic set if some elements of
]{.underline}**[A]{.underline}**[ are zero; it depends on the pattern of
zeros within ]{.underline}**[A]{.underline}**[ (cf. Caswell 2001 for a
discussion of the conditions under which matrices are
ergodic).]{.underline}

[The model is very general, allowing some sites to be dispersal sources
and others to be dispersal targets, allowing any spatial pattern of
dispersal or spatially correlated local growth rates, allowing any
pattern of temporal correlation amongst local growth rates and allowing
any combination or pattern of habitat sizes of local sites.]{.underline}

#### [Using random theory to understand the model's statistical behavior]{.underline}

Together, Eqns. and describe a quite generic model of a declining
metapopulation with density-independent local dynamics. From a viability
analysis perspective, one might ask the question: "Can one predict the
viability of the total metapopulation?" In more precise terms, this is
asking, what are the statistical properties of the metapopulation
trajectories of this type of connected collection of local populations
(of the form in Eqns. and )? Clearly, the matrix **A**(*t*) has a large
number of parameters that would be difficult, if not impossible, to
estimate for any given metapopulation of conservation concern. However,
using random theory, it can be shown that the long-term dynamics can be
described by only two parameters and that it is unnecessary to know the
multitude of other parameters for the purpose of projecting long-run
dynamics.

[To use this random theory, we first need to recognize that this
stochastic metapopulation model falls into the class of random processes
that involve products of ergodic random matrices, in this case products
of ]{.underline}**[A]{.underline}**[(]{.underline}*[t]{.underline}*[),
which can be seen by using Eqn. to project the vector of local
population sizes forward:]{.underline}

[ \[\]]{.underline}

where **N**(*t*) is the column vector of *N~i~*'s at time *t*, in Eqn. .
Products of random ergodic matrices have a well-established theoretical
foundation and have certain well-studied asymptotic statistical
properties. A brief review of two of the key results from this theory is
provided in Box 23.1 and a simulated example is shown to illustrate
these results. As described in Box 23.1, the theory demonstrates that
this stochastic, density-independent, metapopulation will have an
asymptotic growth rate and that the metapopulation, , the individual
*N~i~*(*t*)'s and sets of *N~i~*(*t*)'s representing the units of
habitat comprising a semi-independent local population will be
lognormally distributed with the same parameters:

[ \[\]]{.underline}

[Figure 23. shows an example of this behavior. A metapopulation is
simulated (described in Box 23.1) and over time, the metapopulation
declines at a constant rate and all
]{.underline}*[N]{.underline}~[i]{.underline}~*[(]{.underline}*[t]{.underline}*[)'s
have the same long-term fate. When viewed over short time frames,
]{.underline}*[t]{.underline}*[ small in the figure, the local sites
show different growth rates with some declining more or less than the
long-term rate, but over the long-term their rate of decline is the
same.]{.underline}

[ The model studied here approximates the local dynamics by a simple
exponential growth (or decline) model. However, it has been shown that
the results from random theory (presented in Box 23.1) also apply to a
more complicated metapopulation model where local dynamics are described
by stochastic age-structured Leslie matrices (Sanz and Bravo de la Parra
1998, pgs. 286-294). Essentially, this occurs because even when the
local dynamics are described by a local matrix model, the system can
still be described by products of random matrices.]{.underline}

**[23.3. Diffusion approximation]{.underline}**

[The asymptotic distribution of log
]{.underline}*[M]{.underline}*[(]{.underline}*[t]{.underline}*[) in Eqn.
has the same properties as the distribution of a diffusion process with
drift; it is normal and the mean and variance of the distribution of log
]{.underline}*[M]{.underline}*[(]{.underline}*[t]{.underline}*[)
increase linearly with time, ]{.underline}*[t]{.underline}*[. This
observation in the context of age-structured matrix population models
(Lande and Orzack 1988; Dennis et al. 1991) led to the use of a
diffusion approximation to enable parameterization using simple time
series and to enable calculation of extinction probabilities. Diffusion
approximation methods for single population populations are an important
and established method for approximating stochastic trajectories (Lande
and Orzack 1988; Dennis et al. 1991; chapter 3 in Morris and Doak 2002;
chapter 5 in Lande et al. 2003). The models for single populations are
mathematically analogous to the models used here for metapopulations
with a stochastic process involving products of random matrices.
However, in single population models, the matrix represents a life
history matrix rather than a growth and dispersal matrix, and the
]{.underline}**[N]{.underline}**[(]{.underline}*[t]{.underline}*[)
vector (in Eqn. ) represents different age or stage classes while in the
metapopulation matrix, it represents different local sites and
populations. ]{.underline}

[ A diffusion approximation with drift is a stochastic process with the
following properties (cf. Karlin and Taylor 1981, chap. 7):]{.underline}

[ \[\]]{.underline}

[For any non-overlapping pair of time periods,
]{.underline}*[t]{.underline}*~[1]{.underline}~[ \<
]{.underline}*[t]{.underline}*~[2]{.underline}~[ and
]{.underline}*[t]{.underline}*~[3]{.underline}~[ \<
]{.underline}*[t]{.underline}*~[4]{.underline}~[,
]{.underline}*[X]{.underline}*[(]{.underline}*[t]{.underline}*~[2]{.underline}~[)-]{.underline}*[X]{.underline}*[(]{.underline}*[t]{.underline}*~[1]{.underline}~[)
and
]{.underline}*[X]{.underline}*[(]{.underline}*[t]{.underline}*~[4]{.underline}~[)-]{.underline}*[X]{.underline}*[(]{.underline}*[t]{.underline}*~[3]{.underline}~[)
are independent random variables.
]{.underline}*[X]{.underline}*[(]{.underline}*[t+]{.underline}[t]{.underline}*[)
is a random variable with distribution . Correspondingly the probability
density function for
]{.underline}*[X]{.underline}*[(]{.underline}*[t]{.underline}*[+]{.underline}*[t]{.underline}*[)
given log
]{.underline}*[X]{.underline}*[(]{.underline}*[t]{.underline}*[)
is]{.underline}

[ \[\]]{.underline}

*[Behavior of metapopulation trajectories versus diffusion
trajectories]{.underline}*

[The diffusion approximation is based on the behavior of log
]{.underline}*[M]{.underline}*[(]{.underline}*[t]{.underline}*[) as
]{.underline}*[t]{.underline}*[ goes to infinity, however in PVA
settings the time frame of interest is substantially less than infinity
and typically in the range of 25 to 100 years. How well does the
diffusion approximation do over these finite time periods? To explore
this, a collection of 50 local populations were simulated that were
connected by global dispersal ranging from 0.1 to 5% per year and that
had correlated local dynamics,
]{.underline}*[z]{.underline}~[i]{.underline}~*[(]{.underline}*[t]{.underline}*[),
drawn from a Normal(mean = -0.05, variance = 0.09) and a temporal
covariance of (0.2)(0.09) between the
]{.underline}*[z]{.underline}~[i]{.underline}~*[(]{.underline}*[t]{.underline}*[)'s
of local populations in any given year. ]{.underline}

[If the log metapopulation trajectories behave like a diffusion process,
then if we repeatedly generate a large sample of replicate
metapopulation trajectories, the mean and variance of from those
trajectories, should be a constants over the time period of interest.
Additionally, should be normally distributed. To examine whether the
metapopulation trajectories had these properties, the simulations were
started from a distribution of local population sizes selected from the
equilibrium set of local population distributions and then run forward
for 200 years. This was repeated (using the same initial distribution of
local populations) 1000 times to estimate the distribution of . This
process was repeated for four randomly chosen initial distributions of
local population sizes. The mean and variance of is denoted as ,
respectively, in the figure and discussion below.]{.underline}

[Figure 23. illustrates the results. For dispersal levels 1% or higher,
the trajectories behaved like a diffusion process with roughly constant
and the distributions approximately normal according to a
Kolmogorov-Smirnov test at P = 0.05. For low dispersal, 0.1%, the
trajectories did not behave like a diffusion process for
]{.underline}*[t]{.underline}*[ less than 200 at least. Neither nor were
constant except for ]{.underline}*[t]{.underline}*[ \> 150 and the
normality assumption was generally violated except again at large
]{.underline}*[t]{.underline}*[. This means that when dispersal is very
low, diffusion approximations for this metapopulation would more
approximate than for metapopulations with higher dispersal.]{.underline}

[ Figure 23. illustrates results from one particular model. Repeating
this process for a number of different models indicated some general
behaviors. The higher the dispersal levels, the more trajectories
behaved like a diffusion process. Global dispersal levels of at least 2
to 5% were generally high enough to result in diffusion-like behavior
within a short time frame. Note that localized dispersal has the effect
of lowering the effective dispersal rates. The higher the amount of
temporal covariance between local populations in terms of their yearly
growth rates, the more the trajectories behaved like a diffusion
process. The simulations were done with the local population sizes
within the equilibrium set of local population distributions -- indeed
the theory is predicated on the local populations being near
equilibrium. For metapopulations with 2 to 5% dispersal, the local
populations equilibrated fairly quickly starting from all local
populations with equal numbers. However, at very low dispersal,
equilibration took thousands of time steps. This suggests that the
assumption of equilibrium should be viewed cautiously for
metapopulations that have very low dispersal rates between local
populations.]{.underline}

**23.4. Estimating the parameters**

Maximum likelihood estimates of and can be calculated using the
diffusion approximation for log *M*(*t*). Denote the observed time
series as *M* = *M*(0), *M*(1), *M*(2),..., *M*(*n*). If we approximate
log *M*(*t*) as a diffusion process, the likelihood function is given by
the product of the probability function distributions for the
transitions from log *M*(*t*+1) to log *M*(*t*), which is Eqn. with *t*
= 1, over *t* = 0, 1, 2, ..., *n*. Thus the log likelihood function is:

[ \[\]]{.underline}

[The maximum likelihood estimates are obtained by solving for the and
which maximize Eqn. :]{.underline}

[ \[\]]{.underline}

[Note that the unbiased estimator for would use
(]{.underline}*[n]{.underline}*[-1) rather than
]{.underline}*[n]{.underline}*[. The and are analogous to the estimates
of mean and variance from ]{.underline}*[n]{.underline}*[ samples from a
normal distribution, and confidence intervals on
]{.underline}*[m]{.underline}~[m]{.underline}~*[ and are
analogous:]{.underline}

[ \[\]]{.underline}

[where
]{.underline}*[t]{.underline}~[a]{.underline}~~[,q]{.underline}~*[ is
the critical value of a ]{.underline}*[t]{.underline}*[ distribution at
P = a and ]{.underline}*[q]{.underline}*[ degrees of freedom and is the
critical value of a chi-square distribution at P = a and
]{.underline}*[q]{.underline}*[ degrees of freedom. See Dennis et al.
(1991), page 120, for a more in-depth discussion of maximum likelihood
estimates for diffusion processes. Following Dennis et al.'s monograph,
parameter estimation based on the diffusion approximation has been
widely used for the analysis of single population trajectories. For a
discussion of parameter estimation that is not based on the diffusion
approximation, the reader is referred to Heyde and Cohen (1985).
]{.underline}

[The maximum likelihood estimates assume that the metapopulation has
reached a stochastic equilibrium and thus that the diffusion
approximation is reasonable. When exploring these methods using
simulations, it's important to allow the system to equilibrate, after
starting the simulation with something peculiar like all local
populations at the same size. Equilibruim can be monitored by waiting
for the variance of
(log(]{.underline}*[N]{.underline}*[(]{.underline}*[t]{.underline}*[)]{.underline}*~[i]{.underline}~*[)-log\[mean(]{.underline}*[N]{.underline}~[i]{.underline}~*[(]{.underline}*[t]{.underline}*[))\])
to stabilize. In simulations done for this paper, the distribution
stabilized relatively quickly when dispersal was non-zero. If dispersal
is zero, however, the distribution never stabilizes and the variance of
(log(]{.underline}*[N]{.underline}~[i]{.underline}~*[(]{.underline}*[t]{.underline}*[))-log\[mean(]{.underline}*[N]{.underline}~[i]{.underline}~*[(]{.underline}*[t]{.underline}*[))\])
continually increases. For an actual metapopulation, for which one wants
to conduct a PVA, it is also critical to test the appropriateness of the
diffusion approximation for one's time series data. Dennis et al.
(1991), page 121-122, and Morris and Doak (2002), pages 73-79, review
how to do this which is based on diagnostic procedures for evaluating
the appropriateness of linear models.]{.underline}

#### [Parameter bias]{.underline}

The estimators are unbiased maximum likelihood estimators for the
diffusion approximation, *X*(*t*). It is important to understand whether
and how these estimates are biased when working with short time series
of metapopulation trajectories, *M*(*t*), as opposed to an actual
diffusion process. In particular, is certain to be biased to some degree
since it relies on the diffusion approximation holding for *t* = 1 in
log *M*(*t+t*)/*M*(*t*), regardless of the length of the time series
used for estimation. This is not the case for, , which is also an
unbiased predictor for *M*(*t*) given a long time series (Heyde and
Cohen 1985).

[To numerically explore parameter bias from short time series,
simulations were used to look at the difference between and from a
20-year time series versus their true values
]{.underline}*[m]{.underline}~[m]{.underline}~*[ and . An example
metapopulation of 50 local sites was simulated with global dispersal and
correlated local growth rates,
]{.underline}*[z]{.underline}~[i]{.underline}~*[(]{.underline}*[t]{.underline}*[),
drawn yearly from a normal distribution with mean =
]{.underline}*[m]{.underline}~[i]{.underline}~*[, variance = f and
covariance of 0.2\*f between any two local growth rates Two versions of
the simulation were run; one to model uniform site quality (spatially
uniform ]{.underline}*[m]{.underline}~[i\ ]{.underline}~*[= -0.05) and
one to model highly variable site quality (spatially variable
]{.underline}*[m]{.underline}~[i]{.underline}~*['s). To explore biases
over a range of different dispersal and variability levels, the models
were run with dispersal between 0.1 and 5% per year and local
variability, f, between 0.1 and 0.5. These parameters translated to
metapopulation-level rates,
]{.underline}*[m]{.underline}~[m]{.underline}~*[, in the range of 0.01
to --0.05 and metapopulation-level variability, , in the range of 0.001
to 0.08. For each dispersal and local variability pair, 1000 replicate
metapopulation trajectories were simulated, each with an initial
distribution of local population sizes selected randomly from the
equilibrium set.]{.underline}

[The mean difference between and over the dispersal and local
variability parameter space was very low, \< 0.0015, for both the
uniform and variable ]{.underline}*[m]{.underline}~[i]{.underline}~*[
simulations. Overall the lack of bias in supports metrics that rely
primarily on this parameter, such as the metapopulation
]{.underline}*[l]{.underline}*[ (next section). For most of the
parameter space explored, , representing a 0 to 20% under- or
overestimation of . Larger biases, , representing a \>20% under- or
overestimation, were seen for some parameter combinations. The impact of
this bias depended on where was used. For instance, the effect on the
estimated confidence intervals on (Eqn. ) was minimal with the width of
the interval changing by a median 0.002. The effect on estimated passage
probabilities was higher although not dramatic. For example, the
estimated probability that the metapopulation will be ≤10% of current
levels at the end of 50 years was decreased by 0 to 0.04 (on a scale
from 0 to 1) for the uniform
]{.underline}*[m]{.underline}~[i]{.underline}~[
]{.underline}*[simulation and increased by 0 to 0.04 for the variable
]{.underline}*[m]{.underline}~[i]{.underline}~*[ simulation. The
estimated probability that the metapopulation will pass below 10% of
current levels at any point during the next 50 years was changed by 0 to
0.09. Overall, the effect of bias was low in these simulations, but this
will depend on the particular metapopulation and will need to be
investigated for individual cases of interest.]{.underline}

[In practical applications, one must contend with other factors which
can lead to parameter bias, but which are outside the scope of this
paper. In particular observation error, non-equilibrium local population
distributions and temporal autocorrelation can lead to parameter bias.
Such problems are being studied in the context single population PVA.
Much of this work is likely to be relevant for metapopulation PVA. See
Morris and Doak (2003), chapters 3-5, for a review and discussion of
current work in this area. ]{.underline}

#### [23.5. Metapopulation viability metrics]{.underline}

One of the most basic viability metrics is the long-term geometric rate
of decline (or growth) of a population, termed generally *l* in the PVA
literature. If *l* is less than 1.0, the population ultimately declines
to extinction and 100\*(1-*l*) is roughly the average yearly percent
decline. The metapopulation *l* is exp(*m~m~*) and its estimate is then

[. \[\]]{.underline}

[This definition of ]{.underline}*[l]{.underline}*[ follows Caswell's
use of the symbol ]{.underline}*[l]{.underline}*~[s]{.underline}~[ as
the long-term average stochastic growth rate:
]{.underline}*[l]{.underline}~[s]{.underline}~*[ =
(]{.underline}*[N]{.underline}*[(]{.underline}*[t]{.underline}*[)/]{.underline}*[N]{.underline}*[(0))]{.underline}^[1/]{.underline}^*^[t]{.underline}^*[
as (Caswell 2001, Eqn. 14.44). This is the long-run geometric growth
rate that would be observed in almost every trajectory. Defined this
way, if ]{.underline}*[l]{.underline}*[ \< 1, the population goes
extinct with certainty, eventually. This differs from Dennis et al.'s
use of the symbol ]{.underline}*[l]{.underline}*[ where
]{.underline}*[l]{.underline}*[ is used for and the long-term average
geometric growth rate is instead denoted by
]{.underline}*[a]{.underline}*[ = exp(]{.underline}*[m]{.underline}*[).
The maximum likelihood estimate of ]{.underline}*[l]{.underline}*[ is a
biased estimator; since is normally distributed, the median value of is
but the mean value is not. Dennis et al. (1991, Eqn. 79) give an
unbiased estimator (mean() = ]{.underline}*[l]{.underline}*[) based on
Shimizu and Iwase (1981), although Dennis et al. found negligible
differences between the biased and unbiased estimators in their
examples.]{.underline}

[From the asymptotic distribution of log
]{.underline}*[M]{.underline}*[(]{.underline}*[t]{.underline}*[), Eqn. ,
the probability that the metapopulation is below a threshold
]{.underline}*[b]{.underline}*[ at the end of
]{.underline}*[y]{.underline}*[ years can be calculated:]{.underline}

[ \[\]]{.underline}

Although this uses the asymptotic distribution, this is mitigated by the
fact that is used for the distribution at the end of *y* years but not
at any time before that. The estimate of replaces and by their estimates
and . Like the estimate of *l*, the median estimate of is equal to the
true value, but not the mean.

[Some metapopulations can have a low long-term risk of being below a
threshold due to a ]{.underline}*[l]{.underline}*[ near 1.0, but high
short-term risks of hitting that threshold due to high variability. Such
quasi-extinction or extinction probabilities are commonly used and very
important PVA metrics. The diffusion approximation for log
]{.underline}*[M]{.underline}*[(]{.underline}*[t]{.underline}*[) can be
used to estimate these probabilities for the metapopulation. The
probability of that the diffusion process,
]{.underline}*[X]{.underline}*[(]{.underline}*[t]{.underline}*[),
experiences a decline below a threshold log
]{.underline}*[b]{.underline}*[ at some time
]{.underline}*[T]{.underline}*[ less than
]{.underline}*[y]{.underline}*[ years is calculated by integrating over
the probability density function for first passage times for a diffusion
process with drift (Karlin and Taylor 1981, Theorem 5.3). Lande and
Orzack (1988), page 7419, goes through the calculation, which leads
to:]{.underline}

\[\]

[F( ) is the cumulative distribution function for a standard Normal(mean
= 0, variance =1). The estimate of for the metapopulation uses and with
]{.underline}*[X]{.underline}*[(0) = log
]{.underline}*[M]{.underline}*[(0). The estimated probability of
extinction (to 1 individual) is calculated using Eqn. and setting
]{.underline}*[b]{.underline}*[ equal to 1. The reader is cautioned that
estimates of extinction are problematic and estimates of
quasi-extinction (e.g. some threshold greater than 1 individual) are
more robust (cf. Morris and Doak 2002, pages 43-44). Also Eqn. uses the
diffusion approximation over short time scales since it calculates the
probability of hitting a threshold at any time, including short times,
before ]{.underline}*[y]{.underline}*[ years. This makes Eqn. more
approximate than other metrics.]{.underline}

[Other viability metrics based on the diffusion approximation such as
the mean time to extinction and median time to extinction are discussed
in Lande and Orzack (1988) and Dennis et al. (1991). Karlin and Taylor
(1981) ]{.underline}

*[Risk metric uncertainty]{.underline}*

The 100(1-a)% confidence intervals are often used as characterizations
of uncertainty. These can be calculated for the risk metrics using the
estimated distributions of and . The confidence intervals for are

[. \[\]]{.underline}

[where
]{.underline}*[t]{.underline}~[a]{.underline}~~[,q]{.underline}~*[ is
the critical value of a ]{.underline}*[t]{.underline}*[ distribution at
P = a and ]{.underline}*[q]{.underline}*[ degrees of freedom. The
corresponding significance level, a, for a hypothesis test, such as "Is
]{.underline}*[l]{.underline}*[ \< ]{.underline}*[b]{.underline}*[" is
the a such that:]{.underline}

[. \[\]]{.underline}

[Confidence intervals on and can be calculated by parametric
bootstrapping from the estimated distributions of and : Normal() and
Gamma(shape = (]{.underline}*[n-1]{.underline}*[)/2, scale = ). A large
number of (,) pairs are randomly generated by sampling from these
distributions and the risk metric Y is calculated (Eqn. or ) for each
pair. The range of Y over the (,) bootstrapped pairs, for which both
parameters are within their respective 100(1-a)% confidence intervals,
defines the 100(1-a)% confidence interval for Y. This and other methods
for calculating confidence intervals for diffusion approximation risk
metrics are discussed in Dennis et al. (1991), page
127-128.]{.underline}

An alternate way to present the level of uncertainty is to estimate the
data support for different values of a risk metric. There are both
frequentist and Bayesian approaches for this (see Wade 2001 for a review
geared towards conservation applications). Holmes (2003) presents a
Bayesian approach, which uses posterior probability distributions to
illustrate data support. That method is adapted here for estimating the
level of data support for the metapopulation risk metrics. Let Y be a
risk metric. The probability that Y is greater than some threshold *j*
given the data is

[ \[\]]{.underline}

[where is the likelihood function given , is the likelihood function
given , are the priors on
]{.underline}*[m]{.underline}~[m]{.underline}~*[ and , and the
normalizing constant is ]{.underline}

[ \[\]]{.underline}

[The posterior distribution of Y is over all
]{.underline}*[j]{.underline}*[. Examples of this calculation for
]{.underline}*[l]{.underline}*[ and the probability of being below
thresholds at the end of 25 years are shown in the salmon examples.
Holmes (2003) supplies Splus code for these calculations.]{.underline}

**23.6. A simulated example**

In this example, a collection of 49 local populations in a 7 x 7 grid
was simulated with neighborhood dispersal. Local populations were
specified with variable mean local growth rates; thus, some *m~i~*'s
were much larger than others. The local growth rates in any given year
were slightly correlated between sites. Thus all sites were more likely
than random to have good and bad years together. Dispersal was variable
between 5 and 10% from year to year and was mainly to the four nearest
neighbors (or two and three for corner or edge sites). In specific
terms, **A**(*t*) was specified with *z~i~*(*t*)'s drawn from a normal
distribution with mean = *m~i~* and variance of 0.0625. The *m~i~* were
different for each local population and were chosen randomly between
-0.22 and -0.01. Each year, new *z~i~*(*t*)'s were selected from the
normal distribution for that local population. The *z~i~*(*t*)'s were
correlated among the local populations such that the covariance of
*z~i~*(*t*) and *z~j~*(*t*) was (0.1)(0.0625). The *d~i~*(*t*) varied
from year to year. Each year and separately for each local population,
*d~i~*(*t*) was selected from a uniform random distribution between 0.05
and 0.1; thus the dispersal varied from year to year and between local
populations in any given year. Most of this dispersal, 80%, was to
nearest neighbors. Thus for nearest neighbors, *a~ji~* = 0.80
*d~j~*(*t*)/*nn,* where *nn* is the number of nearest neighbors, and for
non-nearest neighbors, *a~ji~* = 0.2 *d~j~*(*t*)/*nnn; *where *nnn* is
the number of non-nearest neighbors.

The simulation was started from a set of local population sizes randomly
drawn from the stochastic equilibrium, and the starting sizes were drawn
anew from this distribution for each replicate of the simulation. For
each replicate, a 25-year time series was generated, and from this time
series, *m~m~* and were estimated using Eqn. . From the estimates, the
probability that the metapopulation would be below different thresholds
(50% or 75% of starting levels) at the end of 25 years was predicted and
compared to the actual probabilities obtained by repeatedly (1000 times)
running the simulation for 25 years starting from the point where the
initial 25-year time series stopped. This simulation was replicated 500
times to generate the distribution of estimated probabilities of 50 and
75% decline in 25 years versus the true probability. Also, from each
25-year simulation, the metapopulation *l* was estimated and compared to
the actual value calculated by running a 10000-year simulation. For each
estimated risk metric, confidence intervals were estimated via the
methods in section 23.5.

[Figure 23. shows the distribution of ]{.underline}*[l]{.underline}*[
estimates and the estimated probabilities of 50 and 75% decline versus
the true values. As expected, the median estimate of
]{.underline}*[l]{.underline}*[ was equal to the true value ( is an
unbiased estimator of ]{.underline}*[m]{.underline}~[m]{.underline}~*[).
The median estimate of ]{.underline}*[l]{.underline}*[ was 0.97 compared
to the true value of 0.97. The median estimates of 50 and 75% decline
were 0.63 and 0.14 compared to the true values of 0.62 and 0.12,
respectively. Although the median estimates were very close to the true
values, the estimates were variable. The estimates of
]{.underline}*[l]{.underline}*[ ranged between 0.9 and 1.0. The
estimates of declines to thresholds were also variable. The variability
depended on the threshold and the time frame. In this example, there was
low variability around the estimate of 50% decline in 25 years, but high
variability in the estimate of 75% decline. The true values for each of
the metrics are shown by the solid lines in the middle of the
distributions. ]{.underline}

[The variability of the estimates is due to the stochastic nature of the
process and is not a fault of the estimation methods per se; by chance
short trajectories will appear to have underlying parameters that are
different than the true underlying parameters and this leads to
variability in the estimated viability metrics. When estimates are
inherently variable, it is critical that the confidence intervals for
the estimates be correctly estimated. Figure 23. (right panels) confirms
that the estimated confidence intervals properly characterize the
uncertainty for the estimate risk metrics: e.g. 100(1-a)% of the time
the 100(1-a)% confidence intervals contain the true values.]{.underline}

# [23.7. Salmon as metapopulations]{.underline}

[Salmonid populations (Oncorhynchus spp.) show strong spatial
structuring and they have often been referred to as metapopulations
(Reiman and McIntyre 1995; Policansky and Magnuson 1998; Cooper and
Mangel 1999; Hill et al. 2002). Spawning and rearing habitats of
different salmon stocks occur on discrete and physically separated river
or stream sections. Salmon have a well-known and strong tendency to
return to their natal streams with a low (1 to 20%) dispersal to other
stocks (Fulton and Pearson 1981; Mathews and Waples 1991; Quinn 1993).
Within the U.S. Pacific Northwest, collections of anadromous salmon
stocks have been divided into "Evolutionary Significant Units" or ESUs
(Waples 1991) which represent substantially reproductively isolated
conspecific groups that can be distinguished based on their coherence on
a genetic level and known dispersal between the stocks. Salmon within a
stock spawn on individual streams or river sections and the majority of
offspring return to spawn in their natal stream or river. Straying of
returning adults to non-natal streams is spatially structured and occurs
more frequently within subbasins. The stocks within an ESU have some
level of synchrony due to exposure to common migratory corridors between
the ocean and natal stream and also due to exposure to similar
large-scale ocean dynamics (Pearcy 1992; Ware 1995; Mantua et al 1997).
However stocks also show a great deal of asynchrony due to exposure to
their independent spawning and juvenile rearing habitats and variability
in migration timing between stocks (e.g. PSTRT 2001). Throughout the
Pacific Northwest, most salmonid populations show regional decline with
the majority of individual stocks showing steady declines with densities
well below historical levels (Rieman and Dunham 2000; McClure et al.
2003)]{.underline}

#### [23.8. Snake River Spring/Summer Chinook ESU]{.underline}

The Snake River spring/summer chinook ESU (Fig. 23.) includes all spring
and summer chinook spawning within the subbasins of the Tucannon River,
Grande Ronde River and the South, Middle and East Fork Salmon Rivers,
which flow into the Snake River below the Hells Canyon Dam (Mathews and
Waples 1991). Juvenile fish rear in the mountain streams and then
migrate down the Snake and Columbia Rivers to the ocean. After maturing
in the ocean, adult fish return to spawn at variable ages between 3 to 5
years (mean = 4.5 years). Tagging experiments in Columbia River basin
(which the Snake River basin is a part of) have found that the
proportion of individuals that disperse and spawn away from their natal
sites is on the order of 1-3% for wild-born individuals (Quinn 1993).

[The Snake River spring/summer chinook ESU was listed as threatened
under the U.S. Endangered Species Act in 1992. The stocks within this
large and complex basin, like salmon stocks throughout the Pacific
Northwest, are negatively impacted by a variety of factors (Wissman et
al. 1994) and many have experienced substantial declines (Myers et al.
1998; McClure et al. 2003). There is habitat degradation in many areas
related to forestry, grazing, mining and irrigation practices, resulting
in lack of pools, high temperatures, low flows, poor overwintering
conditions, and high sediment loads in many areas. At the same time, a
substantial portion of the ESU is protected as part of federally
designated wilderness (Mathews and Waples 1991). The official ESU
designation does not include salmon in the Clearwater basin since
chinook in this subbasin originate from hatchery fish that were stocked
in the subbasin after the original natural fish were extirpated in the
1940s. However from a metapopulation dynamic perspective, the current
stocks in the Clearwater River basin interact with the stocks within
other subbasins. Thus, in this analysis all stocks in the entire Snake
River basin were analyzed together. ]{.underline}

[A total metapopulation-level time series was available for this ESU
from counts of the total number of wild-born spawners returning through
the Ice Harbor dam at the downstream end of the ESU (Fig. 23.).
Returning spawners can be either wild-born or hatchery-born as
hatcheries have been operating in the basin since the early 1970s.
McClure et al. (2003) discuss the effects of hatchery production on
viability analyses. Effectively by focusing on the wild-born spawner
time series and not incorporating a correction for hatchery production,
the in-stream viability metrics assume that hatchery-born fish all
return to the hatchery and do not spawn in-stream (which would produce
wild-born offspring). As discussed by McClure et al., this means our
viability metrics are optimistic upper bounds, since some unknown
fraction of hatchery fish do stray to the wild spawning grounds and
potentially reproduce. ]{.underline}

[In addition to the metapopulation-level dam count, time series of
redds-per-mile (rpm), which are indices of the density of gravel egg
nests made by spawning females, were available for the majority of
stocks within the Snake River basin. Redds-per-mile are an index of the
redds (and consequently returning spawners) trend within a stock, but
the total redds are unknown since the total spawning habitat is not
surveyed. The majority of the rpm and dam data are available in the
digital appendices of McClure et al. (2003).]{.underline}

#### [Parameter estimation]{.underline}

Our Ice Harbor dam time series starts in 1962 and ends in 1999. The
wild-born component of the dam count is denoted *M*(0), *M*(1), *M*(2),
... *M*(37) where *M*(0) is the 1962 count and *M*(37) is the 1999
count. The maximum likelihood estimates presented in Eqn. assume that
the data do not contain sampling error or other non-process error,
however salmon data typically have high levels of sampling error and
boom-bust cycles that confound estimation of *m~m~* and especially
(Holmes 2001). An alternate approach uses data filtering and examination
of the rate at which variance increases within the time series to
improve parameter estimation and separate out sampling error variance
from the time series (Holmes 2001; Holmes and Fagan 2002; cf also Morris
and Doak 2002 pgs. 171-179). These methods have been extensively
cross-validated with salmon data (Holmes and Fagan 2002; Fagan et al.
2003) and are used here to estimate parameters. First, the data are
transformed using a running-sum:

[ \[\]]{.underline}

[The estimates of ]{.underline}*[m]{.underline}~[m]{.underline}~*[ and
are then]{.underline}

[ \[\]]{.underline}

[The estimate of uses the property that the variance of the underlying
stochastic process should increase linearly with time: . The confidence
intervals for and using are slightly different than Eqn. (Holmes and
Fagan 2002):]{.underline}

[ . \[\]]{.underline}

[where ]{.underline}*[L]{.underline}*[ is the number of counts summed
together for the running sum and ]{.underline}

*[df]{.underline}*[ ≈ 0.333 + 0.212 (]{.underline}*[n]{.underline}*[+1)
-- 0.387 ]{.underline}*[L]{.underline}*[ = 6.84
(]{.underline}*[L]{.underline}*[ = 4 and ]{.underline}*[n]{.underline}*[
= 38 here). The estimated 95% confidence intervals on
]{.underline}*[m]{.underline}~[m]{.underline}~*[ and are (-0.133, 0.020)
and (0.017, 0.111), respectively.]{.underline}

#### [Metapopulation viability metrics]{.underline}

The estimate of *l* for the Snake River spring/summer chinook ESU is .
To the extent that long-term trends continue, the expected population
size in 25 years is 21% of current levels (=). The point estimate of the
probability of that the ESU drops below 10% of current levels at any
time over the next 25 years is

[ \[\]]{.underline}

[The corresponding estimate of 90% decline over the next 50 years is
0.74. The probability of extinction was not estimated since this
requires an estimate of the total population size. The number of
returning spawners is not the total population size since non-mature
fish remain in the ocean. However, if the true
]{.underline}*[l]{.underline}*[ of the metapopulation is less than 1.0,
the population will eventually go extinct.]{.underline}

[ The posterior probability density functions (Eqn. ) for the estimated
metrics are shown in Fig. 23.. The posterior probability distributions
give an indication of the degree to which the data support different
risk levels. The distribution for ]{.underline}*[l]{.underline}[
]{.underline}*[shows considerable data support for a
]{.underline}*[l]{.underline}*[ \< 1, indicating a declining
metapopulation. There is also strong data support for a high risk of 90%
decline over the next 50 years, however the estimate of 90% decline over
25 years is very uncertain. The mean value is 0.23, but the probability
distribution is very broad over the 0 to 1 range. This illustrates that
uncertainty in estimates of probabilities of quasi-extinction can vary
widely depending on the time frame over which one is
interested.]{.underline}

#### [23.9. Puget Sound Chinook ESU]{.underline}

The Puget Sound ESU is a subset of the major chinook salmon group in
Washington's northern coastal basins and Puget Sound. The ESU (Fig. 23.)
includes all spring-, summer-, and fall-runs in the Puget Sound region
from the North Fork Nooksack River to the Elwha River on the Olympic
Peninsula (Myers et al. 1998). The Elwha and Dungeness coastal basins of
the Strait of Juan de Fuca, Hood Canal, and the Puget Sound area north
to the northern Nooksack River Basin and the U.S. Canadian Border are
all a part of the Puget Sound ESU. Basin-to-basin dispersal rates have
been observed at between 0.1 and 6% based on recoveries of tagged
juveniles returning as adults (PSTRT 2001). Fish in this ESU typically
mature at ages 3 and 4 and are coastally-oriented during the ocean phase
of their life history. The Puget Sound ESU does not include Canadian or
coastal Washington populations. The Puget Sound ESU was listed as
threatened under the Endangered Species Act in March of 1999. Trends in
abundance throughout the ESU are predominantly downward with several
populations exhibiting severe short-term declines. Degraded spawning and
rearing habitats as well as access restrictions to spawning grounds and
migration routes have all likely contributed to population declines.
Salmon in this ESU do not migrate through a hydropower system as the
Columbia River ESUs do.

[Our data for this ESU consist of yearly estimates of the total
returning spawners (wild- plus hatchery-born) to the 44 separate river
and creek systems feeding into the Puget Sound (Fig. 23.). These time
series were compiled by the National Marine Fisheries Service (Seattle,
WA) based on a variety of data: redd counts, carcass counts, in-stream
harvest records, weir counts and hatchery return counts. An independent
metapopulation-level count was not available; unlike spawners returning
to the Columbia River basin, spawners here do not pass through a
hydropower system where they can be enumerated. Instead a 1979-1997
index of the metapopulation was constructed by added together the 29
time series for the local populations with data over the 1979-1997
period. As for the Snake River analyses, our viability metrics
implicitly assumes that hatchery fish have not been reproducing and will
be optimistic if some hatchery fish do not return to the hatchery and
instead successfully spawn in the wild.]{.underline}

*[Metapopulation viability metrics]{.underline}*

[Parameters were estimated as for the Snake River. The parameter
estimates are , and . The estimate of ]{.underline}*[l]{.underline}*[
for the Puget Sound Chinook ESU is . The point estimate of the
probability that the ESU drops below 10% of current levels at any time
over the next 25 years is 0.000 and over the next 50 years is 0.001.
]{.underline}

[ The posterior probability distributions (Fig. 23., right panels)
illustrate the high uncertainty, given the data, as to whether this ESU
is declining, stable, or increasing. The most that can be said from the
data is that there is low data support for a severely declining
(]{.underline}*[l]{.underline}*[ \< 0.9) or increasing
(]{.underline}*[l]{.underline}*[ \> 1.1) metapopulation. Interestingly,
the low support for small ]{.underline}*[l]{.underline}*['s translates
into high data support for a low risk of 90% decline in the short-term
(over 25 years). Over the longer term, however, the uncertainty as to
whether the metapopulation is declining or increasing gives rise to a
U-shaped distribution meaning that the data give the most support to a
probability of 0 or 1, reflecting that ]{.underline}*[l]{.underline}*[
could be either less than or greater than 1.0. This example illustrates
that while the data may be equivocal on some questions of conservation
concern, such as "Is ]{.underline}*[l]{.underline}*[ \< 1?", the data
may still give information on other questions, such as "Are the
short-term risk of severe decline high?".]{.underline}

**[23.10. Using the stochastic metapopulation model to investigate
effects of management]{.underline}**

[Determining how to distribute effort in order to recover an at-risk
species is a routine, and challenging, task of conservation managers.
For salmon, management actions tend as a generality to affect an entire
ESU or multiple ESUs or to affect individual stocks. Management actions
such as harvest reductions or increases to survival during migration
(between spawning areas and the ocean) or improvements to estuarine
environments are examples of actions that will tend to improve
conditions for all stocks within an ESU or multiple ESUs. Habitat
improvements or protections that affect spawning areas and management of
in-stream water levels are examples of actions that tend to affect
individual stocks. Without knowing the local stock dynamics or dispersal
rates, one can still give certain types of guidance about how much
effort is required for recovery of a declining metapopulation and about
how effort should be distributed across all local populations.
]{.underline}

#### [Metapopulation-level actions]{.underline}

[When management actions affect all local populations roughly equally,
it can be estimated how change would change the metapopulation
]{.underline}*[l]{.underline}*[. Mathematically, this means that all
]{.underline}*[m]{.underline}~[i]{.underline}~*['s increase by some
d]{.underline}*[m]{.underline}*[. An absolute
d]{.underline}*[m]{.underline}*[ change in all
]{.underline}*[m]{.underline}~[i]{.underline}~*['s is equivalent to
multiplying all elements in
]{.underline}**[A]{.underline}**[(]{.underline}*[t]{.underline}*[) by a
constant = exp(d]{.underline}*[m]{.underline}*[). The mean of the
distribution of log
]{.underline}*[M]{.underline}*[(]{.underline}*[t]{.underline}*[)/]{.underline}*[M]{.underline}*[(0)
becomes . Thus . The change, d]{.underline}*[m]{.underline}*[, can be
translated into currency that is more meaningful from a management
standpoint by using the relationship ]{.underline}*[l]{.underline}*[ =
]{.underline}*[R]{.underline}*~[0]{.underline}~^[1/T]{.underline}^[,
between ]{.underline}*[l]{.underline}*[, the net reproductive rate,
]{.underline}*[R]{.underline}*~[0]{.underline}~[, and the mean
generation time, ]{.underline}*[T]{.underline}*[ (Caswell 2001, Eqn.
5.73). This is illustrated here for harvest and hydropower effects on
salmon in the Snake River spring/summer chinook ESU (cf. McClure et al.
2003).]{.underline}

#### [Harvest: In the Pacific Northwest, harvest rates for salmon are generally expressed in terms of the fraction of spawners that did not return to the spawning grounds but that would have without harvest, e.g. a harvest rate of 0.8 indicates that the actual number of returning spawners is 20% of what it would have been if there had been no harvest. Harvest rates are expressed in this way so that harvest that occurs in the stream versus in the ocean can be compared via a common currency. We can write the net reproductive rate using fecundity and age-specific survival (cf. Caswell 2001, Eqn. 5.61)]{.underline}

[ \[\]]{.underline}

[where ]{.underline}*[h]{.underline}*[ is the harvest rate,
]{.underline}*[s]{.underline}~[i]{.underline}~*[ is the survival from
age ]{.underline}*[i]{.underline}*[-1 to
]{.underline}*[i]{.underline}*[,
]{.underline}*[F]{.underline}~[i]{.underline}~[ ]{.underline}*[is the
fraction of spawners that return at age ]{.underline}*[i]{.underline}*[,
and ]{.underline}*[f]{.underline}*[ is the mean offspring per spawner.
Using Eqn. , the change in the ]{.underline}*[l]{.underline}*[ from a
change in ]{.underline}*[h]{.underline}*[ alone is:]{.underline}

[ \[\]]{.underline}

#### [Hydropower: Juvenile salmon from the Snake River basin must migrate through the mainstem of the Snake River, enter the Columbia River and descend down the Columbia River on their journey to the ocean. This migration, and the return migration of spawning adults, involves passage through four large hydropower dams on the Columbia River and four Snake River hydropower dams. Improving the survival of both juvenile and adult fish migrating through the Columbia and Snake River hydropower systems has been the focus of much effort and is one of the human impacts that has been relatively well-quantified. ]{.underline}

Following a strategy similar to that used for harvest, the effect of
changes in survival through the hydropower system on the rate of decline
at the ESU-level can be estimated. Denoting by *c~d~* and *c~u~* the
proportional increase in down- and upstream passage survival due to
improvement in the hydropower system, the improved net reproductive rate
is

[. \[\]]{.underline}

[Thus, for the assessing the impacts of increased survival through the
hydropower system:]{.underline}

[ \[\]]{.underline}

#### [Estimates of the impacts of harvest and hydropower changes to the Snake River ESU: The mean ocean and in-river 1980-1999 harvest rate for the Snake River spring/summer chinook ESU was ]{.underline}*[h]{.underline}*[ = 0.08 (McClure et al. 2003). By setting ]{.underline}*[h]{.underline}~[new]{.underline}~*[ = 0, we can examine the effect of successful selective harvest management that would substantially eliminate harvest impacts on salmon in this ESU. Using Eqn.  and a mean generation time of 4.5 years, the estimated increase in ]{.underline}*[l]{.underline}*[ with ]{.underline}*[h]{.underline}*~[new]{.underline}~[ = 0 is roughly 2%. NMFS (2000) has recently required that agencies operating the Federal Columbia River Power System implement a variety of activities, including increased spill, improved passage facilities, and increased barging of salmon around the dams as a means of improving survival through the system. The estimated improvement in passage survival from the improvements proposed by NMFS are on the order of 5-6% (i.e. ]{.underline}*[c]{.underline}~[d]{.underline}~[c]{.underline}~[u]{.underline}~*~[\ ]{.underline}~[= 1.05-1.06) for the Snake River spring/summer chinook (McClure et al. 2003). This translates into a 1% improvement in ]{.underline}*[l]{.underline}*[ for this particular ESU, using Eqn. . Thus if the combined effects of substantially reduced harvest and the proposed passage improvements are additive, then roughly a 3% increase in ]{.underline}*[l]{.underline}*[ is estimated for these actions. If the true ]{.underline}*[l]{.underline}*[ is less than 0.97, a 3% increase would not be sufficient to achieve ]{.underline}*[l]{.underline}*[ \> 1. Figure 23. indicates that the data cannot rule out that the ]{.underline}*[l]{.underline}*[ in this is ESU is greater than 0.97, but the data certainly give more support to a lower ]{.underline}*[l]{.underline}*[. This suggests that other recovery actions such as improvements at the stock-level will also be necessary.]{.underline}

#### [Local population-level actions]{.underline}

The effects of changes to individual units of habitat are harder to
quantify than the effects of metapopulaion-level changes. The change in
*l* achieved by a change at the level of a specific unit of habitat
depends on the level of dispersal, the spatial pattern of dispersal,
whether that habitat is connected to source or sink habitat, the level
and pattern of synchrony between sites, etc. In other words, it depends
on the type of detailed information that has traditionally been
difficult to obtain for metapopulations of conservation concern.
Interestingly although it is difficult to determine how much change in
*l* can be achieved, it appears possible to estimate where the largest
d*l* from a given d*m* change (per unit of habitat) in the local growth
rate is achieved, even though the size of the resultant d*l* cannot be
determined.

[Recall that each row of ]{.underline}**[A]{.underline}**[ represents a
unit of habitat, and that a local population is comprised of some set of
units of habitat with high connectivity. When the intrinsic growth rate,
]{.underline}*[m]{.underline}~[j]{.underline}~*[, in a unit of habitat
]{.underline}*[j]{.underline}*[ is changed by
d]{.underline}*[m]{.underline}*[ , to
exp(]{.underline}*[m]{.underline}~[j]{.underline}~*[
+d]{.underline}*[m]{.underline}*[), all the
]{.underline}*[a]{.underline}~[ij]{.underline}~*[ elements of column
]{.underline}*[j]{.underline}*[ in matrix
]{.underline}**[A]{.underline}**[(]{.underline}*[t]{.underline}*[) are
multiplied by exp(d]{.underline}*[m]{.underline}*[). The goal is to
calculate the total change in ]{.underline}*[l]{.underline}*[ from this
d]{.underline}*[m]{.underline}*[ change to all elements in column
]{.underline}*[j]{.underline}*[ by summing over rows
]{.underline}*[i]{.underline}*[:]{.underline}

[ \[\]]{.underline}

[The term is the elasticity of ]{.underline}*[l]{.underline}*[. Caswell
(2001), section 14.4, presents the calculation for the elasticity of
]{.underline}*[l]{.underline}*[ for products of stochastic
matrices:]{.underline}

[ \[\]]{.underline}

[where ]{.underline}*[R]{.underline}*[(]{.underline}*[t]{.underline}*[)
is the relationship between the right eigenvector and
]{.underline}**[A]{.underline}**[(]{.underline}*[t]{.underline}*[), .
Thus, the d]{.underline}*[l]{.underline}*[ from a
d]{.underline}*[m]{.underline}*[ change in a unit of habitat
]{.underline}*[j]{.underline}*[ can be solved for by summing Eqn. over
]{.underline}*[i.]{.underline}*

[ \[\]]{.underline}

[This can be translated into matrix form: ]{.underline}

[ \[\]]{.underline}

[where "°" denotes the scalar product. Using the relationship between
the left eigenvector and
]{.underline}**[A]{.underline}**[(]{.underline}*[t]{.underline}*[),
,]{.underline}

[ \[\]]{.underline}

[Thus d]{.underline}*[l]{.underline}*[ from a change in a unit of
habitat ]{.underline}*[j]{.underline}*[ is a weighted temporal average
of the reproductive value of local population
]{.underline}*[j]{.underline}*[ times it's density:]{.underline}

[ \[\]]{.underline}

[where ]{.underline}*[c]{.underline}*[(]{.underline}*[t]{.underline}*[)
is a constant that depends on ]{.underline}*[t]{.underline}*[ but not
]{.underline}*[j]{.underline}*[. A local population
]{.underline}*[a]{.underline}*[ is comprised of units of habitat in the
set ]{.underline}*[a]{.underline}*[ =
{]{.underline}*[a]{.underline}*~[1]{.underline}~[,
]{.underline}*[a]{.underline}*~[2]{.underline}~[,
]{.underline}*[a]{.underline}*~[3]{.underline}~[, ...,]{.underline}*[
a]{.underline}~[m]{.underline}~*[}, where
{]{.underline}*[a]{.underline}*~[1]{.underline}~[,
]{.underline}*[a]{.underline}*~[2]{.underline}~[,
]{.underline}*[a]{.underline}*~[3]{.underline}~[, ...,]{.underline}*[
a]{.underline}~[m]{.underline}~*[} denotes which rows of
]{.underline}**[A]{.underline}**[ correspond to the units of habitat in
local population ]{.underline}*[a]{.underline}*[. The
d]{.underline}*[l]{.underline}*[ per d]{.underline}*[m]{.underline}*[
per unit habitat for a particular local population
]{.underline}*[a]{.underline}*[ is , where
]{.underline}*[m]{.underline}*[ is the number of units of habitat in
local population ]{.underline}*[a]{.underline}*[. In words, this means
the change in ]{.underline}*[l]{.underline}*[ is proportional to the
product of the 'average' density of individuals in a particular local
population times the 'average' reproductive value of its units of
habitat. ]{.underline}

Although reproductive values are unknown, there are many cases where the
product *v~j~w~j~* is a positive function of *v~j~* as long as dispersal
is not too unidirectional (meaning, dispersal from A to B but not B to
A). This can be shown analytically in three extreme cases: a) 100%
uniform and equal dispersal, b) all *m~j~*'s equal, or c) dispersal
extremely low. In cases a) and b), the reproductive values are all equal
and *v~j~w~j~* = (a constant) x *v~j~*. In case c), and . However, this
positive relationship was also found in simulations with variable local
growth rates, neighborhood dispersal, and dispersal sources and targets.
An obvious exception to this positive relationship is if dispersal is
unidirectional, for example a linear chain of local populations with
dispersal via a steady directional wind or ocean current. However, as
the simulations below illustrate, the general relationship can still
hold even when dispersal is strongly, although not strictly,
directional.

*[Density and ]{.underline}[l]{.underline}[ sensitivity]{.underline}*

[Three different types of metapopulation models were used to look at the
relationship between average densities in units of habitat versus the
d]{.underline}*[l]{.underline}*[ from a small increase in the local
growth rate in each unit of habitat. In each model, dispersal was
non-uniform among the local populations so that some sites were
dispersal sources (more dispersal out than in) and others dispersal
targets (more dispersal in than out). In the first model, local growth
rates were equal among all sites and dispersers were globally
distributed among all sites. In the second model, local growth rates
were variable so that some sites had much higher local growth rates than
others and dispersal was mainly to nearest neighbors. In the third
model, local growth rates were again variable and dispersal mainly to
the two south and east neighbors, however a small proportion of
dispersers were distributed globally. Thus the three examples illustrate
global, local and directional dispersal.]{.underline}

A hundred randomly generated matrix models in each of these three
categories were made and d*l~j~* calculated via Eqn. . Figure 23. shows
the relationship between the average density of a local population and
the d*l* from increasing the local growth rate in that unit of habitat.
The *x*-axis ranks the d*l~j~*, thus "1" indicates the local population
with the highest d*l~j~* in any simulation and "49", the lowest. The
*y*-axis shows the corresponding mean density rank of that local
population; "1" indicates the population had the highest density amongst
the 49 sites and "49", the lowest. The results from the 100 randomly
generated models are summarized by showing a box plot, which shows the
median and range of all density ranks for the sites with a given
d*l~j~*/d*m* rank. Thus, the box plot at the x-axis position "1" shows
the range of density ranks for the units of habitat with the highest
d*l*/d*m* in each model. The model results show a strong positive
relationship between the relative density rank within a unit of habitat
and which unit of habitat produced the largest increases in the
metapopulation *l* for a given d*m*. The two to three units of habitat
with the highest average densities were consistently the units that
produced the largest d*l* for a given d*m*. This suggests that plotting
the distribution of the relative densities within local populations in a
metapopulation could give a rapid indication of the sensitivity of the
metapopulation to changes to individual local populations.

One application of this would be to estimate where negative impacts
would lead to the greatest decrease in *l *, thus suggesting where
protection in most critical. It would also suggest where local
improvements would be most effective for a given increase in the local
growth rate. However in actual management situations where improvements
are being sited, one is generally trying to maximize the "bang per
buck", . The cost, d\$, is actual monetary cost or some combination of
monetary, logistical and political costs and d*m~j~*/d\$ is the cost of
a unit improvement to a unit of habitat *j*. Thus d*l~j~*/d*m* is one
part of the equation, and the other part, the cost of a unit improvement
in different habitats, would have to come from a specific analysis of
the costs and estimated effects of management actions on different local
populations.

#### [Example using the Snake River ESU: The overall level of salmon dispersal between and among stocks within this ESU is known to be fairly low and spatially localized (Mathews and Waples 1991; Quinn 1993). In addition, there is high variability in the habitat quality between stocks, with some stocks relatively pristine and protected within wilderness areas while others are exposed to high and multiple impacts (such as stream degradation and disturbance, pollution, in-stream harvest, and irrigation impacts). Figure 23. (top) shows the distribution of average normed redds-per-mile per mile for 50 Snake River spring/summer chinook stocks. For each year between 1980-1995, the redds-per-mile count for each stock was divided by the maximum count among the 50 stocks in that year. The average over the 16 years was then used as an estimate of the average normed redds-per-mile. The long-tailed distribution is the expectation from theory given low dispersal and high variability in stock habitat quality. ]{.underline}

Estimation of the average normed redds-per-mile was repeated using a
variety of different time periods. Regardless of the time period or
number of years used for averaging, six stocks consistently appeared
among the top five stocks with the highest density of redds: Johnson
Cr., Poverty Cr. and Secesh R. in the South Fork of the Salmon R., the
Lostine R. in the Grande Ronde subbasin, Marsh Cr. in the Middle Fork of
the Salmon R., and the Imnaha R. Perhaps not surprisingly, all of these
are in relatively low impacted regions of the ESU. At a subbasin level,
the overall highest redd density was in the South Fork of the Salmon R.
where summer-run chinook primarily occur (Fig. 23., bottom). The other
regions are primarily spring-run chinook. The South Fork of the Salmon
R. is relatively pristine and few hatchery fish have been released into
this subbasin, and the stocks presumably have experienced relatively low
interbreeding with hatchery-reared stocks. In addition, the later
run-timing may somehow be associated with less straying, lower harvest
or lower hydropower impacts.

This analysis predicts that the *l* of the Snake R. spring/summer
chinook ESU would be most sensitive to changes to the summer-run stocks
in the South Fork of the Salmon River and to the spring-run stocks, the
Lostine R., Imnaha River, and Marsh Creek and should be preferentially
protected from impacts. This can be counterintuitive in some situations.
For example, imagine making choices about where to allow a limited
catch-and-release sport fishery. Sites with the highest density would
seem to be the prime candidates whereas the analysis of d*l*/d*m*
indicates just the opposite. In terms of determining where to direct
improvements, the d*l*/d*m* suggests that these pristine sites are where
a given d*m* would produce the greatest metapopulation *l*, however the
regions where d*l*/d*m* is the highest are not necessarily the regions
where *m* is most easily improved. Indeed a given unit of improvement
may be more difficult in pristine sites. Choosing where to direct stock
improvements requires consideration of the cost and difficulty of a
given d*m* for different stocks in combination with the estimate of the
sensitivity of *l* to local changes.

# [23.11. PVA in practice]{.underline}

The purpose of this chapter is to present a theoretical framework for
metapopulation PVA using time series data and diffusion approximations.
These methods are then illustrated using data from two salmon
metapopulations. The salmon analyses are intended as an example of how
to calculate the diffusion parameters and metrics. An actual PVA must
grapple with other important issues that are outside the scope of this
chapter, but which anyone contemplating an actual PVA must be aware.
Morris and Doak (2002, chapter 12) give a review of the criticisms and
caveats surrounding the use of PVA and outline general recommendations
and cautions when conducting a PVA. In the context of diffusion
approximation methods in particular, Holmes (2003) outlines an approach
using matrix models to conduct sensitivity analyses in order to chose
among different parameterization methods and metrics for a specific PVA
application.

[One of the issues that is especially pertinent for our chapter is the
issue of variability in estimated risk metrics. A number of recent PVA
cross-validations using actual data on a large number of different
populations have shown that careful PVA analyses give unbiased risk
estimates (Brook et al. 2000; Holmes and Fagan 2002; Fagan et al. 2003).
Although this is very encouraging, a difficult issue is the high
inherent variability associated with estimated probabilities (such as
probability of extinction) -- even though they may be unbiased (Ludwig
1996, 1999; Fieberg and Ellner 2000; Holmes 2001; Ellner et al. 2002).
How to properly use risk metrics that have high variability is currently
being debated within the field with arguments ranging from "don't use
them" (Ludwig 1996, 1999; Fieberg & Ellner 2000), to "use to estimate
risks within collections of populations" (Fagan et al. 2001; Holmes and
Fagan 2002), to "use where data are extensive and high quality" (Coulson
et al. 2001), to "PVA metrics based on data, even if variable, are
better than the alternatives" (Brook et al. 2002). An encouraging aspect
of diffusion approximation methods is that cross-validations using real
time series data have indicated that the uncertainty in the estimated
metrics appears to be properly characterized (Holmes and Fagan 2002).
Nonetheless, how to use and present metrics with high variability,
albeit well characterized, is not an easy question to answer.
Presentation of 100(1-a)% is an oft used approach, but experience in the
forum of salmon recovery planning in the Pacific Northwest has shown
that it is easy to misinterpret confidence intervals. For example, it is
easy to interpret 95% confidence intervals for
]{.underline}*[l]{.underline}*[ that overlap 1.0 as an indication that
the data are equivocal as to whether the population is declining or
increasing; whereas, there may be considerable data support for a
declining population. Graphical presentations of data support for
different risk levels have been more compelling and informative,
although translating levels of data support into numbers that
policymakers can use to take uncertainty into account in policy
decisions has been challenging.]{.underline}

# [23.12. Discussion]{.underline}

This chapter has focused on the calculation of metapopulation PVA
metrics, however there are other more general PVA insights from an
examination of stochastic metapopulations and of this specific class of
declining density-independent metapopulations. First, by definition the
trajectory of a stochastic metapopulation is subject to random processes
and thus the metapopulation trajectory observed in any one snippet of
time is unlikely to capture the long-term dynamics. The shorter the time
frame, the farther the observed trend is likely to be from the long-term
trend. Thus the trends in any two adjacent time periods are unlikely to
be identical, and the difference indicates not necessarily a change in
the underlying rate of decline, can be due simply to chance. The
variability of observed rates of decline can be estimated from the level
of the variability driving the long-term dynamics, and thus statistical
tests performed to determine the likelihood that an apparent change in
trend occurred due to the stochastic nature of the process rather than
an underlying change in conditions.

[Second, the local populations within a metapopulation are linked and
experience the same long-term growth rates, regardless of the underlying
difference in local population conditions (i.e. whether they are
'sources' or 'sinks'). However, measured over a short time period, there
will be differences in the observed local population trends due to
chance and local conditions. This means that over a given time period,
local populations will appear to be declining at different rates, but
this is not an indication the long-term trends and not necessarily
related local conditions being better or worse than other areas. That
the long-term trends of the individual local populations are the same as
the metapopulation has a direct impact on PVA for local populations
within a metapopulation. The rate of decline observed among the
different local populations will differ as will the apparent level of
variability in the local time series. Thus if an individual viability
analysis is done using parameters estimated from local population time
series alone, it will appear that there is tremendous variability among
the local populations risk levels when in fact their long-term risks are
similar. When looking at the long-term risks, use of
metapopulation-level parameters leads to better estimates of the
long-term local population risks. Short-term risks, however, are still
strongly influenced by local conditions. Clearly estimates of both
short-term and long-term risks are needed to capture the whole viability
picture for a metapopulation. Although local populations within the type
of metapopulations modeled here will be eventually repopulated by
dispersal if they undergo extreme declines, the resulting loss of
genetic diversity leads to a gradual erosion of the genetic health of
the metapopulation. Indeed this has happened for salmon species
throughout the Pacific Northwest. ]{.underline}

[Recovery planning for endangered and threatened species typically
requires determining where to put the most effort. Rarely is it the case
that maximal effort can be applied everywhere. Using the stochastic
metapopulation model, a sensitivity analysis was used to look for local
characteristics that predict where local changes would produce the
biggest change in the metapopulation growth rate. Interestingly, local
density (not absolute numbers) was a strong predictor of the where a
unit change in local growth rates led to the largest metapopulation
growth rate. This relationship was observed even in simulations with
dispersal sources and targets and strongly directional dispersal,
although it will break down when dispersal is strictly unidirectional.
Determining which local populations are best suited for restoration
efforts also requires assessing the feasibility, cost, and acceptance of
restoration efforts. Indeed when it comes to actually implementing
recovery actions, optimizing the efficiency of effort in terms affecting
recovery requires solving a complex function of biological, economic and
political information. However, understanding the population dynamics of
the species of concern and gaining insight regarding how the demography
of the species will respond to alternative management actions are
fundamental and primary components of this conservation
equation.]{.underline}

# [23.13. Acknowledgements]{.underline}

The graduate students in the Winter 2000 seminar on Spatial Ecology and
Conservation Biology at the University of Washington were the initial
impetus of this work. Special thanks are given to Amanda Stanley who did
many of the simulations that eventually led to the realization of strong
patterns within the dynamics of stochastic metapopulations. Jon Hoekstra
(NMFS) provided us with data on the Puget Sound chinook ESU. Discussions
with Michelle McClure, John McElhany, and Mary Ruckleshaus (leads of the
NMFS Technical Recovery Teams for Pacific NW salmon) have been
instrumental at various stages of this work. Thanks to Michael McCarthy
whose detailed review and critiques led to many improvements in the
presentation of the theory.

# [Literature Cited]{.underline}

Akçakaya, H. R. 1997. RAMAS ® Metapop: viability analysis for
stage-structured metapopulations. Version 2. Applied Biomathematics,
Setauket, New York, USA.

Beissinger, S.R. 2002. Population viability analysis: past, present,
future. In S.R. Beissinger and D.R. McCullough. Editors. Population
viability analysis. University of Chicago Press, Chicago.

Beissinger, S. R., and M. I. Westphal. 1998. On the use of demographic
models of population viability in endangered species management. Journal
of Wildlife Management 62:821-841

Boyce, M. S. 1992. Population viability analysis. Annual Review of
Ecology and Systematics 23:481-506.

Boyce, M. 1996. Review of RAMAS/GIS. Quarterly Review of Biology
71:167-168.

[Brook, B.W., Burgman, M.A., Akcakaya, H.R., O\'Grady, J.J., and
Frankham, R. 2002. Critiques of PVA ask the wrong questions: Throwing
the heuristic baby out with the numerical bath water. Conservation
Biology 16: 262-263.]{.underline}

[Brook, B. W., J. J. O\'Grady, A. P. Chapman, M. A. Burgman, H. R.
Akçakaya, and R. Frankham. 2000. Predictive accuracy of population
viability analysis in conservation biology. Nature.
404:385-387.]{.underline}

[Caswell, H. 2001. Matrix population models: construction, analysis, and
interpretation. Sinauer Press, Sunderland, MA, USA.]{.underline}

[Cooper, A. B. and M. Mangel. 1999. The dangers of ignoring
metapopulation structure for the conservation of salmonids. Fisheries
Bulletin 97:213--226.]{.underline}

[Coulson, J.C., G. M. Mace, E. Hudson, and H. Possingham. 2001. The use
and abuse of population viability analysis. Trends in Ecology and
Evolution. 16: 219-221.]{.underline}

[Dennis, B., P.L. Munholland and J.M. Scott. 1991. Estimation of growth
and extinction parameters for endangered species. Ecological Monographs
64:205-224.]{.underline}

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

[Fagan, W. F., E. E. Holmes, J. J. Rango, A. Folarin, J. A. Sorensen, J.
E. Lippe, and N. E. McIntyre. 2003. Cross-validation of quasi-extinction
risks from real time series: an examination of diffusion approximation
methods. Submitted; 2]{.underline}^[nd]{.underline}^[ revision in
review.]{.underline}

Fagan, W.F., E. Meir, J. Prendergast, A. Folarin, and P. M. Kareiva.
2001. Characterizing vulnerability to extinction for 758 species.
Ecology Letters 4:132-138.

[Fieberg, J., and S.P. Ellner. 2000. When is it meaningful to estimate
an extinction probability? Ecology 81: 2040-2047.]{.underline}

[Fulton, L. A., and R. E. Pearson. 1981. Transplantation and homing
experiments on salmon, Oncorhynchus spp., and steelhead trout, Salmo
gairdneri, in the Columbia River system: Fish of the 1939-44 broods.
U.S. Dep. Commer., NOAA Tech. Memo. NMFS F/NWC-12, 97p.]{.underline}

[Furstenberg, H. and Kesten, H. 1960. Products of random matrices.
Annals of Mathematical Statistics 31:457-469.]{.underline}

[Gilpin, M.E. and M.E. Soule. 1986. Minimum viable populations of
species extinction. Pp. 19-34 In M.E. Soule (ed.). Conservation biology:
the science of scarcity and diversity. Sinauer Associates, Inc.,
Sunderland, Massachusetts.]{.underline}

[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.]{.underline}

[Hill, M. F., A. Hastings, and L. W. Botsford. 2002. The effects of
small dispersal rates on extinction times in structured metapopulation
models. American Naturalist 160:389-402.]{.underline}

[Holmes, E. E. 2001. Estimating risks in declining populations with poor
data. Proceedings of the National Academy of Science
98:5072-5077.]{.underline}

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

[Holmes, E. E. 2003. Beyond theory to application and evaluation:
diffusion approximations for population viability analysis. Manuscript
conditionally accepted, Ecological Applications.]{.underline}

IUCN. 1994. IUCN Red List Categories, Prepared by the IUCN Species
Survival Commission As approved by the 40th Meeting of the IUCN Council,
Gland, Switzerland.

[Karlin, S. and H. M. Taylor. 1981. A second course in stochastic
processes. Academic Press, New York, NY, USA.]{.underline}

[Lacy, R.C. 1993. VORTEX: A computer simulation model for Population
Viability Analysis. Wildlife Research 20:45-65.]{.underline}

[Lande, R. 1993. Risks of population extinction from demographic and
environmental stochasticity and random catastrophes. American Naturalist
142:911--927.]{.underline}

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

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.

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

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

Mantua, N. J., S. R. Hare, Y. Zhang, J. M. Wallace, and R. C. Francis.
1997. A Pacific interdecadal climate oscillation with impacts on salmon
production. Bulletin of the American Meteorological Society
78:1069-1079.

Mathews, G. M. and R. S. Waples. 1991. Status review of Snake River
spring and summer chinook salmon. U.S. Dept. Commer., NOAA Tech. Memo.
NMFS F/NWC-200.

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

Morris, W.F., P. L. Bloch, B. R. Hudgens, L. C. Moyle, and J. R.
Stinchcombe. 2002. Population viability analysis in endangered species
recovery plans: past use and future improvements. Ecological
Applications 12:708-712.

[Morris, W. F., and D. F. Doak. 2002. Quantitative conservation biology:
theory and practice of population viability analysis. Sinauer
Associates, Inc, Sunderland MA.]{.underline}

Myers, J.M., R.G. Kope, G.J. Bryant, D. Teel, L.J. Lierheimer, T.C.
Wainwright, W.S. Grand, F.W. Waknitz, K. Neely, S.T. Lindley, and R.S.
Waples. 1998. Status review of chinook salmon from Washington, Idaho,
Oregon, and California. U.S. Dept. Commer., NOAA Tech. Memo.
NMFS-NWFSC-35, 443pp.

NMFS (National Marine Fisheries Service). 2000. Operation of the federal
Columbia River power system including the juvenile fish transportation
program and the Bureau of Reclamation\'s 31 projects, including the
entire Columbia Basin Project. U.S. Department of Commerce, NOAA,
National Marine Fisheries Service, Northwest Region, NW Region, Seattle.

Pearcy, W.G. 1992. Ocean Ecology of North Pacific Salmon. University of
Washington Press, Seattle, WA. 179p.

Policansky, D. and J. J. Magnuson. 1998. Genetics, metapopulations, and
ecosystem management of fisheries. Ecological Applications 8:S119-S123.

Possingham, H. P., and I. Davies. 1995. ALEX: a model for the viability
of spatially structured populations. Biological Conservation 73:143-150.

Puget Sound Technical Recovery Team (PSTRT). 2001. Independent
populations of chinook salmon in Puget Sound. Public review draft.
April, 2001. Northwest Fisheries Science Center, NOAA, Seattle, WA, USA.

Quinn, T. P. 1993. A review of homing and straying of wild and
hatchery-produced salmon. Fisheries Research 18:29-44.

Ralls, K., S. R. Beissinger, and J. F. Cochrane. 2002. Guidelines for
using population viability analysis in endangered-species management.
Pages 521-550 in Population Viability Analysis (S. R. Beissinger and D.
R. McCullough, Eds.). University of Chicago Press, Chicago, IL.

Reed, J. M., Mills, L. S.; Dunning Jr, J. B., Menges, E. S., McKelvey,
K. S., Frye, R., Beissinger, S. R., Anstett, M. C., Miller, P. 2002.
Emerging issues in population viability analysis. Conservation Biology
16:7-19.

Rieman, B. E. and J. B. Dunham. 2000. Metapopulations and salmonids: a
synthesis of life history patterns and empirical observations. Ecology
of Freshwater Fish 9:51-64.

Rieman, B. E. and J. D. McIntyre. 1995. Occurrence of bull trout in
naturally fragmented habitat patches of varied size. Transactions of the
American fisheries Society 124:285-296.

Sanz, L. and R. Bravo de la Parra. 1998. Variables aggregation in time
varying discrete systems. Acta Biotheoretica 46:273-297.

Shimizu, K. and K. Iwase. 1981. Uniformly minimum variance unbiased
estimation in lognormal and related distributions. Communications in
Statistics, Part A -- Theory and Methods 10:1127-1147.

Shaffer, M.L. 1981. Minimum population sizes for species conservation.
BioScience 31:131-134.

Tuljapurkar, S. D. 1990. Population dynamics in variable environments.
Lecture Notes in Biomathematics. Springer-Verlag, New York, NY, USA.

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

[Wade, P. R. 2001. The conservation of exploited species in an uncertain
world: novel methods and the failure of traditional techniques. In
Conservation of exploited species, (J. D. Reynolds and G. Mace, eds.),
Cambridge University Press.]{.underline}

Waples, R. S. 1991. Definition of \"species\" under the Endangered
Species Act: Application to Pacific salmon. U.S. Department of Commerce,
NOAA Tech. Memo. NMFS F/NWC-194, 29 pp.

Ware, D. M. 1995. A century and a half of change in the climate of the
NE Pacific. Fisheries Oceanography 4:267-277.

Wilcove, D. S., D. Rothstein, J. Dubow, A. Phillips and E. Losos. 1998.
Quantifying threats to imperiled species in the United States.
Bioscience 48:607-615.

Wissmar, R. C., J. E. Smith, B. A. McIntosh, H. W. Li, G. H. Reeves and
J. R. Sedell. 1994. A history of use and disturbance in riverine basins
of eastern Oregon and Washington (Early 1800s - 1900s). Northwest
Science 68:1-35.

**[Figures]{.underline}**

[Figure [23.]{#figBox1}. Illustration of two of the main results from
random theory: (A) All local population go towards the same long-term
rate of population growth (or decline) as
]{.underline}*[t]{.underline}*[ gets large. (B) The distribution of log
]{.underline}*[M]{.underline}*[(]{.underline}*[t]{.underline}*[) is a
normal distribution with mean given by the long-term rate of growth (or
decline) multiplied by ]{.underline}*[t]{.underline}*[ and the variance
given by ]{.underline}*[t]{.underline}*[ multiplied by the rate that
variance increases in an individual trajectory, i.e. for
]{.underline}*[t]{.underline}*[ not overly small. Here the variance was
estimated from one time series using ]{.underline}*[t]{.underline}*[ =
10 and this used to predict the distribution at
]{.underline}*[t]{.underline}*[ = 100;]{.underline}

Figure [23.]{#figDiffTest}. Illustration of the performance of a
diffusion approximation for modeling the behavior of a metapopulation
with 50 local populations and uniform 0.1%, 1%, or 5% yearly dispersal.
The diffusion approximation performs well for a given time frame when
and are constants over that time frame and when log *M*(*t*)*/M*(0) is
normal.

[Figure [23.]{#figEstMetricPerf}. Estimated viability metrics and their
estimated confidence intervals versus the true values for a 49-site
metapopulation in a 7 x 7 grid with 5-10% dispersal to the closest
neighboring sites. The left panels compared the true metrics to the
distribution of estimated metrics from 500 simulations starting from the
same initial conditions. The right panels show the performance of the
estimated confidence intervals by looking at the fraction of estimated
100(1-a)% confidence intervals that contain the true
values.]{.underline}

Figure [23.]{#figSnakeMap}. Map of the Snake River spring/summer chinook
ESU. The ESU includes the stocks from the Snake River and its
tributaries between the Ice Harbor and Hells Canyon dams. The Hells
Canyon hydropower dam has no passage facilities and blocks the migration
of salmon into their historical habitat in the upper Snake River basin.

Figure [23.]{#figESURiskMetrics}. Estimated posterior probability
distributions for *l*, and the probability of 90% decline in 25 and 50
years. The posterior probability distributions, which were calculated
using uniform priors on and , indicate the relative levels of data
support for different risk metric values. The left and right panels show
the distributions for the Snake River Spring/Summer chinook and Puget
Sound ESUs, respectively.

Figure [23.]{#figPugetSound}. Map of the Puget Sound chinook ESU.

Figure [23.]{#figBangforBuck}. Relationship between the influence of a
given habitat unit on the metapopulation *l* versus the average density
in that habitat unit. One hundred 7x7 metapopulations with spatially
variable dispersal rates (some sites dispersal sinks and others targets)
were randomly generated in each of three classes: 1) spatially uniform
growth rates and global dispersal, 2) spatially variable growth rates
with neighborhood dispersal and 3) spatially variable growth rates with
directional neighborhood dispersal to the S and E two neighbors only.
The x-axis shows the rank in terms of d*l*/d*m* and the y-axis shows a
box plot of the distribution of density ranks for sites with a given
d*l*/d*m* rank across all 100 models in each class. Thus the box plot at
x = 1 shows the distribution of ranks for the sites with the highest
d*l*/d*m* in each model. The line in each box shows the median density
rank for the sites with a given d*l*/d*m* rank, the box encloses 50% of
the ranks, and the whiskers show the range from all 100 randomly
generated models.

Figure [23.]{#figESUdists}. Distribution of densities of redds in the
Snake River spring/summer chinook ESU at a stock and subbasin level. The
average normed redds densities (upper plot) are shown for the 50 stocks
with 1980-1995 data (the years were chosen to maximum the number of
stocks with data). For each stock the normed redd density was averaged
over the 16 years to get an estimate of the normed average density. In
the lower plot, relative average densities over all stocks within
different basins are shown (with the number of stocks in each basin
shown above the bars). The basin designations are: GR = Grande Ronde, I
= Imnaha, SFS = South Fork Salmon, MFS = Middle Fork Salmon, US = Upper
Salmon, C = Clearwater. The redds due to hatchery fish released into
stocks were removed before doing these analyses; otherwise the density
will be artificially high simply due to hatchery fish releases. This
correction could not be done for the Upper Salmon or Clearwater regions
since the fraction of spawners that are hatchery strays were unknown;
however, the hatchery releases are very high in these basins and thus
the corrected relative densities would be much lower than shown.

**Box 23.1 Key results from random theory**

Here are reviewed two of the fundamental results from the theory of
products of random matrices and interpreted in the context of our
metapopulation model. The reader is referred to Caswell (2001) section
14.3.3 and Tuljapurkar (1990) section 2.1 for other reviews interpreted
in the context of demographic, single population models.

*The metapopulation and local populations decline at the same rate*

One of the basic results from Furstenberg and Kesten's "Products of
random matrices" (1960), is that the product of ergodic random matrices
asymptotically goes to an equilibrium. Say that **X**~t~ is an ergodic
random 'k x k' matrix and that **Y** (also a *k* x *k* matrix) denotes
the product of *n* of the **X** matrices: **X**~1~,** X**~2~,**
X**~3~,...** X**~n~. Then Furstenberg and Kersten's results say that,
**Y** goes an equilibrium state such that:

[= a constant which is the same for all ]{.underline}*[a]{.underline}*[
\[1\]]{.underline}

[We can use this result to show that the long-run exponential growth
rate of the metapopulation and the local populations will be the same.
]{.underline}

[Thus from Eqn. 1,]{.underline}

[== a constant = ]{.underline}*[m]{.underline}~[m]{.underline}~*

*[The distribution of local population and metapopulation sizes is
lognormally distributed]{.underline}*

[One of the most powerful results, for our purposes at least, concerns
the statistical distribution of the metapopulation and local
trajectories. This tells us what distribution of sizes we would see if
we ran our model over and over again, and allows us to make population
viability analyses for metapopulations since we have a prediction about
the likelihood of different metapopulation futures. Random theory
(Furstenberg and Kersten 1960, Tuljapurkar and Orzack 1980) shows that
any sum of the
]{.underline}**[N]{.underline}***~[i]{.underline}~*[(]{.underline}*[t]{.underline}*[)'s,
such as the total metapopulation (all
]{.underline}*[i]{.underline}*['s), a single local population (one
]{.underline}*[i]{.underline}*[), or any other subset, goes to the same
distribution:]{.underline}

[ \[2\]]{.underline}

[where the sum of local populations is denoted in matrix terms as and
]{.underline}**[c]{.underline}**[ is a column vector with 0's and 1's to
show which local populations to sum together.]{.underline}

# [Example]{.underline}

These results are simple to see with simulations. Here is shown an
example of a linear chain of ten local populations connected via 2%
yearly dispersal to their nearest neighbors and 0.2% to non-nearest
neighbors. The local dynamics were where z*~i~* is a normally
distributed random variable, Normal(*m~i~,*). The local growth rates,
*m~i~*'s, for local populations 1 to 10 were respectively: 0.97, 1.00,
0.96, 0.83, 0.88, 1.00, 1.00, 0.89, 0.99, 0.81. Figure 1A shows that the
long-run growth rate of the local population and metapopulations is
equal to the same constant. Figure 1B shows that the distribution of
metapopulation size after 100 years is . The expected distribution was
estimated using the maximum likelihood (ML) estimates for *m~m~* and
(Eqn. ) from a single 1000-year time series of metapopulation counts.
The ML estimate for relies on an assumption of normality for *t* = 1,
although strictly speaking normality only holds for *t* large. However,
it does quite well as can be seen in Figure 1B.