Incorporating Multiple Sources of Stochasticity into Dynamic Population Models


Calder, C.; Lavine, M.; Muller, P.

Ecology 84(6): 95-402

2003


The application of standard Bayesian techniques to time-series models of population dynamics is discussed. Many standard statistical models used to assess population dynamics ignore significant sources of stochasticity. Standard time-series models for population dynamics can be extended to include both observational and process error. They can also be used to perform inference on parameters in a Bayesian setting. Analysis of simulated data demonstrates that ignoring observation error can be misleading.

epEcology,
84(6),
2003,
pp.
1395-1402
C)
2003
by
the
Ecological
Society
of
America
INCORPORATING
MULTIPLE
SOURCES
OF
STOCHASTICITY
INTO
DYNAMIC
POPULATION
MODELS
CATHERINE
CALDER,
1,4
MICHAEL
LAVINE,
1
PETER
MCILLER,
1,2
AND
JAMES
S.
CLARK
S
linstitute
of
Statistics
and
Decision
Sciences,
Duke
University,
Box
90251,
Durham,
North
Carolina
27708
USA
2
Department
of
Biostatistics,
University
of
Texas
M
D.
Anderson
Cancer
Center,
Houston,
Texas
77030-4009
USA
3
Department
of
Biology
and
Nicholas
School
of
the
Environment,
Duke
University,
Durham,
North
Carolina
27708
USA
Abstract.
Many
standard
statistical
models
used
to
examine
population
dynamics
ig-
nore
significant
sources
of
stochasticity.
Usually
only
process
error
is
included,
and
un-
certainty
due
to
errors
in
data
collection
is
omitted
or
not
directly
specified
in
the
model.
We
show
how
standard
time-series
models
for
population
dynamics
can
be
extended
to
include
both
observational
and
process
error
and
how
to
perform
inference
on
parameters
in
these
models
in
the
Bayesian
setting.
Using
simulated
data,
we
show
how
ignoring
observation
error
can
be
misleading.
We
argue
that
the
standard
Bayesian
techniques
used
to
perform
inference,
including
freely
available
software,
are
generally
applicable
to
a
variety
of
time-series
models.
Key
words:
Bayesian;
Markov
chain
Monte
Carlo
(MCMC)
simulation;
nonlinear
models;
normal
dynamic
linear
models;
observation
error;
population
dynamics;
state-space;
time-series;
uncertainty.
V3A
1V1
03c
IS
INTRODUCTION
Population
biologists
use
time-series
data
to
infer
the
factors
that
regulate
natural
populations
(Stenseth
1995,
Stenseth
et
al.
1998,
Bjornstad
et
al.
1999)
and
to
determine
when
populations
may
be
at
risk
of
ex-
tinction.
Inference
on
the
dynamics
of
natural
popu-
lations
is
based
on
estimated
model
parameters
and
should
incorporate
the
uncertainty
in
these
parameters.
Many
analyses,
however,
ignore
the
multiple
sources
of
stochasticity
that
commonly
impact
population
time-
series
(Dennis
and
Taper
1994,
Bjornstad
and
Grenfell
2001).
These
analyses
include
process
error,
accounting
for
the
fact
that
population
dynamics
is
not
a
deter-
ministic
process
and
that
the
model
for
the
process
may
be
misspecified,
but
they
often
ignore
uncertainty
due
to
errors
in
the
observations.
In
this
way,
the
true
abun-
dances
through
time
is
an
unobserved
("latent")
state
variable
on
which
the
data
provides
imperfect
infor-
mation.
DeValpine
and
Hastings
(2002)
extend
stan-
dard
population
models
to
include
both
observation
error
and
process
error.
Using
maximum
likelihood
techniques,
they
show
that
ignoring
observation
error
can
lead
to
inaccurate
estimates
of
population
param-
eters.
Their
approach
provides
a
more
flexible
method
than
classical
error-in-variables
models
(Carpenter
et
al.
1994).
We
apply
an
alternative
method
for
incor-
porating
observation
error
in
models
of
population
dy-
namics,
which
is
based
on
the
Bayesian,
rather
than
frequentist,
paradigm
of
statistics.
Frequentist
meth-
ods,
sometimes
referred
to
in
the
statistical
literature
Manuscript received
3
December
2001;
revised
2
July
2002;
accepted
18
July
2002;
final
version
received
21
October
2002.
Corresponding
Editor:
0.
N.
Bjornstad.
For
reprints
of
this
Spe-
cial
Feature,
see
footnote
1,
p.
1349.
E-mail:
as
classical
methods,
treat
parameters
as
fixed,
un-
known
constants.
Alternatively,
under
the
Bayesian
paradigm,
parameters
are
viewed
as
random
variables;
the
data
are
used
to
update
one's
beliefs
about
the
dis-
tribution
of
the
model
parameters.
The
coherence
of
the
Bayesian
method
provides
a
straightforward
way
to
account
for
observation
error
in
addition
to
process
error.
The
state-space
model
framework
provides
a
struc-
ture
for
extending
time-series
models
to
include
both
observation
and
process
error.
The
data
are
assumed
to
arise
from
an
unobserved
state
variable
that
represents
the
"true"
dynamic
process.
This
underlying
variable
evolves
over
time
by
a
process
model
that
explicitly
models
process
error.
The
model
for
the
relationship
between
the
actual
data
and
the
state
variable
incor-
porates
observation
error.
Bayesian
state-space
models
with
linear
structure
and
normal
error
distributions
al-
low
entirely
analytic
results;
we
review
the
relevant
expressions.
We
discuss
extensions
and
alternative
ap-
proaches
for
models
with
nonlinear
structure
and
non-
normal
error
distributions,
including
Markov
chain
Monte
Carlo
(MCMC)
posterior
simulation.
The
co-
herence
of
the
Bayesian
paradigm
allows
for
more
straightforward
inference
than
do
maximum
likelihood
methods
and
numerical
quadrature
techniques
required
for
classical
inference.
In
the
Bayesian
setting,
state
variables
can
be
treated
as
parameters,
and
full
pos-
terior
inference
can
be
performed
as
if
they
are
time-
invariant
parameters.
We
describe
how
to derive
full
posterior
distributions
for
the
state
variables/parame-
ters
and
time
invariant
parameters.
Throughout,
we
use
models
for
density-dependent
population
dynamics
as
illustrative
examples.
Bayesian
statistics
and
Bayesian
state-space
models
are
not
new
to
the
ecology
literature.
Soudant
et
al.
1395
SPECIAL
F
me
1396
CATHERINE
CALDER
ET
AL.
Ecology,
Vol.
84,
No.
6
(1997)
developed
a
Bayesian
dynamic
model
for
phy-
toplankton
time
series
that
allows
for
time-varying
in-
fluence
of
the
covariates.
Cottingham
and
Schindler
(2000)
use
a
Bayesian
dynamic
linear
model
to
model
how
phytoplankton
respond
to
pulsed
nutrient
loading.
State-space
models
have
been
used
in
the
fisheries
lit-
erature
(e.g.,
see
references
in
Millar
and
Meyer
[2000a]).
For
example,
Millar
and
Meyer
(2000a)
and
Meyer
and
Millar
(2000b)
introduce
a
Bayesian
non-
linear
state-space
model
to
incorporate
more
realism
into
fish
stock
assessment
models.
In
this
article,
we
show
how
standard
population
models
can
be
extended
to
the
state
space
framework
in
order
to
include
multiple
sources
of
error.
Rather
than
focusing
on
a
particular
model,
we
concentrate
on
general
techniques
necessary
to
fit
Bayesian
state-space
models.
The
relative
simplicity
of
the
methods
de-
scribed
here
have
broad
application
and
lead
to
state-
ments
of
uncertainty
that
take
a
more
comprehensive
accounting
of
variability.
STATE-SPACE
MODELING
The
standard
models
used
to
describe
changes
in
population
size
have
the
form
y
t
=
G(D
t
_
i
)
+
w„
where
y
t
is
a
function
of
the
size
of
the
population
at
time
t
and
D
t
represents
all
of
the
observed
population
counts
up
until
time
t,
i.e.,
D
t
=
fYi,
y
,
,
.
The
function
2
GO
provides
a
deterministic
relationship
between
the
size
of
the
population
at
time
t
and
its
size
in
the
past.
GO
will
typically
also
be
a
function
of
unknown
pa-
rameters
that
can
be
estimated
from
data.
The
term
w
t
represents
the
process
error
and
depends
intrinsically
on
the
function
GO;
it
accounts
for
the
variability
in
the
size
of
the
population
that
cannot
be
captured
by
GO.
Population
census
data
come
from
trapping,
counts,
photographs,
and
so
forth
and
are
typically
observed
with
measurement
error.
Counts
are
rarely
equal
to
the
true
size
of
the
population.
The
error
term
w
t
in
the
model
for
population
dynamics
described
in
the
pre-
vious
paragraph
subsumes
both
process
error
and
mea-
surement
error.
It
is
often
useful
to
extend
this
model
to
a
state-space
model
that
explicitly
separates
multiple
sources
of
stochasticity.
DeValpine
and
Hastings
(2002)
describe
how
inference
can
be
performed
for
state-space
models
within
the
frequentist
framework.
We
focus
on
Bayesian
state-space
models.
The
Bayes-
ian
state-space
model
is
based
on
the
Kalman
filter
(Kalman
1960),
which
is
a
popular
technique
used
in
engineering
and
statistical
quality
control.
While
not
inherently
a
Bayesian
technique,
the
Kalman
filter
pro-
vides
a
method
for
forecasting
that
is
consistent
with
the
theory
of
Bayesian
inference.
Harrison
and
Stevens
(1976)
discuss
the
principles
of
Bayesian
forecasting
and
its
relationship
to
the
Kalman
filter.
Meinhold
and
Singpurwalla
(1983)
present
a
less
technical
version
of
these
issues.
Yt,1
X
I+
FIG.
1.
This
diagram
represents
the
conditional
indepen-
dence
structure
of
a
state-space
model.
Each
of
the
x's,
given
the
values
of
the
surrounding
nodes,
is
conditionally
inde-
pendent
of
the
rest
of
the
graph.
The
standard
framework
for
a
Bayesian
state-space
model
is
as
follows:
Observation
equation
y
t
=
F(x
t
)
+
v„
vt
N[0,
V]
Process
equation
x
t
=
G(x
t
_
i
)
+
w„
wt
N[O,
x
o
=
N[m
o
,
C
o
].
(1)
In
Eq.
1,
y
t
is
a
function
of
the
value
of
the
time
series
at
time
t;
x
t
is
an
unknown
underlying
state
variable,
e.g.,
the
log
of
population
density
at
time
t,
that
is
propagated
through
time
by
the
function
GO.
It
rep-
resents
the
true
size
of
the
population
at
time
t.
The
function
FO
models
the
deterministic
relationship
be-
tween
the
underlying
state
variable
and
the
observa-
tions,
y
t
.
If
the
data
are
believed
to
be
unbiased
esti-
mates
of
y,
FO
can
be
taken
to
be
the
identity
(Coulson
et
al.
2001).
The
variable
v
t
represents
the
measurement
or
observation
error.
It
is
modeled
using
a
normal
dis-
tribution
with
mean
0
and
a
known
variance,
V.
This
normality
assumption
is
not
necessary;
it
is
possible
under
the
Bayesian
framework
to
perform
inference
assuming
any
error
distribution.
Generalization
to
un-
known
V
is
straightforward
(West
and
Harrison
1997).
Bjornstad
et
al.
(1999)
suggest
a
Poisson
counting
pro-
cess
for
the
observation
error.
The
variable
w
t
repre-
sents
the
process
error.
It
is
also
usually
assumed
to
come
from
a
normal
distribution
with
mean
0
and
var-
iance
W.
In
the
Bayesian
framework,
it
is
necessary
to
specify
the
prior
distribution
of
the
initial
latent
state
variable
x
o
.
We
assume
that
it
comes
from
a
normal
distribution
with
known
mean,
m
o
,
and
variance,
C
o
.
Standard
state-space
models
further
assume
that
the
error
terms
are
conditionally
independent.
Fig.
1
shows
the
structure
of
a
state-space
model.
Each
of
the
x's
given
the
values
of
the
surrounding
nodes
is
condi-
tionally
independent
of
the
rest
of
the
graph.
Normal
dynamic
linear
models
When
the
functions
FO
and
GO
are
linear
and
the
error
distributions
are
normal,
the
Bayesian
state-space
model
is
termed
a
normal
dynamic
linear
model
(NDLM).
The
functions
FO
and
GO
are
replaced
by
the
constants
f
and
g
that
premultiply
x
t
in
the
obser-
vation
equation
and
x
t
_
1
in
the
process
equation,
re-
X
t-i
x
t
June
2003
ECOLOGICAL
INFERENCE
AND
FORECASTS
1397
a
0.6
-
0.4
-
0.2
-
0.0
-
b
2
4
6
10
0.6
0.4
0.2
00
0
2
4
6
8
10
C
0.6
0.4
0.2
0.0
2
4
6
8
10
x
FIG.
2.
This
figure
demonstrates
how
the
forward
filtering
backward
smoothing
(FFBS)
algorithm
updates
posterior
distributions
as
it
processes
data
sequentially.
The
details
of
the
model
and
data
are
given
in
State-space
modeling:
Normal
dynamic
linear
models.
The
two
data
points,
y,
and
y
2
,
are
represented
by
the
black
diamonds
in
(a)
and
(b),
respectively.
The
lines
represent
the
following
distributions:
(a)
dashed
line,
p(x);
solid
line,
p(x
l
iy,);
(b)
dashed
line,
p(x
2
1y
i
);
solid
line,
/*AY,
Y2);
(c)
solid
line,
p(xik,
Y2).
3N
111V3.J
1V1
03c
IS
spectively.
In
the
case
when
either
the
observation
equation
or
the
process
equation
is
linear
with
an
in-
tercept
term,
the
underlying
state
parameter
is
replaced
by
the
vector
{1,
xj.
The
posterior
distribution
of
the
x
t
's
in
an
NDLM
can
be
found
analytically
by
taking
advantage
of
the
model's
conditional
independence
structure.
As
before,
Dt
=
Y2,
yt)
The
desired
posterior
distributions
are
p(x,ID
T
)
for
t
=
1,
2,
. . .
T,
but
for
the
moment
we
consider
finding
p(x
t
iD
t
)
for
all
t.
Using
Bayes'
The-
orem,
p(4D
t
)
PO'bcp
Dt-OP(xtlpt-)
(2)
where
p(yjx„
D")
is
the
likelihood
of
y
t
given
all
of
the
past
data
and
p(x
t
iD
t
,)
is
prior
density
of
x
t
con-
ditional
on
the
past
data
values.
Given
the
linear
structure
of
the
process
equation,
the
fact
that
p(x"ID
t
_
i
)
is
normal
allows
p(x
t
ID"),
the
second
factor
on
the
right
side
of
Eq.
2,
to
be
found
in
closed
form.
This
is
done
simply
by
updating
the
moments
of
p(x"ID")
according
to
the
process
equa-
tion.
Once
the
next
data
point,
y,,
is
processed,
the
prior
of
x
t
given
D"
is
updated
to
the
posterior
distribution
P(xtiDt).
In
this
manner,
the
distributions
p(x
t
iD
t
)
are
com-
puted
sequentially.
We
begin
with
p(x
0
),
which
we
as-
sume
to
be
normal
with
mean
m
o
and
variance
C
o
.
This
procedure,
known
as
the forward
filtering
algorithm,
due
to
Carter
and
Kohn
(1994)
and
Friihwirth-Schnatter
(1994),
is
Theorem
4.1
in
West
and
Harrison
(1997).
The
forward
filtering
algorithm
defines
a
procedure
for
sequentially
determining
the
posterior
distribution
of
each
of
the
x
t
's
given
D
1
.
But
we
desire
the
posterior
distributions
p(x
t
ID
T
)
for
t
=
1,
2,
. . .
T,
i.e.,
the
pos-
terior
distributions
of
the
latent
variables
given
all
of
the
data.
These
distributions
can
be
found
by
recur-
sively
updating
the
moments
ofp(x,ID
T
).
The
formulas
for
this
backward
smoothing
algorithm
are
given
in
Theorem
4.4
of
West
and
Harrison
(1997).
Together,
these
two
algorithms
are
know
as
the
forward
filtering
backward
smoothing
(FFBS)
algorithm
and
are
due
to
Carter
and
Kohn
(1994)
and
Friihwirth-Schnatter
(1994).
It
is
a
two-step
algorithm.
First,
the
forward
filtering
algorithm
is
run
on
the
data.
Then,
using
the
final
posterior
distribution
derived
by
the
first
algo-
rithm,
p(x
T
ID
T
),
the
backward
smoothing
algorithm
is
used
to
recursively
find
p(x
t
ID
T
)
for
all
t.
Both
of
these
algorithms
simply
update
the
moments
of
normal
dis-
tributions
using
the
structure
of
the
state-space
model
(Eq.
1).
For
details
on
the
algorithm,
see
the
appendix.
Fig.
2
demonstrates
how
the
FFBS
algorithm
updates
posterior
distributions
as
it
processes
data
sequentially.
For
simplicity,
we
assume
that
the
data
arise
from
a
NDLM
as
in
Eq.
1
with
the
function
FO
and
GO
taken
to
be
the
identity
and
V
=
W
=
1.
1398
CATHERINE
CALDER
ET
AL.
Ecology,
Vol.
84,
No.
6
SPECIAL
FEATURE
Assume
we
have
two
data
points,
y
i
=
3
and
y
2
=
8,
and
that
our
prior
on
x
0
is
9V(5,
3).
The
dashed
line
in
Fig.
2a
shows
the
prior
distribution
of
x
1
given
no
data,
i.e.,
p(x
1
)
=
9V(5,
4).
The
moments
of
this
distri-
bution
are
computed
using
the
formulas
in
step
b
of
the
forward
filtering
(FF)
algorithm.
This
distribution
can
be
interpreted
as
the
prior
distribution
for
x
1
.
The
point
on
this
graph
represents
the
first
data
point,
y
i
=
3,
and
the
solid
line
represents
the
posterior
distribution
of
x
1
given
the
first
data
point,
p(x
i
ly
i
)
=
N(3.4,
0.8).
The
moments
of
this
distribution
are
calculated
using
step
d
of
the
FF
algorithm.
Since
y
i
is
less
than
the
prior
mean
for
x
1i
and
FO
is
the
identity,
y
1
pulls
the
distribution
of
x
1
to
the
left.
The
dashed
line
in
Fig.
2b
represents
the
prior
dis-
tribution
of
x
2
given
only
the
first
data
point,
p(x
2
ly
j
)
=
N(3.4,
1.8).
Its
moments
are
calculated
by
process-
ing
the
posterior
distribution
of
x
1
given
y
i
through
the
process
equation
(using
step
b
of
the
FF
algorithm).
Since
GO
is
the
identity,
this
process
stretches
out
the
distribution
represented
by
the
solid
line
in
the
first
graph
but
keeps
it
centered
at
3.4.
This
distribution
can
be
thought
of
as
the
prior
distribution
of
x
2
given
only
y
i
and
can
be
updated
to
the
posterior
distribution
of
x
2
given
both
y
i
and
y
2
,
p(x
2
k,
y
2
)
using
step
d
of
the
FF
algorithm.
Since
y
2
=
8,
this
posterior
distribution
is
to
the
right
of
the
prior
distribution.
We
have
computed
the
posterior
distribution
for
x
2
given
all
the
data
in
this
simple
example,
but
we
must
use
the
backward
smoothing
(BS)
algorithm
to
derive
the
posterior
distribution
of
x
1
given
all
of
the
data.
Using
the
BS
algorithm,
p(x
i
ly
i
,
y
2
)
=
N(4.7,
0.571).
If
we
had
more
data,
we
would
follow
the
same
pro-
cedure
adding
an
additional
step
in
both
the
FF
and
BS
algorithms
for
each
data
point.
Nonlinear
dynamic
models
Inference
on
the
latent
state
variables
in
nonlinear
dynamic
models
is
slightly
more
complicated.
The
pos-
terior
distributions
of
the
x
t
's
cannot
be
found
in
closed
form
as
they
are
above.
Instead,
we
rely
on
numerical
techniques
to
explore
their
posteriors.
The
standard
method
for
performing
this
inference
is
the
Markov
chain
Monte
Carlo
(MCMC)
algorithm
that
approxi-
mates
the
desired
posterior
distributions
(Gilks
et
al.
1996).
The
algorithm
constructs
a
Markov
chain
with
the
posterior
distribution
as
its
stationary
distribution.
The
chain
is
run
for
sufficient
time
allowing
it
to
ap-
proach
its
stationary
distribution.
Following
conver-
gence,
samples
approximate
the
posterior
distribution.
The
Gibbs
sampler
is
a
version
of
MCMC
simulation
that
we
use
to
approximate
the
posterior
distributions
of
the
x
t
's.
The
algorithm
works
by
iteratively
sampling
from
each
of
the
full
conditional
distributions
of
each
of
the
states
given
the
current
value
of
the
other
var-
iables.
For
example,
if
we
want
to
sample
from
the
joint
distributionp(A,
B,
C),
we
could
construct
a
Gibbs
sampler
as
follows:
1)
sample
Am
p(A1B
0-1
),
Ci_
1)
)
2)
sample
B(0
p(BIA(°,
Co
-
'))
3)
sample
0"
p(CIA(0,
B(0)
4)
repeat
steps
1-3
many
times,
incrementing
i
after
each
iteration.
When
sampling
from
these
distributions,
we
condition
on
the
currently
imputed
values,
i.e.,
the
most
recent
parameter
values.
In
the
nonlinear
dynamic
model,
each
of
the
states,
given
the
next
state
and
the
previous
state,
is
indepen-
dent
of
the
other
states.
The
conditional
independence
structure
shown
in
Fig.
1
still
applies.
Thus,
we
need
only
condition
on
neighboring
states.
The
Gibbs
sam-
pler
for
the
nonlinear
dynamics
model
works
as
fol-
lows:
1)
Choose
initial
values
for
each
of
the
latent
state
variables
and
denote
them
as
40),
xr,
. . .
40).
2)
For
each
1
t
T,
sequentially
draw
a
sample
from
P
(xt
x11
)
,
x
.
4
0
1,
xM
l)
,
. .
.
,
4
-1
),
D
t
)
(3)
P(xt
I
4
2
1,
xV
)
,
yr)
(4)
P(YtIxt)P(xtlx
V,71
1)
xt)
(5)
for
each
t.
3)
Repeat
step
2,
I
times.
The
Gibbs
sampler
above
requires
simulation
from
the
distribution
that
is
proportional
to
P(yt
I
xt)P(xt
I
x
V
2
1)/Xx
`+
-
1
1)
1xt).
This
can
be
done
using
a
Metropolis-Hastings
(M-H)
step,
or
in
some
cases
a
special
case
of
M-H
known
as
a
Metropolis
step
(Gelfand
and
Smith
1990,
Tierney
1994).
See,
for
example,
Gilks
et
al.
(1996)
for
a
gen-
eral
review
of
the
algorithm.
Carlin,
Polson,
and
Stoffer
(1992)
develop
an
efficient
M-H
algorithm
for
sam-
pling
the
latent
state
parameters
in
a
nonlinear
dynamic
model.
The
samples
xp
for
i
>
B„
where
B
x
is
the
number
of
iterations
until
the
Markov
chain
converges,
will
be
samples
from
p(x,ID
T
).
Inference
for
parametric
state-space
models
The
posterior
distribution
of
parameters
of
the
model
for
population
dynamics
captured
in
GO
are
usually
the
focus
of
a
time-series
analysis.
An
extended
Gibbs
sampler
can
be
used
to
find
the
posterior
distributions
of
model
parameters
and
the
latent
state
variables
si-
multaneously.
After
assigning
initial
values
to
all
pa-
rameters,
we
can
sample
from
the
full
conditional
dis-
tribution
of
the
model
parameters
given
the
currently
imputed
values
of
the
x
t
's.
We
then
sample
from
the
full
conditional
distributions
of
the
x
t
's
conditioning
on
the
currently
imputed
values
of
the
model
parameters.
The
observation
and
process
error
variances,
V
and
W,
are
often
unknown.
Their
posterior
distributions
can
1
-
—V
T
-
XT-1
2
W
-
a
June
2003
ECOLOGICAL
INFERENCE
AND
FORECASTS
1399
be
recovered
by
sampling
from
their
full
conditional
distributions
within
the
Gibbs
sampler.
If
a
is
a
param-
eter
of
the
function
GO,
the
Gibbs
sampler
iterates
over
the
following
steps:
1)
Sample
from
p(x
t
lx,_„
x„„
a,
V,
W,
y)
for
t
=
1,
2,
. . .
,
T
using
the
M-H
algorithm
or
directly
using
FFBS.
2)
Sample
from
p(alxi,
x2,
,
xT,
V,
W,
D
T
)
using
the
M-H
algorithm
or,
if
possible,
directly.
3)
Sample
fromp(Vixi,
x2,
,
xT
,
a,
W,
D
T
)
directly.
4)
Sample
fromp(Wlxi,
x2,
,
xT,
a,
V,
D
T
)
directly.
In
each
of
the
steps,
the
values
of
the
variables
to
the
right
of
the
conditional
bar
are
assumed
to
be
the
currently
imputed
values
of
each
of
the
parameters.
It
is
also
necessary
to
specify
prior
distributions
for
a,
V,
and
W
if
they
are
assumed
to
be
unknown.
These
prior
distributions
must
be
included
when
determining
each
of
the
full
conditional
distributions.
In
the
next
section,
we
show
the
results
of
a
simulation
study
for
a
state-space
model.
SIMULATION
STUDY
We
use
a
simulation
study,
similar
to
the
one
in-
cluded
in
deValpine
and
Hastings
(2002),
to
demon-
strate
the
importance
of
including
observational
error
and
additionally
to
demonstrate
how
this
particular
model,
when
framed
as
a
Bayesian
dynamic
model,
can
be
fit.
We
focus
on
two
models,
the
observation-error
model
and
the
no-observation-error
model.
Both
mod-
els
are
based
on
a
common
model
for
density
depen-
dence,
the
Ricker
model
(e.g.,
Dennis
and
Taper
1994).
The
observation-error
model
explicitly
incorporates
observation
error
into
the
Ricker
model,
while
the
no-
observation-error
model
lumps
process
error
and
ob-
servation
error
together.
The
observation-error
model
can
be
written
as:
log()
=
x,
+
v,
v
t
-
N[0,
V]
x,
=
x,
1
+
a
+
bey
-1
+
w,
wt
-
9L[0,
W]
x
o
-
Admo,
Co]
V
-
IG[a
v
,
13,,]
W
-
IG
[a
w
,
o
w
]
p(a,
b)
-
NAP,
']•
IG
represents
the
inverse
gamma
distribution,
and
9V
2
represents
the
bivariate
normal
distribution.
The
no-observation-error
model
can
be
written
as:
Yt
i
`
=
.4-1
±
a
±
best
s
-1
+
w
t
wt
-
N[O,
W]
Y:
-
N[mo,
Co]
W
-
IG[a
w
,
o
w
]
P(a
,
b)
-
NAP,
]
and
y;K
is
defined
to
be
log(y).
First,
we
simulate
data
from
the
observation-error
mod-
el
with
the
parameter
a
fixed
at
0.1
and
the
parameter
b
fixed
at
-0.01.
This
is
the
"correct"
model
for
the
data.
We
then
find
the
posterior
distributions
of
the
model
pa-
rameter
b
under
both
models
and
compare
them.
For
the
observation-error
model,
we
use
a
Gibbs
sampler
to
iteratively
sample
from
the
full
conditional
distributions
of
each
of
the
parameters
in
the
model.
Each
of
the
x
t
's
are
sampled
using
a
M-H
step
and
the
parameters
a,
b,
V,
and
W
are
sampled
directly
given
the
current
values
of
the
x
t
's
using
standard
linear
model
theory.
The
full
conditional
distributions
are
[
1
exPL
-
2
c(x0
mo>
o
2
1
,
-
ex
i
-
x
o
-
a
-
bex
0
)
2
p(xt1
—)
cc
1
expl
--
2v
[log(Yt)
xt]z
-
2W
(xt
-
xt-1
-
a
-
1
best
-1)z
bext)
2
}
1
(
2
Te
t±i
xt
a
for
t
=
1,
2,
..
.
,
T
-1
1
Axil
—)
cc
expl
flog(.YT)
xTP
-
bexT-92}
P[(I'
,
b)i
-
]
cc
Nz[alt
Si
where
S
=
(
TT
,
YX
+
-1
)
1
m
=
S
(
1
Tv
X'Z
+
-1
1,1,)
and
X
and
Z
are
defined
to
be
a
e
xo
0
5
e
xi
1
-
]
x=
ii.
.
0
If
:
1
-
]
,
e
xT_I,
p(VI-)
oc
IGIct,
+
',
0,
+
p(WI
-)
«
IG[
a
w
+
2,
1
T
13
w
+
- E
(x
t
x
t
_
i
a
-
best-1)z
2
t=i
We
can
sample
the
posterior
distributions
for
a,
b,
and
W
directly
under
the
no-observation-error
model
without
the
need
for
MCMC
using
standard
results
from
linear
regression.
We
generate
samples
from
these
posteriors
of
the
same
size
as
the
samples
generated
from
the
MCMC
algorithm
used
to
analyze
the
obser-
vational
error
model.
We
use
three
different
values
of
V.
For
each
value
of
V,
we
simulate
10
random
data
sets
and
summarize
the
resulting
posterior
distributions
for
b
under
both
P(x01
)cc
n:
1
-
]
xv_i
El
1
T
2
E
[log(y
t
)
-
x
t
]
2
}
t-i
3N
111V3.J
1V1
03c
IS
0.001
0.01
0.05
0.001
0.01
1400
CATHERINE
CALDER
ET
AL.
Ecology,
Vol.
84,
No.
6
SPECIAL
FEATURE
0.05
0.0
Observation
-error
model
0.05
-
0.0
-
No-observation
-error
model
-0.05
-0.05
-
-0.1
-0.1
-
-0.15
-0.15
-
-0.2
-0.2
-
models
(Fig.
3).
Because
we
know
that
the
observation-
error
model
is
the
correct
model
for
the
data,
differ-
ences
in
the
posterior
distributions
demonstrate
mis-
leading
inference
that
results
from
combining
the
two
types
of
error.
In
this
case,
the
posterior
distributions
under
the
observation-error
model
are
mostly
tighter
than
under
the
no-observation-error
model.
In
other
words,
in
the
no-observation-error
model,
the
two
types
of
error
are
inappropriately
combined,
leading
to
overly
broad
credible
intervals.
It
also
appears
from
Fig.
3
that
the
posterior
median
becomes
a
considerably
more
biased
estimate
of
b
as
the
true
value
of
V
used
to
simulate
the
data
increases.
These
results
are
not
meant
TABLE
1.
Summary
of
Bayes
factors
for
50
trials.
Bayes
Number
factor
Evidence
against
of
trials
<1
none
1-3
not
worth
more
than
a
bare
mention
3-20
positive
0-150
strong
>150
very
strong
Notes:
The
true
values
of
the
parameters
in
the
Ricker
model
are
V
=
0.01,
W
=
0.01,
a
=
0.1,
b
=
-0.01.
The
Bayes
factor
is
computed
by
taking
the
ratio
of
the
p
(datalH
th
,)
to
p
(datalH,,"„
)
.
The
interpretations
of
the
Bayes
factors
are
taken
from
Kass
and
Raftery
(1995).
H
th
model
with
observation
error;
1-1,,
oth
model
without
observation
er-
ror.
0.05
to
be
generalized;
if
we
use
a
different
model
or
dif-
ferent
parameter
values,
we
expect
different
results.
The
simulation
study
presented
here
demonstrates
that
for
one
particular
example
it
makes
a
difference
wheth-
er
or
not
observation
error
is
modeled
explicitly.
In
addition
to
comparing
posterior
parameters
under
the
two
models,
we
also
computed
Bayes
factors,
which
assess
the
weight
of
the
evidence
in
favor
of
one
model
as
opposed
to
the
other
model.
(See
Kass
and
Raftery
[1995]
for
an
overview
of
the
technique.)
If
H
1
rep-
resents
the
hypothesis
that
the
data
are
generated
by
model
1
and
H2
represents
the
hypothesis
that
the
data
are
generated
by
model
2,
the
Bayes
factor
B12
is
de-
fined
as
follows:
B12'
=
p
(data
I
H2)*
(6)
P
(data
I
H
i
)
The
Bayes
factor
represents
the
ratio
of
the
posterior
odds
of
hypothesis
H
1
to
its
prior
odds.
Kass
and
Raf-
tery
(1995)
illustrate
how
to
calculate
and
interpret
Bayes
factors
in
various
settings.
A
recipe
for
approx-
imating
Bayes
factors
using
output
from
an
MCMC
simulation
is
due
to
Newton
and
Raftery
(1994).
The
numerator
and
denominator
of
the
Bayes
factor,
Eq.
6,
can
be
approximated
as
-1
"`
fi(data
I
IL)
=
1
1
P[dataI0(
0
,H]
-1
.
This
is
the
harmonic
mean
of
the
likelihood
of
the
data,
0
0
0
1
49
FIG.
3.
The
vertical
lines
represent
95%
credible
intervals
of
the
posterior
distribution
of
the
parameter
b
under
the
two
different
models
for
10
data
sets
simulated
from
the
observation-error
model.
(A
95%
Bayesian
credible
interval
is
the
range
from
the
0.025
quantile
to
the
0.975
quantile
of
a
posterior
distributions.
It
can
be
approximated
using
samples
from
a
posterior
distribution.)
The
posterior
medians
are
denoted
by
the
dash
through
the
vertical
lines.
This
simulation
study
was
repeated
for
V
=
0.001,
0.01,
and
0.05,
and
the
corresponding
posterior
summaries
are
grouped
according
to
the
value
of
V
used
to
simulate
the
data
set.
The
true
value
of
b
is
represented
by
the
horizontal
line
at
b
=
-0.01.
For
each
of
the
trials,
the
value
of
W
used
to
simulate
the
data
was
0.01
and
the
length
of
the
simulated
data
sets,
7',
was
60.
June
2003
ECOLOGICAL
INFERENCE
AND
FORECASTS
1401
given
the
samples
of
the
parameters
taken
in
the
MCMC
algorithm.
It
can
be
shown
that
p(datallt)
con-
verges
almost
surely
to
p(datalH)
as
m
00.
Table
1
summarizes
the
Bayes
factors
computed
in
this
way
for
50
trials;
the
two
models
are
fit
and
the
Bayes
factor
comparing
them
is
computed
for
each
trial.
Also
pro-
vided
in
the
table
are
interpretations
of
the
Bayes
fac-
tors.
Table
1
shows
that
all
but
one
of
the
50
trials
provided
strong
positive
evidence
against
the
null
hy-
pothesis,
i.e.,
against
the
hypothesis
that
the
data
were
generated
from
the
no-observation-error
model.
DISCUSSION
State-space
models
provide
a
framework
for
incor-
porating
observational
error
into
dynamic
models
of
population
size.
Bayesian
state-space
models
and
pro-
cedures
used
to
perform
inference
on
them
are
general
enough
to
be
used
on
more
complicated
models
than
the
one
discussed
here,
including
those
having
higher
embedding
dimension
and
thresholds.
We
need
only
to
iteratively
sample
from
the
full
conditional
distribu-
tions
of
each
model
parameter,
which
can
be
done
di-
rectly
using
the
Metropolis-Hastings
algorithm.
When
performing
a
statistical
analysis
on
data
con-
taining
observation
error,
we
recommend
considering
a
state-space
model
because
it
explicitly
separates
pro-
cess
error
from
observation
error
and
therefore
affords
a
more
detailed
understanding
of
the
posterior
distri-
bution.
Our
simulation
study
illustrates
that
ignoring
observational
error
can
result
in
misleading
inference,
especially
in
estimates
of
posterior
uncertainty
of
mod-
el
parameters.
The
results
of
this
specific
simulation
study
are
not
general.
For
others
models,
the
results
of
ignoring
observation
error
will
be
different;
we
simply
show
that
explicitly
modeling
observation
error
can
make
a
difference.
There
is
freely
available
software
that
can
be
used
to
fit
Bayesian
state-space
models.
BUGS
(Spiegel-
halter
et
al.
1996), Bayesian
inference
using
Gibbs
sam-
pling,
takes
advantage
of
the
acyclical
graphical
struc-
ture
of
a
Bayesian
model
to
determine
the
full
condi-
tional
distributions
required
to
construct
a
Gibbs
sam-
pler.
See
Meyer
and
Millar
(1999)
for
an
illustration
of
the
BUGS
software
in
the
context
of
fitting
a
non-
linear
Bayesian
state-space
model
to
analyze
fisheries
stock
assessments.
Also
freely
available
is
the
BATS
software,
which
fits
a
variety
of
Bayesian
dynamic
models
as
described
in
West
and
Harrison
(1997).
In-
structions
for
this
package
are
available
in
Pole
et
al.
(1994).
General
purpose
statistical
software
packages
can also
be
used
to
fit
state-space
models.
We
propose
fitting
Bayesian
state-space
models
as
opposed
to
classical
state-space
when
it
is
necessary
to
explicitly
include
observation
error
in
a
time-series
model.
The
coherence
of
Bayesian
paradigm
allows
straightforward
inference;
both
the
time-invariant
pa-
rameters
and
the
underlying
variables
can
be
treated
in
the
same
way.
This
is
true
regardless
of
how
the
anal-
ysis
is
performed.
LITERATURE
CITED
Bjornstad,
0.
N.,
J.-M.
Fromentin,
N.
C.
Stenseth,
and
J.
Gjosmter.
1999.
Cycles
and
trends
in
cod
populations.
Ecology
96:5066-5071.
Bjornstad,
0.
N.,
and
B.
T.
Grenfell.
2001.
Noisy
clockwork:
time
series
analysis
of
population
fluctuations
in
animals.
Science
293:638-643.
Carlin,
B.
P.,
N.
G.
Polson,
and
D.
S.
Stoffer.
1992.
A
Monte
Carlo
approach
to
nonnormal
and
nonlinear
state-space
modeling.
Journal
of
the
American
Statistical
Association
87:493-500.
Carpenter,
S.
R.,
K.
L.
Cottingham,
and
C.
A.
Stow.
1994.
Fitting
predator-prey
models
to
time
series
with
observa-
tion
errors.
Ecology
75:1254-1264.
Carter,
C.
K.,
and
R.
Kohn.
1994.
On
Gibbs
sampling
for
state
space
models.
Biometrika
81(3):541-553.
Cottingham,
K.,
and
D.
E.
Schindler.
2000.
Effects
of
grazer
community
structure
on
phytoplankton
response
to
nutrient
pulses.
Ecology
81:183-200.
Coulson,
T.,
E.
A.
Catchpole,
S.
D.
Albon,
B.
J.
T.
Morgan,
J.
M.
Pemberton,
T.
R.
Clutton
Brock,
M.
J.
Crawley,
and
B.
T.
Grenfell.
2001.
Age,
sex,
density,
winter
weather,
and
population
crashes
in
soay
sheep.
Science
292:1528-
1531.
Dennis,
B.,
and
M.
L.
Taper.
1994.
Density
dependence
in
time
series
observations
of
natural
populations-estimation
and
testing.
Ecological
Monographs
64:205-224.
deValpine,
P.,
and
A.
Hastings.
2002.
Fitting
population
mod-
els
with
process
noise
and
observation
error.
Ecological
Monographs
72:57-76.
Frithwirth-Schnatter,
S.
1994.
Data
augmentation
and
dy-
namic
linear
models.
Journal
of
Time
Series
Analysis
15:
183-202.
Gelfand,
A.
E.,
and
A.
F.
M.
Smith.
1990.
Sampling-based
approaches
to
calculating
marginal
densities.
Journal
of
the
American
Statistical
Association
85:129-133.
Gilks,
W.
R.,
S.
Richardson,
and
D.
J.
Spiegelhalter.
1996.
Markov
chain
Monte
Carlo
in
practice.
Chapman
and
Hall,
New
York,
New
York,
USA.
Harrison,
P.
J.,
and
D.
F.
Stevens.
1976.
Bayesian
forecasting
(with
discussion).
Journal
of
the
Royal
Statistical
Society,
Series
B
38:205-247.
Kalman,
R.
E.
1960.
A
new
approach
to
linear
filtering
and
prediction
problems.
Journal
of
Basic
Engineering,
Series
B
82:34-45.
Kass,
R.
E.,
and
A.
E.
Raftery.
1995.
Bayes
factors.
Journal
of
the
American
Statistical
Association
90(430):773-795.
Meinhold,
R.
J.,
and
N.
D.
Singpurwalla.
1983.
Understand-
ing
the
Kalman
filter.
American
Statistician
37(2):123-127.
Meyer,
R.,
and
R.
B.
Millar.
1999.
BUGS
in
Bayesian
stock
assessment.
Canadian
Journal
of
Fisheries
and
Aquatic
Sci-
ences
56:1078-1086.
Millar,
R.
B.,
and
R.
Meyer.
2000a.
State-space
modeling
of
non-linear
fisheries
biomass
dynamics
using
the
Gibbs
sam-
pler.
Applied
Statistics
49:327-342.
Millar,
R.
B.,
and
R.
Meyer.
2000b.
Bayesian
state-space
modeling
of
age
structured
data:
fitting
a
model
is
just
the
beginning.
Canadian
Journal
of
Fisheries
and
Aquatic
Sci-
ences
57:43-50.
Newton,
M.
A.,
and
A.
E.
Raftery.
1994.
Approximate
Bayes-
ian
inference
by
the
weighted
likelihood
bootstrap
(with
discussion).
Journal
of
the
Royal
Statistical
Society,
Series
B
56:3-48.
Pole,
A.,
M.
West,
and
P.
J.
Harrison.
1994.
Applied
Bayesian
forecasting
and
time
series
analysis.
Chapman
and
Hall,
New
York,
New
York,
USA.
3N
111V3.J
1V1
03c
IS
1402
CATHERINE
CALDER
ET
AL.
Ecology,
Vol.
84,
No.
6
Soudant,
D.,
B.
Beliaeff,
and
G.
Thomas.
1997.
Dynamic
linear
Bayesian
models
in
phytoplankton
ecology.
Ecolog-
ical
Modelling
9:161-169.
Spiegelhalter,
D.
J.,
A.
Thomas,
N.
Best,
and
W.
R.
Gilks.
1996.
BUGS
0.5,
Bayesian
inference
using
Gibbs
sam-
pling.
Manual
version
ii.
Medical
Research
Council
Bio-
statistics
Unit,
Institute
of
Public
Health,
Cambridge,
UK.
Stenseth,
N.
C.
1995.
Snowshoe
hare
populations—squeezed
from
below
and
above.
Science
269:1061-1062.
Stenseth,
N.
C.,
W.
Falck,
K.
Chan,
0.
N.
Bjornstad,
M.
O'Donoghue,
H.
Tong,
R.
Boonstra,
S.
Boutin,
C.
Krebs,
and
N.
G.
Yoccoz.
1998.
From
patterns
to
processes:
phase
and
density
dependencies
in
the
Canadian
lynx
cycle.
Pro-
ceedings
of
the
National
Academy
of
Sciences
(USA)
95:
15430-15435.
Tierney,
L.
1994.
Markov
chains
for
exploring
posterior
dis-
tributions.
Annals
of
Statistics
22:1701-1728.
West,
M.,
and
J.
Harrison.
1997.
Bayesian
forecasting
and
dynamic
models.
Second
edition.
Springer-Verlag,
New
York,
New
York,
USA.
APPENDIX
The
forward
filtering
backward
smoothing
algorithm
is
available
in
ESA's
Electronic
Data
Archive:
Ecological
Archives
E084-033-A1.