
https://cpw.cvlcollections.org/files/original/2fd585dbd66cb8975b18c0fa5e0e94e8.pdf
1c089511e42775f13baa36e749768ac5
PDF Text
Text
The research in this publication was partially or fully funded by Colorado Parks and Wildlife.
Dan Prenzlow, Director, Colorado Parks and Wildlife • Parks and Wildlife Commission: Marvin McDaniel, Chair • Carrie Besnette Hauser, ViceChair
Marie Haskett, Secretary • Taishya Adams • Betsy Blecha • Charles Garcia • Dallas May • Duke Phillips, IV • Luke B. Schafer • James Jay Tutchton • Eden Vardy
�Special Invited Article
Received: 12 May 2016,
Revised: 8 June 2016,
Environmetrics
Accepted: 9 June 2016,
Published online in Wiley Online Library: 19 July 2016
(wileyonlinelibrary.com) DOI: 10.1002/env.2402
Hierarchical animal movement models for
populationlevel inference�
b
b
Mevin B. Hootena*, c Frances E. Buderman
,
Brian
M.
Brost
,
Ephraim M. Hanks and Jacob S. Ivand
New methods for modeling animal movement based on telemetry data are developed regularly. With advances in telemetry
capabilities, animal movement models are becoming increasingly sophisticated. Despite a need for populationlevel inference, animal movement models are still predominantly developed for individuallevel inference. Most efforts to upscale
the inference to the population level are either post hoc or complicated enough that only the developer can implement the
model. Hierarchical Bayesian models provide an ideal platform for the development of populationlevel animal movement
models but can be challenging to fit due to computational limitations or extensive tuning required. We propose a twostage
procedure for fitting hierarchical animal movement models to telemetry data. The twostage approach is statistically rigorous and allows one to fit individuallevel movement models separately, then resample them using a secondary MCMC
algorithm. The primary advantages of the twostage approach are that the first stage is easily parallelizable and the second
stage is completely unsupervised, allowing for an automated fitting procedure in many cases. We demonstrate the twostage procedure with two applications of animal movement models. The first application involves a spatial point process
approach to modeling telemetry data, and the second involves a more complicated continuoustime discretespace animal
movement model. We fit these models to simulated data and real telemetry data arising from a population of monitored
Canada lynx in Colorado, USA. Copyright © 2016 John Wiley & Sons, Ltd.
Keywords: hierarchical model; resource selection model; spatial statistics; telemetry data; trajectories
1. INTRODUCTION
The field of movement ecology is booming, in large part, because of the increased availability of telemetry data sources (Cagnacci et al.,
2010). Contemporary telemetry data are acquired via satellite communication devices affixed to individual animals. These devices often
collect many types of data, but most studies are focused on the position data, primarily to learn about environmental influences on individuallevel movement. Many new statistical models for animal trajectories have been proposed in recent years, and they vary in form depending on
the motivation for the project and type of inference desired (Hooten et al., in press). For example, most individualbased statistical models
for telemetry data fall into one of three classes: point process models, discretetime models, or continuoustime models, with each being
appropriate in certain settings (McClintock et al., 2014).
Statistical inference arising from fitting animal movement models to telemetry data is sometimes focused on the individual level. For
example, a movement ecologist might ask how a specific individual animal responded to environmental cues while migrating between
summer and winter home ranges (e.g., Hooten et al., 2010a). However, many animal movement studies are concerned with populationlevel inference. That is, for several individuals, is there evidence of consistent behavioral responses to environmental variables? To obtain
populationlevel inference, the wellaccepted approach is to use a hierarchical model with random effects for individuals that are pooled at
the population level. For example, consider the Bayesian hierarchical model
yj � Œyj jˇ j ; � j �
(1)
*
Correspondence to: M. B. Hooten, U.S. Geological Survey, Colorado Cooperative Fish and Wildlife Research Unit; Departments of Fish, Wildlife, & Conservation
Biology and Statistics, Colorado State University, Fort Collins, CO 80523 U.S.A. Email: mevin.hooten@colostate.edu
a
U.S. Geological Survey, Colorado Cooperative Fish and Wildlife Research Unit; Departments of Fish, Wildlife, & Conservation Biology and Statistics, Colorado
State University, Fort Collins, CO 80523 U.S.A.
b
Department of Fish, Wildlife, and Conservation Biology, Colorado State University
c
Department of Statistics, Pennsylvania State University
d
Colorado Parks and Wildlife
†
This material is associated with the President’s Invited Lecture from the 26th Annual Conference of the International Environmetrics Society (TIES), presented
18 July 2016.
322
Environmetrics 2016; 27: 322–333
Copyright © 2016 John Wiley & Sons, Ltd.
�POPULATIONLEVEL INFERENCE FOR ANIMAL MOVEMENT
Environmetrics
ˇ j � Œˇ j j�ˇ ; † ˇ �
(2)
�ˇ � Œ�ˇ �
(3)
i
h
�1
† �1
ˇ � †ˇ
(4)
� j � Œ� j �
(5)
where yj are measurements associated with each individual j (j D 1; : : : ; J ) and we use ‘Œ: : :�’ to denote a probability distribution
or mass/density function as necessary (Gelfand and Smith, 1990). The priors in (3)–(5) are for the auxiliary datalevel parameters � j ,
populationlevel coefficients �ˇ , and precision matrix † �1
, forming the familiar threelevel hierarchical model (Berliner, 1996). The hierˇ
archical model in (1)–(5) provides a straightforward and intuitive means for obtaining inference for �ˇ , which is the ultimate goal of many
animal movement studies. Similar hierarchical models have become popular, and now standard, tools for obtaining upscaled inference in
many other fields such as atmospheric science (Cressie and Wikle, 2011), ecology (Hobbs and Hooten, 2015), and sociology (Gelman and
Hill, 2006).
The complexity of modern animal movement models makes implementation challenging. Furthermore, increases in the quantity of data
resulting from newer telemetry devices has outpaced computational methods for fitting animal movement models. Animal ecologists may
wish to extend individuallevel models to provide statistically rigorous populationlevel inference, but, in many cases, the algorithms required
to fit such models become prohibitively challenging to program or are too slow in settings with large data sets and/or many individuals. For
example, Hanks et al. (2011) performed a post hoc metaanalysis to obtain populationlevel inference for northern fur seals (Callorhinus
ursinus) because the implementation of a full hierarchical movement model was not computationally feasible. Furthermore, in the Bayesian
setting, Markov Chain Monte Carlo (MCMC) algorithms for most animal movement models require tuning from the user due to lack of
conjugacy. In cases where data sets from tens or hundreds of individuals are available, it may not be feasible to tune individuallevel
MetropolisHastings updates for all parameters.
We present a statistically rigorous twostage procedure for economizing hierarchical animal movement models to provide exact
populationlevel inference using a sequence of algorithms that are fast, stable, and require little or no tuning by the user. Our approach is
simple. First, we fit individuallevel models (1) independently using a preferred stochastic sampling algorithm. Independent model fits in the
first stage allow for parallel processing, leading to an improvement in computational efficiency that scales with the number of processors.
Second, we obtain exact populationlevel inference using a secondary MCMC algorithm that requires no tuning. The secondary algorithm
is based on a littleknown technique for Bayesian metaanalysis proposed by Lunn et al. (2013). We found that our twostage procedure
provides substantial computational improvements in both speed and ease of use in cases with large data sets and/or complicated data models.
In what follows, we present a general twostage procedure for fitting a broad class of hierarchical animal movement models. We then
demonstrate the approach for a basic point process model for telemetry data (i.e., resource selection function model) and verify it using simulation. In our second application, we show how the approach can be applied to a continuoustime discretespace (CTDS) animal movement
model using telemetry data with complicated error structure. We apply the CTDS model to satellite telemetry data from a population of
Canada lynx (Lynx canadensis) in Colorado, USA. Finally, we close with a summary and discussion of the approach and future directions.
2. TWOSTAGE PROCEDURE
Many animal movement models have been constructed solely for individuallevel inference (e.g., Jonsen et al., 2005; Johnson et al., 2008b;
Hooten et al., 2010a; Brost et al., 2015; Buderman et al., 2016). However, the desired scientific inference is usually at the population
level to assess if the population, as a whole, is responding to certain environmental cues. Hierarchical statistical models provide a natural
framework for obtaining upscaled populationlevel inference (Gelman and Hill, 2006; Hobbs and Hooten, 2015). As the complexity of animal
movement models increases, hierarchical models that include nonlinear components become challenging to implement due to computational
limitations and user supervision requirements. It is often much simpler to fit individuallevel models to data, as long as individuals are
assumed independent. Following Lunn et al. (2013), we propose a simple twostage procedure for obtaining populationlevel inference under
the full hierarchical model. The twostage procedure only requires independent individuallevel model fits and an unsupervised resampling
algorithm to obtain populationlevel inference without any user tuning.
The first stage in the procedure involves fitting a data model like (1) independently for each individual j (j D 1; : : : ; J ). In addition to
the prior for auxiliary datalevel parameters � j from (5), we also specify a prior for the individuallevel parameters ˇ j as ˇ j � Œˇ j � (where
the priors for � j and ˇ j can differ by individual). The priors for ˇ j are only used in the first stage of the twostage procedure and do not
affect the final inference. The posterior distribution for individual j is
Œ� j ; ˇ j jyj � D ’
Œyj jˇj ; �j �Œˇj �Œ�j �
Œyj jˇj ; �j �Œˇj �Œ�j �d ˇj d �j
:
(6)
In principle, any stochastic sampling algorithm can be used to obtain samples from the posterior distribution in (6), but those relying on
MCMC are most commonly applied in the animal movement literature. However, because we treat the models in (6) for all J individuals
independently in the first stage, they can be fit in parallel using readily available software (e.g., the ‘parallel’ R package; R Core Team,
2016). Additionally, if we choose a sampling algorithm for fitting the models in (6) that is unsupervised (i.e., requiring no supervised tuning),
then the entire twostage procedure can be automated. An unsupervised fitting procedure will be used much more often by ecologists in
323
Environmetrics 2016; 27: 322–333
Copyright © 2016 John Wiley & Sons, Ltd.
wileyonlinelibrary.com/journal/environmetrics
�Environmetrics
M. B. HOOTEN ET AL.
situations where data exist for a large number of individuals. Thus, automatic MCMC algorithms like BUGS (Lunn et al., 2009), JAGS
(Plummer, 2003), or STAN (Carpenter et al., 2016) can be used to fit the individuallevel models in (6), or alternatively, importance sampling
or particle filtering (e.g., LibBi; Murray, 2013) methods can also be employed. Finally, the choice of priors Œˇ j � can also lead to fully
automatic and parallelizable firststage algorithms. For example, if the data model (1) is Poisson (i.e., yj � Pois.exp.Xj ˇ j //, where Xj is
a design matrix of covariates for the j th individual), then � j is empty because the Poisson does not have a separate dispersion parameter.
A multivariate loggamma prior distribution (Crooks, 2010; Bradley et al., 2015) for ˇ j facilitates the use of a Monte Carlo sampler to
obtain posterior samples from (6). For nonconjugate priors, adaptively tuned MCMC algorithms (e.g., Givens and Hoeting, 2012) are
straightforward to implement and provide a way to obtain unsupervised stageone samples for ˇ j .
The second stage in the twostage procedure involves an MCMC algorithm resembling that used to fit the full hierarchical model but with
a critical simplification. To hfit the ˇfull
i hierarchical model in (1)–(5), we sequentially sample from the fullconditional distributions Œˇ j j�� for
ˇ
�
, using an MCMC algorithm. In our second stage algorithm, we use the MCMC algorithm for the full
j D 1; : : : ; J , Œ�ˇ j��, and † �1
ˇ ˇ
hierarchical model as a template but modify the updates for ˇ j . Updates for the individuallevel auxiliary parameters, � j , are automatically
coupled with those from ˇ j but are only necessary if we desire inference for � j . In fact, if � j are considered nuisance parameters, it is not
necessary to store samples for them in our twostage procedure.
in the second stage model remain the same as in the
The fullconditional distributions for populationlevel parameters �ˇ and † �1
ˇ
MCMC algorithm to fit the full hierarchical model in (1)–(5):
0
1
J
Y
Œˇ j j�ˇ ; † ˇ �A Œ�ˇ � ;
(7)
Œ�ˇ j�� / @
j D1
h
ˇ i
ˇ
† �1
ˇ ˇ�
0
/@
J
Y
1
h
i
:
Œˇ j j�ˇ ; † ˇ �A † �1
ˇ
(8)
j D1
is Wishart, then the fullconditional distributions in
If the model for ˇ j and prior for �ˇ are multivariate Gaussian and the prior for † �1
ˇ
(7) and (8) are multivariate Gaussian and Wishart, respectively. These specific distributions are commonly used in many animal movement
models for populationlevel parameters and permit conjugate Gibbs updates in our secondstage algorithm.
The joint fullconditional distribution for the datalevel auxiliary parameters, � j , and individuallevel parameters, ˇ j , is
Œ� j ; ˇ j j�� / Œyj jˇ j ; � j �Œˇ j j�ˇ ; † ˇ �Œ� j �
(9)
which, depending on the form of data model Œyj jˇ j ; � j �, would normally require a MetropolisHastings update. In this case, the MetropolisHastings ratio for the joint update of � j and ˇ j is
ˇ
h ˇ
ih ˇ
ih ih
i
ˇ
ˇ
ˇ � �
yj ˇ ˇ �j ; � �j ˇ �j ˇ �kˇ ; † kˇ � �j � k�1
; ˇ k�1
ˇ � j ; ˇj
j
j
ˇ
ˇ
ih
ih
ih
i
(10)
rj D h ˇ
ˇ
ˇ k k
ˇ
yj ˇ ˇ k�1
ˇ k�1
� �j ; ˇ �j ˇ � k�1
; � k�1
; ˇ k�1
ˇ �ˇ ; † ˇ � k�1
j
j
j
j
j
j
to the
for the k
where, the ‘�’ superscript represents the proposal for ˇ j and the ‘k’ and ‘k � 1’ superscripts correspond
ˇ MCMC sample
i
h
ˇ
k�1 , is chosen to
;
ˇ
or k � 1 iteration of the MCMC algorithm (for k D 2; : : : ; K). Typically, the proposal distribution, � �j ; ˇ �j ˇ � k�1
j
j
�
��
�
�0
�0
�
�
k�1
k�1
Q
be a multivariate Gaussian random walk such that � j ; ˇ j � N � j ; ˇ j
; † j that requires tuning for each individual j by
Q j using trial and error or an adaptive MCMC approach (e.g., Roberts and Rosenthal, 2009).
adjusting †
However, if we use the posterior samples for � j and ˇ j from the first stage (6) as the proposal in the second stage update for ˇ j , then the
proposal distribution is
h ˇ
i� �h i
ˇ
ˇ
i
h
yj ˇ ˇ �j ; � �j ˇ � � �j
�
� ˇ k�1 k�1
� RR
(11)
� j ; ˇj ˇ � j ; ˇj
Œyj jˇ j ; � j �Œˇ j �Œ� j �d ˇ j d � j
which does not depend on the previous � k�1
and ˇ k�1
. The MetropolisHastings ratio from (10) simplifies to
j
j
ˇ
ih
i
ˇ
ˇ �j ˇ �kˇ ; † kˇ ˇ k�1
j
ˇ
ih i
rj D h
ˇ k k
ˇ k�1
ˇ �ˇ ; † ˇ ˇ �j
j
h
(12)
324
remain unchanged. Thus, we keep the samples for � �j and ˇ �j , from the first stage, with probability
while the updates for �ˇ and † �1
ˇ
min.rj ; 1/. However, we only need to explicitly save samples for the auxiliary individuallevel parameters (� j ) in the first or second
stages if we desire inference on them because rj , from (12), does not depend
ˇ Lunn et
.h
i al. (2013) note that, when the
h ˇ on � j . iFurthermore,
ˇ
ˇ k k
�
ˇ k�1
, a mere quotient involving the
stage one priors for ˇ j are diffuse, the ratio simplifies further to rj D ˇ �j ˇ �kˇ ; † kˇ
;
†
ˇ
j
ˇ
ˇ
wileyonlinelibrary.com/journal/environmetrics
Copyright © 2016 John Wiley & Sons, Ltd.
Environmetrics 2016; 27: 322–333
�POPULATIONLEVEL INFERENCE FOR ANIMAL MOVEMENT
Environmetrics
individuallevel process distributions. However, we retain the form in (12) so that we can use prior information when available. Because
there is no Markov dependence in the proposal for ˇ j , we select ˇ �j (and � �j , if desired) uniformly at random from the output resulting
from the firststage model fits. More importantly, the MetropolisHastings ratios (rj , for j D 1; : : : ; J ) in (12) do not contain a tuning
parameter, resulting in unsupervised updates. Paired with the Gibbs updates for �ˇ and † �1
, the secondstage algorithm is fully automatic,
ˇ
and samples from the fullconditional for ˇ j can be obtained in parallel (within the broader secondstage MCMC algorithm) creating the
potential for additional computational efficiency. Critically, the MetropolisHastings ratio, rj in (12), is not a function of the data. Therefore,
complicated data models do not need to be reconsidered in the secondstage algorithm. The utility of the simple twostage procedure is that
it is intuitive, facilitates parallelization, and can result in algorithms that are fully automatic.
In what follows, we provide two example applications where the twostage procedure for obtaining populationlevel animal movement
inference is valuable. The first application involves a spatial point process modeling approach for telemetry data commonly referred to as
“resource selection function” (RSF) analysis (e.g., Manly et al., 2007). The second application involves a continuoustime discretespace
animal movement model proposed by Hooten et al. (2010a) and Hanks et al. (2015a).
3. APPLICATIONS
3.1. Hierarchical point process models
Perhaps the most common model fit to temporally independent telemetry data is the RSF model. The RSF model is a heterogeneous point
process model that conditions on the number of telemetry observations. Assuming there is no measurement error associated with the telemetry data sij (typically a 2 � 1 vector) for observations i D 1; : : : ; nj and individuals j D 1; : : : ; J , the data model takes the form of a
weighted distribution (Patil and Rao, 1977) such that sij � Œsij jˇ j � and
g.x.sij /; ˇ j /f .sij /
Œsij jˇ j � � R
g.x.s/; ˇ j /f .s/d s
(13)
where g.x.s/; ˇ j / is the “selection” function and f .s/ is the “availability” function. Thus, the animal movement interpretation of (13) is
that inference for ˇ j provides insight about how individual j selects resources (i.e., covariates, x) from those available to it. The selection
function is often chosen to be exponential (i.e., g.x.sij /; ˇ j / � exp.x.sij /0 ˇ j /) and the availability function is typically assumed to be
uniform on the support of the point process (i.e., f .sij / � unif.S/ for sij 2 S � < � <).
Warton and Shepherd (2010) and Aarts et al. (2012) showed that the RSF model in (13) can be fit using a variety of approaches, including
a Poisson likelihood. The Poisson likelihood can be considered by first preprocessing the data such that yj � .y1;j ; : : : ; ym;j /0 represents
counts of telemetry locations in grid cells corresponding to a discretization of the support S. As the grid cell size decreases with respect to
the resolution of the covariates x, a Poisson data model coincides with the point process model. Thus, the corresponding hierarchical model
yj � Pois.exp.Xj ˇ j //
(14)
ˇ j � N.�ˇ ; † ˇ /
(15)
�ˇ � N.�0 ; † 0 /
(16)
�1
† �1
ˇ � Wish..S�/ ; �/
(17)
assumes the same form as (1)–(5) and allows for populationlevel resource selection inference on �ˇ . To fit the full hierarchical model
, sequentially. Standard MetropolisHastings
directly using MCMC, we sample from the fullconditional distributions for ˇ j , �ˇ , and † �1
ˇ
updates for ˇ j require tuning, but the model can be fit using a single MCMC algorithm for moderately sized data sets. Alternatively, the
weighted leastsquares proposal approach of Gamerman (1997) could be used to acquire samples for ˇ j from the posterior distribution.
However, to adequately approximate the point process model, the grid cells often need to be quite small, resulting in a finescale discretization
of the support S and increasing the computational burden.
The twostage procedure we described in the previous section can easily be employed to fit the hierarchical model in (14)–(17). For the
first stage, we can use an MCMC or Hamiltonian Monte Carlo algorithm (via BUGS, JAGS, or STAN; Lunn et al., 2009; Plummer, 2003;
Carpenter et al., 2016) to fit the individual level models in parallel. For our spatial point process setting, the individuallevel models are
yj � Pois.exp.Xj ˇ j //
(18)
ˇ j � N.�0 ; † 0 /
(19)
for j D 1; : : : ; J , independently. Note that the individuallevel parameter model in (19) is an exchangeable prior for all j D 1; : : : ; J . Also,
if the individual data sets yj and Xj are so large that they are difficult to store in memory simultaneously for all J individuals, the firststage
model fitting can be fully distributed among separate machines or performed in sequence. This highlights another primary advantage of the
twostage procedure.
325
Environmetrics 2016; 27: 322–333
Copyright © 2016 John Wiley & Sons, Ltd.
wileyonlinelibrary.com/journal/environmetrics
�Environmetrics
M. B. HOOTEN ET AL.
Figure 1. (a) Simulated animal positions (points) based on a spatial point process (13) with one simulated covariate (background image, dark shading
represents larger values). (b) A zoomed in spatial map (from inset white box in (a)) showing positions from five individual animals as different point types
(i.e., ˘, 4, 5, C, �)
The secondstage algorithm for obtaining populationlevel inference is an MCMC algorithm with Gibbs updates for �ˇ and † �1
as
ˇ
described in the previous section and updates for ˇ j using MetropolisHastings based on the acceptance ratio in (12), which becomes
ˇ
� ˇ
� �
�
ˇ
ˇ
N ˇ �j ˇ �kˇ ; † kˇ N ˇ k�1
ˇ �0 ; † 0
j
ˇ
� � ˇ
� :
rj D �
ˇ k k
ˇ
�
N ˇ �j ˇ �0 ; † 0
N ˇ k�1
;
†
ˇ
j
ˇ
ˇ
(20)
Within the second stage MCMC algorithm, the updates for ˇ j can also be parallelized because they are independent, although this model
is simple enough that parallelization is not necessary in the second stage algorithm. Thus, the data, yj for j D 1; : : : ; J , which could include
counts for 10s or 100s of thousands of grid cells and 100s of individuals, do not appear in the second stage algorithm. The absence of yj
leads to a more computationally efficient second stage algorithm than the original algorithm to fit the full hierarchical model directly.
We simulated point process data from 20 individuals (Figure 1), resulting in approximately 30 simulated telemetry fixes per individual,
and fit the hierarchical RSF model using: 1.) a single MCMC algorithm, and 2.) our twostage procedure. We compared the populationlevel
results from the fits resulting from each procedure.
For the firststage algorithm in our twostage procedure, we fit the individuallevel models independently using an adaptive MCMC
algorithm in parallel using R (R Core Team, 2016) and assumed N.0; 100�I/ priors for ˇ j , a N.0; 100�I/ prior for �ˇ , and a Wish..3�I/�1 ; 3/
. Our firststage algorithm uses a multivariate Gaussian proposal for ˇ j and adapts the tuning using a single variance parameter,
prior for † �1
ˇ
resulting in an unsupervised algorithm for the individuallevel model fits. We could have also used BUGS or JAGS to fit the firststage
models, but our adaptive MCMC algorithm required less computing time.
The single MCMC algorithm to fit the full hierarchical model required 2.62 min to obtain 20,000 MCMC samples in R, whereas the
firststage algorithm required 0.57 min to obtain the same number of samples using an adaptive MCMC algorithm in parallel for the 20
individuals. The secondstage algorithm required only 1.49 min in R, which implies that the total compute time to fit the model using the
twostage procedure was 2.06 min (0.56 min less than the single MCMC algorithm). Also, the twostage procedure requires no tuning and
results in much larger effective MCMC sample sizes for parameters. The effective MCMC sample sizes for �ˇ and ˇ j were 8560 and 1398
(averaged across individuals) for the single MCMC algorithm, but were 17,590 and 15,184 for the twostage algorithm (out of 20,000 total
samples). Thus, to obtain the same effective MCMC sample size using MCMC for all parameters, we would need an order of magnitude
more samples from the single MCMC algorithm.
Figure 2 illustrates the similarities in inference for the slope parameters �ˇ1 and ˇj1 for j D 1; : : : ; 20 when fitting the hierarchical RSF
model using a single MCMC algorithm (black) versus the twostage procedure (gray). Notice that the single MCMC algorithm and the twostage procedure provide very similar inference. In terms of inference, there exists some variability among individuals, but the populationlevel
inference (Figure 2, top) suggests a consistent overall positive population response to the covariate.
3.2. Hierarchical continuoustime discretespace models
The previous application, involving spatial point process models, involves a commonly used model specification and desired type of inference in ecological research, but more contemporary methods have been developed to explicitly model the dynamics of animal movement
based on temporally dependent telemetry data with observations close in time. Among these methods are discretetime and continuoustime
approaches to modeling the individual animal trajectories (McClintock et al., 2014). We focus on the continuoustime class of models in
what follows.
326
wileyonlinelibrary.com/journal/environmetrics
Copyright © 2016 John Wiley & Sons, Ltd.
Environmetrics 2016; 27: 322–333
�POPULATIONLEVEL INFERENCE FOR ANIMAL MOVEMENT
Environmetrics
Figure 2. Posterior means (points) and 50% and 95% credible intervals for �ˇ1 and ˇj1 for j D 1; : : : ; 20. Single MCMC algorithm results are shown in
black, and twostage procedure results are shown in gray
Continuoustime statistical models for animal movement processes have existed for decades (e.g., Dunn and Gipson, 1977; Blackwell,
1997) and are usually based on Brownian motion (i.e., Wiener processes). Up until the late 1990s, most Brownian motion models for
trajectories utilized an Ornstein–Uhlenbeck process (i.e., a Wiener process with attraction to a central position). Johnson et al. (2008b) also
proposed an Ornstein–Uhlenbeck model but for the velocity (i.e., temporally differentiated position) rather than the position process. Hooten
and Johnson (2016a) generalized the continuoustime velocity models of Johnson et al. (2008b) in the context of Gaussian processes with
covariance structure induced by temporal basis functions. Buderman et al. (2016) used a simplified basis function parameterization to model
Canada lynx (Lynx canadensis) movement while accounting for measurement error in the telemetry data. Buderman et al. (2016) refer to
their model as a “functional movement model” and use it to provide inference for the true underlying continuous position process (i.e., �.t /,
for time t ) of an individual.
327
Environmetrics 2016; 27: 322–333
Copyright © 2016 John Wiley & Sons, Ltd.
wileyonlinelibrary.com/journal/environmetrics
�Environmetrics
M. B. HOOTEN ET AL.
The approach developed by Buderman et al. (2016) assumes that the telemetry data sij are observed with error. In fact, for the Canada lynx
in our study, the bivariate measurement error follows an unusual Xshaped pattern because the telemetry data are collected by Service Argos
(Costa et al., 2010) that relies on polar orbiting satellites. Thus, Brost et al. (2015) and Buderman et al. (2016) developed a measurement
error model based on a mixture distribution to account for the Xshaped Argos pattern (see Appendix A for details). Properly accounting for
measurement error adds another level to the hierarchical model in (1)–(5) such that
sij � Œsij j�j .ti /; �j �
(21)
yj � Œyj jˇ j ; � j �
(22)
ˇ j � Œˇ j j�ˇ ; † ˇ �
(23)
�ˇ � Œ�ˇ �
(24)
i
h
�1
† �1
ˇ � †ˇ
(25)
� j � Œ� j �
(26)
�j � Œ�j �
(27)
for j D 1; : : : ; J individuals, and where yj is an mj � 1 vector that represents a latent process that is linked to the true continuous position
process f�j .t /; 8t g by a deterministic functional h such that yj D h.f�j .t /; 8t g/ and �j are measurement error covariance parameters.
Hooten et al. (2010a) developed an individuallevel animal movement model based on (21) and (22) where the latent variables yj represent
a sequential multinomial process indicating transitions among grid cells on a discretization of the spatial support S. The latent process model
in (22) relies on a continuoustime discretespace (CTDS) representation of the position process. However, because the functional h.�/, which
links the position process with the data, is noninvertible in their model, Hooten et al. (2010a) proposed a Bayesian multiple imputation
procedure to account for uncertainty in the true position process when making inference on ˇ j . The multiple imputation procedure used by
Hooten et al. (2010a) differs from the twostage procedure we described herein because it does not allow for feedback from the individuallevel parameters ˇ j to the position process f�j .t /; 8t g or measurement error parameters �j . Hooten et al. (2010a) used an imputation
model to interpolate the position process and then integrated over the uncertainty in the position process while fitting (22) to provide posterior
inference for the individuallevel parameters ˇ j .
Hanks et al. (2015a) showed that the multinomial process of Hooten et al. (2010a) could be reparameterized such that Œyj jˇ j ; � j � can
be modeled using Poisson regression. Specifically, let �cj represent the amount of time individual j remains in a grid cell for the cth
“stay/move”
pair
with the discretization of the individual’s path through a landscape (for c D 1; : : : ; nj ). Then let yclj �
��
�
� associated
where the index l D 1; 2; 4; 5 (l D 3 is not necessary because it corresponds to the middle cell that is captured by
Pois �cj exp x0clj ˇ j
�cj ) denotes moves to neighboring grid cells in each cardinal direction (i.e., north, east, south, west). That is, if individual j moved north for
“stay/move” pair c, then the data point yc1j D 1 and yc2j D yc3j D yc4j D 0 (see Appendix B for details). The Poisson reparametrization
dramatically improves computational efficiency at the individual level because the total number of observations used in the model (4mj ) is a
function of the grid cell size rather than the position process discretization as used in Hooten et al. (2010a). Thus, Hanks et al. (2015a) were
able to fit the CTDS model to large telemetry data sets in a fraction of the time required by the multinomial method developed by Hooten et
al. (2010a). However, neither Hooten et al. (2010a) nor Hanks et al. (2015a) attempted to fit a hierarchical model like that in (21)–(27) to
obtain population level inference for �ˇ .
In our application involving populationlevel inference for Canada lynx,˚ we use the� model developed by Buderman et al. (2016) to obtain
Q j .t /; 8t and hence yQclj for all c, l, and j , while accounting
the imputation distribution for the true individuallevel position process �
for the complicated nature of Argos telemetry error (see Appendix A for details). In what follows, we combine all yQclj into a single vector
representing the latent process yQ j and use yQ j as data in a twostage implementation of the hierarchical model in (21)–(27).
To fit the hierarchical model using the twostage procedure described in Section 2, we apply the same two stages of algorithms as in
the previous application. For the first stage, we use the data model in (22) and specify multivariate Gaussian priors for the individuallevel
parameters ˇ j � N.�0 ; † 0 /. We use an adaptively tuned MCMC algorithm to obtain samples from the posterior distributions
Z
Œˇ j jfsij ; 8i; j g� D
��
�
�
ˇ j jQyj yQ j jfsij ; 8i; j g d yQ j
(28)
�
�
for j D 1; : : : ; J , and where yQ j jfsij ; 8i; j g represents the imputation distribution for the latent Poisson process. To perform the integration
�
�
in (28), we sample yQ kj � yQ j jfsij ; 8i; j g on the kth MCMC iteration and then let the MetropolisHastings update ˇ kj depend on yQ kj as
described in Hooten et al. (2010a) and Hanks et al. (2015a). As in the first application, we can fit the J models for all individuals in parallel,
dramatically reducing the required computational time.
328
wileyonlinelibrary.com/journal/environmetrics
Copyright © 2016 John Wiley & Sons, Ltd.
Environmetrics 2016; 27: 322–333
�POPULATIONLEVEL INFERENCE FOR ANIMAL MOVEMENT
Environmetrics
Figure 3. (a) Colorado, USA, with major highways and the city of Denver shown. The telemetry data spanning a year of time for 18 individual Canada lynx
are shown as points. A shaded relief map is shown as the background image to illustrate the topography of the area. (b) and (c) closeup views of the predicted
paths for two individual Canada lynx. For clarity, only the posterior mean path is shown
For the second stage of the twostage procedure, we use the posterior samples for fˇ j ; 8j g, from the first stage, as proposals in the
sequentially in a completely
MCMC algorithm to fit the hierarchical model in (22)–(25). In doing so, we update fˇ j ; 8j g, �ˇ , and † �1
ˇ
unsupervised secondstage MCMC algorithm. Recall that the MetropolisHastings acceptance ratio for ˇ j is identical to that used in the
previous application (20). As a result of the twostage implementation and the adaptive tuning in the firststage algorithm, the procedure is
completely automatic after the data are preprocessed to obtain the imputation distribution, and populationlevel inference for �ˇ can easily
be obtained.
Using telemetry data from J D 18 individual Canada lynx in Colorado, USA (Figure 3(a)), we applied the twostage procedure to fit the
hierarchical model in (22)–(25). We used the functional movement model of Buderman et al. (2016) to obtain the imputed path distribution
(Figure 3(b,c)) for each individual and used nearly continuous imputed path realizations to create the latent Poisson data realizations yQ kj
(resulting in approximately 450 discretespace transitions per individual, nj
450). Canada lynx are a subalpine species that tend to prefer
forested ecosystems (McKelvey et al., 2000); thus, we focused on two covariates: elevation and distance to forest (Figure 4). Each covariate
was included in the model as a “static” driver, rather than a gradientbased driver of movement (Hanks et al., 2015a). Static drivers can be
interpreted as affecting overall motility in the CTDS model. For priors in the first stage, we used ˇ j � N.0; 100I/ for all j D 1; : : : ; 18. We
� Wish..3 � I/�1 ; 3/ as priors for the populationlevel parameters and precision matrix. See Appendix B
used �ˇ � N.0; 100I/ and † �1
ˇ
for additional details on the CTDS animal movement model.
329
Environmetrics 2016; 27: 322–333
Copyright © 2016 John Wiley & Sons, Ltd.
wileyonlinelibrary.com/journal/environmetrics
�Environmetrics
M. B. HOOTEN ET AL.
Figure 4. Images of covariates with telemetry observations overlaid as black points: (a) Elevation and (b) distance to forest. Light shading corresponds to
larger values
Figure 5. Posterior estimates for the populationlevel parameters �ˇ and individuallevel parameters ˇ j . The posterior mean is shown as a central point,
and the 50% and 95% credible intervals are shown as the thick and thin black lines. (a) Results for the elevation covariate and (b) results for the distance from
forest covariate
We fit the overall hierarchical model using the twostage procedure, and the resulting algorithms required 0.86 min for the first stage
(using an adaptive MCMC algorithm in parallel) and 1.62 min for the second stage. Figure 5 shows the results of the model fit in terms of
posterior means, and 50% and 95% credible intervals for the populationlevel parameters �ˇ and individuallevel parameters ˇ j .
While there exists substantial variability among individual Canada lynx, with some individuals exhibiting clear relationships with the
covariates (e.g., individuals 2, 4, and 5), the posterior distributions for � did not indicate a populationlevel effect for either covariate at the
95% level (but both did at the 50% level). For the individuals that did show evidence of an effect (i.e., 95% credible intervals not overlapping
zero), the negative response to elevation indicates that overall motility decreases at higher elevations, leading to greater residence times in
those regions, as opposed to lower elevations (Figure 5(a)). Similarly, for individuals with significant effects related to distance from forest
we see positive influence on motility implying that those Canada lynx have higher motility (and hence lower residence time) in regions
farther from forest (Figure 5(b)). Thus, the inference in our application involving Canada lynx agrees with that obtained in other studies
(e.g., McKelvey et al., 2000).
4. CONCLUSION
Our findings indicate that the twostage procedure we described herein holds tremendous value for fitting hierarchical animal movement models to telemetry data for populationlevel inference. We applied the twostage procedure to two types of commonly used animal movement
models of varying complexity and found that it worked well in both cases.
The spatial point process modeling approach we described in the first application is a commonly used model, but still fairly simple.
Much more complicated spatiotemporal point process models have been used to model temporally correlated telemetry data (e.g., Johnson
et al., 2008a; Johnson et al., 2013; Brost et al., 2015), and adapting the twostage procedure to those models is the subject of ongoing
research. For example, Brost et al. (2015) developed a model with a timevarying dynamic availability component that depended on an
additional smoothness parameter. Thus, the data model developed by Brost et al. (2015) required substantially more computation time than
330
wileyonlinelibrary.com/journal/environmetrics
Copyright © 2016 John Wiley & Sons, Ltd.
Environmetrics 2016; 27: 322–333
�POPULATIONLEVEL INFERENCE FOR ANIMAL MOVEMENT
Environmetrics
the simulated example we presented in Section 3.1 and would benefit from a twostage implementation where individuallevel models could
be fit independently on separate processors and then recombined using the secondstage MCMC algorithm to yield populationlevel inference
for �ˇ .
In our example involving Canada lynx, the continuoustime discretespace reparameterization developed by Hanks et al. (2015a) already
provides significant improvements in computational efficiency over the motivating model developed by Hooten et al. (2010a). However,
additional computational gains can be achieved using the twostage fitting procedure to provide populationlevel inference.
Despite the wide range of potential applications to many types of hierarchical models, we found it surprising that the twostage fitting
procedure of Lunn et al. (2013) is not more well known. For our situations with large amounts of telemetry data and potentially complicated
data models, we found the twostage procedure works very well and is trivial to implement. We also found it very helpful to be able to use
different data models, firststage fitting algorithms, and easy parallelization. As a potential caveat, the twostage procedure described by Lunn
et al. (2013) may not be very efficient when the population induces extreme amounts of shrinkage in the individuallevel parameters. Thus,
in these cases, more samples would be needed in the firststage algorithm. However, in a preliminary simulation study, we found that the
twostage procedure performs poorly only for data sets with very small amounts of data (i.e., < 20 observations for a subset of individuals).
Animal movement models have also been developed to account for more mechanistic interactions among individuals (e.g., Russell et
al., 2016; Scharf et al., 2015), and while we did not address those specifically, the approach we presented may also be beneficial in those
settings. Furthermore, Bayesian animal movement models have been fit using integrated nested Laplace approximation (INLA; Rue et al.,
2009; Illian et al., 2012; Illian et al., 2013; RuizCárdenas et al., 2012; Jonsen, 2016), and one could use INLA to fit the hierarchical point
process model in our first example. However, the twostage MCMC approach presented herein allows for the following: inference on joint
relationships among model parameters, easy parallelization in the first stage, and the ability to use Bayesian multiple imputation techniques,
such as in our second example involving the CTDS movement model.
Acknowledgements
Support for this research was provided by NSF 1614392, NSF EEID 1414296, CPW T01304, and NOAA AKC 188000. The authors thank
The International Environmetrics Society and the editors of Environmetrics for their support and assistance with this work. The authors
also thank Walt Piegorsch, Mindy Rice, Devin Johnson, Peter Craigmile, Erin Peterson, Ron Smith, and the other organizers of the TIES
2016 annual meeting. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the
U.S. Government.
REFERENCES
Aarts G, Fieberg J, Matthiopoulos J. 2012. Comparative interpretation of count, presenceabsence, and point methods for species distribution models.
Methods in Ecology and Evolution 3:177–187.
Berliner L. 1996. Hierarchical Bayesian time series models. In Maximum Entropy and Bayesian Methods, K Hanson, R Silver (eds). Kluwer Academic
Publishers: Dordrecht, the Netherlands, 15–22.
Blackwell P. 1997. Random diffusion models for animal movement. Ecological Modelling 100:87–102.
Bradley J, Holan S, Wikle C. 2015. Computationally efficient distribution theory for Bayesian inference of highdimensional dependent countvalued data.
arXiv preprint arXiv:1512.07273.
Brost B, Hooten M, Hanks E, Small R. 2015. Animal movement constraints improve resource selection inference in the presence of telemetry error. Ecology
96:2590–2597.
Buderman F, Hooten M, Ivan J, Shenk T. 2016. A functional model for characterizing long distance movement behavior. Methods in Ecology and Evolution
7:264–273.
Cagnacci F, Boitani L, Powell RA, Boyce MS. 2010. Animal ecology meets gpsbased radiotelemetry: a perfect storm of opportunities and challenges.
Philosophical Transactions of the Royal Society of London B: Biological Sciences 365:2157–2162.
Carpenter B, Gelman A, Hoffman M, Lee D, Goodrich B, Betancourt M, Brubaker M, Guo J, Li P, Riddell A. 2016. Stan: a probabilistic programming
language. Journal of Statistical Software.
Costa D, Robinson P, Arnould J, Harrison AL, Simmons SE, Hassrick JL, Hoskins AJ, Kirkman SP, Oosthuizen H, VillegasAmtmann S, Crocker DE. 2010.
Accuracy of Argos locations of pinnipeds atsea estimated using Fastloc GPS. PLoS One 5:e8677, 1–9.
Cressie N, Wikle C. 2011. Statistics for SpatioTemporal Data. John Wiley and Sons: New York, New York, USA.
Crooks G. 2010. The amoroso distribution. arXiv preprint arXiv:1005.3274.
Dunn J, Gipson P. 1977. Analysis of radiotelemetry data in studies of home range. Biometrics 33:85–101.
Gamerman D. 1997. Sampling from the posterior distribution in generalized linear mixed models. Statistics and Computing 7:57–68.
Gelfand A, Smith A. 1990. Samplingbased approaches to calculating marginal densities. Journal of the American Statistical Association 85:398–409.
Gelman A, Hill J. 2006. Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press: Cambridge, United Kingdom.
Givens G, Hoeting J. 2012. Computational Statistics, Vol. 710. John Wiley & Sons: Hoboken, New Jersey.
Hanks E, Hooten M, Alldredge M. 2015a. Continuoustime discretespace models for animal movement. Annals of Applied Statistics 9:145–165.
Hanks E, Hooten M, Johnson D, Sterling J. 2011. Velocitybased movement modeling for individual and population level inference. PLoS One 6:e22795,
1–17.
Hobbs N, Hooten M. 2015. Bayesian Models: A Statistical Primer for Ecologists. Princeton University Press: Princeton, New Jersey, USA.
Hooten M, Johnson D. 2016a. Basis function models for animal movement. Journal of the American Statistical Association. In Revision.
Hooten M, Johnson D, Hanks E, Lowry J. 2010a. Agentbased inference for animal movement and selection. Journal of Agricultural, Biological and
Environmental Statistics 15:523–538.
Hooten M, Johnson D, McClintock B, Morales J. in press. Animal Movement: Statistical Models for Telemetry Data. Chapman & Hall/CRC.
Illian J, Sorbye S, Rue H, Hendrichsen D. 2012. Using INLA to fit a complex point process model with temporally varying effects – a case study. Journal of
Environmental Statistics 3:1–25.
Illian J, Martino S, Sørbye S, GallegoFernández J, Zunzunegui M, Esquivias M, Travis J. 2013. Fitting complex ecological point process models with
integrated nested Laplace approximation. Methods in Ecology and Evolution 4:305–315.
331
Environmetrics 2016; 27: 322–333
Copyright © 2016 John Wiley & Sons, Ltd.
wileyonlinelibrary.com/journal/environmetrics
�Environmetrics
M. B. HOOTEN ET AL.
Johnson D, Hooten M, Kuhn C. 2013. Estimating animal resource selection from telemetry data using point process models. Journal of Animal Ecology
82:1155–1164.
Johnson D, London J, Lea M, Durban J. 2008b. Continuoustime correlated random walk model for animal telemetry data. Ecology 89:1208–1215.
Johnson D, Thomas D, Ver Hoef J, Christ A. 2008a. A general framework for the analysis of animal resource selection from telemetry data. Biometrics
64:968–976.
Jonsen I. 2016. Joint estimation over multiple individuals improves behavioural state inference from animal movement data. Scientific Reports 6:20625, 19.
Jonsen I, Flemming J, Myers R. 2005. Robust statespace modeling of animal movement data. Ecology 45:589–598.
Lunn D, Barrett J, Sweeting M, Thompson S. 2013. Fully Bayesian hierarchical modelling in two stages, with application to metaanalysis. Journal of the
Royal Statistical Society: Series C (Applied Statistics) 62:551–572.
Lunn D, Spiegelhalter D, Thomas A, Best N. 2009. The BUGS project: evolution, critique and future directions. Statistics in Medicine 28:3049–3067.
Manly B, McDonald L, Thomas D, McDonald T, Erickson W. 2007. Resource Selection by Animals: Statistical Design and Analysis for Field Studies.
Springer Science & Business Media: Dordrecht, the Netherlands.
McClintock B, Johnson D, Hooten M, Ver Hoef J, Morales J. 2014. When to be discrete: the importance of time formulation in understanding animal
movement. Movement Ecology 2:1–14.
McKelvey K, Aubry K, Ortega Y. 2000. History and distribution of lynx in the contiguous United States. In Ecology and Conservation of Lynx in the United
States, L Ruggiero, J Squires, S Buskirk, K Aubry, K McKelvey, G Koehler, C Krebs (eds). University Press of Colorado: Boulder, Colorado, USA,
207–264.
Murray L. 2013. Bayesian statespace modelling on highperformance hardware using libBi. arXiv preprint arXiv:1306.3277.
Patil G, Rao C. 1977. The weighted distributions: a survey of their applications. In Applications of Statistics, P Krishnaiah (ed.). North Holland Publishing
Company: Amsterdam, the Netherlands.
Plummer M. 2003. JAGS: a program for analysis of Bayesian graphical models using Gibbs sampling. Proceedings of the 3rd International Workshop on
Distributed Statistical Computing, Vol. 124. Technische Universit at Wien Wien, Austria, 125.
R Core Team. 2016. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing: Vienna, Austria.
Roberts G, Rosenthal J. 2009. Examples of adaptive MCMC. Journal of Computational and Graphical Statistics 18:349–367.
Rue H, Martino S, Chopin N. 2009. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal
of the Royal Statistical Society: Series B (Statistical Methodology) 71(2):319–392.
RuizCárdenas R, Krainski E, Rue H. 2012. Direct fitting of dynamic models using integrated nested Laplace approximations INLA. Computational Statistics
& Data Analysis 56:1808–1828.
Russell JC, Hanks EM, Haran M, Hughes DP. 2016. A spatiallyvarying stochastic differential equation model for animal movement. arXiv preprint
arXiv:1603.07630.
Scharf H, Hooten M, Fosdick B, Johnson D, London J, Durban J. 2015. Dynamic social networks based on movement. arXiv:1512.07607.
Warton D, Shepherd L. 2010. Poisson point process models solve the “pseudoabsence problem” for presenceonly data in ecology. Annals of Applied
Statistics 4:1383–1402.
APPENDIX A. THE IMPUTATION DISTRIBUTION
Buderman et al. (2016) developed a phenomenological statistical model for estimating an individual’s underlying continuoustime path based
on Argos telemetry data and a semiparametric regression using temporal basis functions. We used this model to precalculate an imputation
distribution for the true path. For the j th individual, the FMM developed by Buderman et al. (2016) is
n
o
Figure A1. (a) The telemetry data (large points) and imputation distribution (lines) for the individual’s path �k
j .t /; 8k; j; t using the posterior predictive
distribution of the FMM (Buderman et al., 2016). (b) A closeup view of the imputation distribution showing the temporal discretization of the imputed path
realizations. (c) A closeup view of a single imputed path realization crossing through the firstorder neighborhood of the center grid cell
332
wileyonlinelibrary.com/journal/environmetrics
Copyright © 2016 John Wiley & Sons, Ltd.
Environmetrics 2016; 27: 322–333
�POPULATIONLEVEL INFERENCE FOR ANIMAL MOVEMENT
sij �
N.�j .ti /; † i /
N.�j .ti /; H† i H0 /
with prob. p
with prob. 1 � p
Environmetrics
(A1)
�j .ti / D Wj .ti /˛
(A2)
˛ � N.0; † ˛ /
(A3)
where sij represents the i th telemetry observation, �j .ti / is the true individual position at time ti , † i is an error covariance matrix on the
first axis, and H† i H0 is the error covariance matrix on a rotated axis (H is a rotation matrix). The probability p allows the telemetry data to
arise from a bivariate Gaussian mixture that captures the Xshaped error pattern inherent to Argos data. The matrix Wj .ti / contains basis
vectors (i.e., bspline basis vectors) at time ti for individual j , and ˛ is a set of regression coefficients corresponding to the temporal basis
functions. Buderman et al. (2016) set † ˛ � Diag � 2˛ and tuned � 2˛ to induce regularization in the model and improve predictive ability
(i.e., ridge regression).
The imputed path distribution is obtained by sampling from the posterior predictive distribution of Œ�j .t /jfsij ; 8i; j g� for a large, but
finite, set of times t 2 T to obtain posterior realizations �kj .t / for k D 1; : : : ; K MCMC iterations. Figure A1(a)shows an example set of
path realizations (lines) that could result from fitting the FMM from Buderman et al. (2016) to telemetry data (points).
Figure A1(b) shows a zoomed in section of the path realizations that highlight the temporal discretization. At a finer spatial resolution, we
can see that the path realizations cross through an example grid cell and its associated neighborhood (Figure A1(c)). This idea is critical for
processing the path realizations for use with the CTDS model.
APPENDIX B. CTDS MODEL
For each individual j in the original CTDS model, each segment (between points) in Figure A1(c) served as a multinomial data vector
yij � .y1i ; y2i ; y3i ; y4i ; y5i /0j , where yij � MN.1; pij / (Hooten et al., 2010a). The multinomial vectors were constructed using the
function yij D h.f�j .t /; 8t g/ based on the imputed path realizations by coding a transition as either a stay or a move in a certain direction
according to the schematic in Figure B1.
Figure B1. Discrete set of possible transitions at any time t , used to create the multinomial vector y.t /, based on the function h.f�j .t /; 8t g/. (a) Move
up: y.t / D .1; 0; 0; 0; 0/0 , (b) move right: y.t / D .0; 1; 0; 0; 0/0 , (c) stay: y.t / D .0; 0; 1; 0; 0/0 , (d) move down: y.t / D .0; 0; 0; 1; 0/0 , and (e) move
left: y.t / D .0; 0; 0; 0; 1/0
Hanks et al. (2015a) reparameterized the multinomial imputation data using sufficient statistics. They denoted residence time as �lj
(approximated by �t times the number of consecutive stays in the current grid cell, Figure A1(c)) for l D 1; : : : ; L “stay–move” pairs
�lj =�t
and then defined the probability of staying in the current grid cell for time �lj as p3ij
D .1 � plj;move /�lj =�t , where plj;move is the
probability of moving. Hanks et al. (2015a) let plj;move D �t � �lj;move and �t ! 0 yielding
lim .1 � pj;move /�lj =�t D e ��lj �lj;move
(B1)
�t!0
which implies that �lj � Exp.�lj;move /.
Similarly, Hanks et al. (2015a) showed that the movement probability to neighboring grid cell c is pclj =plj;move D �clj =�lj;move .
Thus, combing the residence probability
Qmodel with the movement probability yields a likelihood for the sufficient statistic
Q
.�lj ; y1lj ; y2lj ; y4lj ; y5lj /0 equal to L
lD1 c¤3 �clj exp.��lj �clj /. The likelihood for this reparameterized CTDS model coincides with
a Poisson where �clj is the movement rate to neighboring cell c and �lj is an offset. Thus, any software capable of fitting a Poisson
generalized linear model with an offset can fit the CTDS model if the true path is observed at a fine enough temporal resolution.
Hanks et al. (2015a) used a multiple imputation approach to account for the uncertainty in the path distribution based on (28). The
movement rates can then be linked to the environmental covariates by a loglinear link �clj D x0clj ˇ j , where the covariates x0clj can
be specified in several meaningful ways to capture either differential movement rates (i.e., motility) or gradientbased directional bias in
movement relative to environmental covariates (see Hanks et al., 2015a for details). The reparameterized CTDS model of Hanks et al.
(2015a) is much more computationally efficient than that of Hooten et al. (2010a) because the dimensionality of the data 4L depends on the
grid cell size instead of the temporal discretization of the path.
333
Environmetrics 2016; 27: 322–333
Copyright © 2016 John Wiley & Sons, Ltd.
wileyonlinelibrary.com/journal/environmetrics
�
Dublin Core
The Dublin Core metadata element set is common to all Omeka records, including items, files, and collections. For more information see, http://dublincore.org/documents/dces/.
Title
A name given to the resource
Journal Articles
Description
An account of the resource
CPW peerreviewed journal publications
Text
A resource consisting primarily of words for reading. Examples include books, letters, dissertations, poems, newspapers, articles, archives of mailing lists. Note that facsimiles or images of texts are still of the genre Text.
Dublin Core
The Dublin Core metadata element set is common to all Omeka records, including items, files, and collections. For more information see, http://dublincore.org/documents/dces/.
Title
A name given to the resource
Hierarchical animal movement models for populationlevel inference
Description
An account of the resource
<span>New methods for modeling animal movement based on telemetry data are developed regularly. With advances in telemetry capabilities, animal movement models are becoming increasingly sophisticated. Despite a need for populationlevel inference, animal movement models are still predominantly developed for individuallevel inference. Most efforts to upscale the inference to the population level are either </span><i>post hoc</i><span> or complicated enough that only the developer can implement the model. Hierarchical Bayesian models provide an ideal platform for the development of populationlevel animal movement models but can be challenging to fit due to computational limitations or extensive tuning required. We propose a twostage procedure for fitting hierarchical animal movement models to telemetry data. The twostage approach is statistically rigorous and allows one to fit individuallevel movement models separately, then resample them using a secondary MCMC algorithm. The primary advantages of the twostage approach are that the first stage is easily parallelizable and the second stage is completely unsupervised, allowing for an automated fitting procedure in many cases. We demonstrate the twostage procedure with two applications of animal movement models. The first application involves a spatial point process approach to modeling telemetry data, and the second involves a more complicated continuoustime discretespace animal movement model. We fit these models to simulated data and real telemetry data arising from a population of monitored Canada lynx in Colorado, USA.</span>
Bibliographic Citation
A bibliographic reference for the resource. Recommended practice is to include sufficient bibliographic detail to identify the resource as unambiguously as possible.
Hooten, M. B., F. E. Buderman, B. M. Brost, E. M. Hanks, and J. S. Ivan. 2016. Hierarchical animal movement models for populationlevel inference. Environmetrics 27:322333. <a href="https://doi.org/10.1002/env.2402" target="_blank" rel="noreferrer noopener">https://doi.org/10.1002/env.2402</a>
Creator
An entity primarily responsible for making the resource
Hooten, Mevin B.
Buderman, Frances E.
Hanks, Ephraim M.
Ivan, Jacob S.
Subject
The topic of the resource
Hierarchical model
Resource selection model
Spatial statistics
Telemetry data
Trajectories
Extent
The size or duration of the resource.
12 pages
Date Created
Date of creation of the resource.
20160719
Rights
Information about rights held in and over the resource
<a href="http://rightsstatements.org/vocab/InCNC/1.0/" target="_blank" rel="noreferrer noopener">In Copyright  NonCommercial Use Permitted</a>
Format
The file format, physical medium, or dimensions of the resource
application/pdf
Language
A language of the resource
English
Is Part Of
A related resource in which the described resource is physically or logically included.
Environmetrics
Type
The nature or genre of the resource
Article