 Research
 Open Access
 Published:
Identifying resting locations of a small elusive forest carnivore using a twostage model accounting for GPS measurement error and hidden behavioral states
Movement Ecology volume 9, Article number: 17 (2021)
Abstract
Background
Studies of animal movement using location data are often faced with two challenges. First, time series of animal locations are likely to arise from multiple behavioral states (e.g., directed movement, resting) that cannot be observed directly. Second, location data can be affected by measurement error, including failed location fixes. Simultaneously addressing both problems in a single statistical model is analytically and computationally challenging. To both separate behavioral states and account for measurement error, we used a twostage modeling approach to identify resting locations of fishers (Pekania pennanti) based on GPS and accelerometer data.
Methods
We developed a twostage modelling approach to estimate when and where GPScollared fishers were resting for 21 separate collar deployments on 9 individuals in southern Oregon. For each deployment, we first fit independent hidden Markov models (HMMs) to the time series of accelerometerderived activity measurements and apparent step lengths to identify periods of movement and resting. Treating the state assignments as given, we next fit a set of linear Gaussian state space models (SSMs) to estimate the location of each resting event.
Results
Parameter estimates were similar across collar deployments. The HMMs successfully identified periods of resting and movement with posterior state assignment probabilities greater than 0.95 for 97% of all observations. On average, fishers were in the resting state 63% of the time. Rest events averaged 5 h (4.3 SD) and occurred most often at night. The SSMs allowed us to estimate the 95% credible ellipses with a median area of 0.12 ha for 3772 unique rest events. We identified 1176 geographically distinct rest locations; 13% of locations were used on > 1 occasion and 5% were used by > 1 fisher. Females and males traveled an average of 6.7 (3.5 SD) and 7.7 (6.8 SD) km/day, respectively.
Conclusions
We demonstrated that if auxiliary data are available (e.g., accelerometer data), a twostage approach can successfully resolve both problems of latent behavioral states and GPS measurement error. Our relatively simple twostage method is repeatable, computationally efficient, and yields directly interpretable estimates of resting site locations that can be used to guide conservation decisions.
Background
Innovations in animal biologging technology in recent years have opened significant new avenues of research into animal behavior and their landscape use. Lightweight biotelemetry devices containing Global Positioning System (GPS) or similar technologies have increasingly been adopted by researchers targeting a variety of terrestrial and aquatic species to collect time series of individual animal locations. These telemetry devices often contain multiple sensors (e.g., thermometer, accelerometer), allowing time series of location data to be paired to matching time series of physiological measurements or other auxiliary information [1]. Innovations in technology have been followed by advancements in statistical methods to analyze telemetry data, resulting in a broad array of multivariate random walktype models that estimate statedependent movement parameters from location time series data [2, 3]. However, there is no generally applicable movement model for animal location data because each analysis must be tailored both to the set of research questions under study and the particular structure of a given dataset (e.g. regular or irregular recording intervals or availability of auxiliary data) [4]. Given these issues, examples of using statistical movement models to directly inform management of animal populations are limited.
Analyses of time series data arising from animal telemetry studies such as this one are often complicated by two competing challenges arising from different kinds of unobservable information. The first is that any sufficiently long time series of animal locations (“track”) is likely to result from multiple behaviors [5]. These behaviors are seldom observed directly, but each is associated with different characteristic movement patterns. Hidden Markov models (HMMs) address this challenge by assuming that observed variables (e.g., step length, turning angle) arise from a mixture of distinct probability distributions associated with latent discrete behavioral states. Serial dependence is explicitly modeled by allowing animals to transition between states at each timestep where the statedependent transition probability is a parameter to be estimated [6]. For example, Morales et al. [7] developed an HMM for elk (Cervus canadensis) where tracks were assigned to “exploratory” or “encamped” behavioral modes based on step lengths and turning angles between successive GPS positions. Because of HMMs’ flexibility and power, there are many recent examples of their application to positional data [8,9,10], and specific software tools have been developed for analysis of movement data with HMMs [11, 12]. HMMs can include or be fit entirely to nonpositional data [13,14,15], and inclusion of auxiliary data with positional data can improve inference from HMMs [1]. However, when applied to positional data, HMMs are only appropriate for time series data recorded with minimal error because measurement error can severely affect model inference when the scale of measurement error is large compared to the scale of movement [16].
Observation error and missing observations comprise the second major challenge to analysis of location data because, in both cases, the true location of an animal is not directly observed. State space models (SSMs) are flexible time series frameworks that, in the context of animal movement, treat location as a latent “state” to accommodate missing observations and estimate observation error about the true location at each timestep [17]. The term SSM is used broadly in the literature and in some cases HMMs may be considered a special case of SSMs [18]. Here, we use the term in a manner similar to Patterson et al. [4] to refer to a model where the latent state is continuous. To avoid confusion of the term “state” as used in the SSM, we refer to discrete latent HMM states as “behavioral states” and continuous latent SSM states as “latent locations”. When process and observation errors are assumed to be Gaussian, the Kalman filter (KF) can be used to efficiently fit models to animal location data [17]. For example, Johnson et al. [19] fit a linear SSM using the KF to estimate the true position from irregularly observed harbor seal (Phoca vitulina) locations. However, while Johnson et al. [19] demonstrate that SSM can be used to estimate behavioral state specific movement parameters, the behavioral state must be known in advance.
Although robust methods exist to deal with these two challenges separately, addressing both simultaneously remains difficult. For example, the HMM of Morales et al. [7] assumed that measurement error was negligible in scale relative to the observed movement, and the SSM of Johnson et al. [19] assumed that behavioral states were sufficiently described by an observed covariate. Jonsen et al. [20] developed a model accounting for both behavior switching and measurement error by allowing multiple observations to be assigned to a true location, thus allowing for the replication necessary to provide an estimate of measurement error. McClintock et al. [1] used a similar approach but included independent information on measurement error. Additionally, while closedform likelihoods can be calculated for HMMs and SSMs separately, no such expression exists for models that propose to combine features of both [21]. This necessitates use of computationally intensive approaches such as Gibbs sampling, which are not guaranteed to converge even after weeks of sampling. For these reasons, McClintock [22] proposed a twostage approach to impute missing and imperfectly observed locations using a SSM followed by a HMM to identify latent behavioral states.
We were interested in using biotelemetry data to characterize resting and movement behavior and to identify resting sites of GPScollared fishers (Pekania pennanti) in the Oregon and California Revested Lands (O&C Lands) managed by the Bureau of Land Management (BLM) in southern Oregon, USA. This is a unique region for western fisher research owing to checkerboard of public and private ownerships that causes a pronounced juxtaposition in land management practices. Fishers are a mediumsized member of the weasel family, typically associated with closedcanopy, older forests and are considered by land managers as a sensitive species when planning management actions. These territorial animals are thought to have at least one resting bout per day, often in structures such as snags or tree hollows [23]. The locations and associated vegetation of such resting locations have been described as representing “the best source of information that resources managers can use to maintain or improve habitat conditions for fishers [24].” While biotelemetry has been applied to this small forestdwelling carnivore in recent years [25, 26], including examples of behavioral state switching models [27, 28], much remains unknown in terms of their basic habitat and space use such as duration, timing and location of resting, and information such as average distance moved per day. Understanding fisher movement and resting patterns is a first step to determining what types of forest management are conducive to the persistence of fisher populations. We deployed GPS collars on 9 individuals on 21 occasions with a goal of describing basic movement ecology to inform management of this carnivore. We were particularly interested in identifying the location of resting sites as this is thought to be a limiting factor for this species [23]. For each deployment, our data consisted of regular time series of GPS fix attempts and accelerometerderived activity summaries. This species is often associated with areas of dense forest cover and resting in cavities formed in trees, snags, logs, or under the snow, which contributed to failure of a sizeable proportion of attempted GPSfixes in each time series, resulting in missing location data for some timesteps. The remaining GPS fixes were known to be observed with error.
Our research question required us to account for both latent behavioral states (i.e., separating periods of resting from periods of movement) and telemetry error to obtain robust estimates of potential resting locations. To resolve both problems, we implemented a twostage approach similar to that of McClintock [22], but in the opposite order. We chose to use the HMM first because our ultimate goal was an estimate of resting site location, which required us to first identify when a fisher was resting in order to model where it was resting. We first identified periods of resting and movement from time series of activity and partially observed step lengths using an HMM. Treating the resulting behavioral state estimates as known, we then fit an SSM to generate estimates of the true location of the animal at each timestep. Our results allowed estimation of valuable natural history and behavior information from an elusive, forest dependent, small carnivore. In particular, our goals were to use these methods to summarize three characteristics valuable to scientists and forest managers: (1) locations of putative resting locations that may be considered as focal spots for conservation; (2) characteristics of visitation and reuse by one or more resting fishers as a potential characteristic of presumed importance; and (3) resting durations and activity patterns (e.g., diurnal vs. nocturnal) to help guide field efficiencies or biases.
Methods
Study animal
Fishers are forestdwelling mesocarnivores whose range in North America and the western states has substantially declined [29, 30]. Fishers are associated with older forest elements such as large live trees, snags, and logs (e.g., [31]) and are a cavity denning obligate, thus they seek structures with specific shelterproviding cavities for parturition and kit development [32]. Fishers rest in similar structures throughout the year [23, 33], and resting locations have been used to monitor habitat availability and trends [34]. Fishers were recently proposed for listing as Threatened under the Endangered Species Act (ESA) throughout the western United States [35], although this listing was determined to be not warranted for the northern California and southern Oregon distinct population segment [36]. Fishers are designated as a Sensitive Oregon Conservation Strategy Species [37]. At the time of study initiation, little was known about the fisher population and habits in the southern Oregon Cascades (but see [38]). These fishers differ from others in the region as they were reintroduced from Minnesota and British Columbia (1961–1981 [39]), and thus are larger than other fishers in Oregon and California (male average 4.5 kg (0.8 SD); female 2.6 kg (0.2 SD)) and our scope of inference is limited to this reintroduced population.
Study area
Land ownership in our study area is a patchy mosaic of federal (BLM and United States Forest Service), state, and private (primarily commerciallyowned timberlands) lands. Thus, our study provided a unique opportunity to monitor fisher behavior across a range of variable forest management practices and across a range of seasonal variation. Pacific fishers are often associated with mixed conifer oak woodlands and are active during all times of year although snow may alter their movement patterns [40]. We conducted our study between 2015 and 2018 and during all four seasons. Elevations in this mountainous region range from 1200 to 2500 m. Forest vegetation types are primarily mixedconifer stands with predominant tree species that include white fir (Abies concolor), Douglasfir (Pseudotsuga menziesii), red fir (Abies magnifica), ponderosa pine (Pinus ponderosa), sugar pine (Pinus lambertiana), lodgepole pine (Pinus contorta), and incense cedar (Calocedrus decurrens). Natural openings include perennial meadows and frozen lakes during winter. Average annual snowfall for the Howard Prairie Dam weather station located at 1396 m elevation within the study area is 349 cm and average annual precipitation is 82 cm, which is an unusually high amount of snow for fisher populations in the western United States. Winter mean annual snow depth for the study period was 37 cm at the Howard Prairie weather station. Snow persisted from November to May each year in much of the core study area. Mean annual precipitation was 77 cm for 2015 and 2016 (Natural Resources Conservation Service, SNOTEL data 2015–2016).
Data description
We first identified areas inhabited by fishers using baited remote cameras at randomly selected locations within a systematic grid and with strategic placement (see [41] for methods). At locations with fisher detections, we trapped, anesthetized, and handled fishers using established protocols [42]. We fitted GPS collars (W500 Wildlink GPS Logger, 60 g, ATS Isanti, MN) to 9 individual fishers. Data used in this study were retrieved by recapturing the animal and downloading the archived data. We also retrieved some data remotely with UHF transmission to ensure some information was gained if an animal was not recaptured, but these remotely sampled data were not used in this study as all collars were successfully recovered. Each individual fisher was tracked for between one and four collar deployments resulting in 21 total deployments. Given an expected collar battery life ranging from approximately 30 to 120 days, we designed our study to obtain telemetry data representative of different seasons to ensure our scope of inference was not limited to a single season in case of seasonal differences. The number of days each collar was operational ranged from 23 to 170 (Table 1). We retain collar deployment within each model to assess the degree of variation between collars and individual behavior.
For each deployment, we programmed the collar to attempt a GPS location fix every 15 or 30 min, with the longer time interval generally corresponding to longerduration deployments. If the collar was unable to obtain a GPS position fix after 180 s the collar was programmed to abort the fix attempt and no location was recorded. The fix success rate varied among collar deployments with between 10 and 55% of attempted position acquisitions failing. However, failed acquisitions more commonly occurred at the end of a deployment when battery power was diminishing. Thus, we truncated each time series by eliminating all observations after the last three consecutive successful location fixes. Fix rate in remaining times series ranged from 48 to 90% (Table 1). Missing location data were interspersed throughout each remaining time series. The average duration of data with consecutive successful GPS fixes before a failure was 6.75 h, but this varied among deployments (range of means: 2.75 to 21.5 h). Fishers are territorial carnivores and observed fisher GPS locations were generally confined to areas consistent with distinct territories with limited overlap between individuals. Territories of individual fishers were consistent across deployments (Fig. 1).
Each GPS collar contained a triaxial lowg accelerometer. The Wildlink collar is programmed to summarize accelerometer measurements by an “activity” variable. To create this variable, the accelerometer is sampled every second to detect a minimum acceleration of 0.77 ms^{− 2} in any direction and activity is summarized as the percent of seconds between each attempted GPS acquisition in which this minimum acceleration was detected. For example, for collars set to a 15min interval, an activity measurement of 0.5 corresponded to 450 of 900 s for which the accelerometer recorded movement (personal communication, ATS Isanti, MN). Activity values ranged from 0.005 to 0.995 and were recorded for each timestep even if the GPS fix failed.
For each collar deployment, we had two related time series: the set of GPS fix attempts including missing GPS positions; and the set of activity measurements representing sequential pairings of consecutive GPS fix attempts. We use t as the time index for both time series with the understanding that data in the timeseries representing summaries over the interval t to t + 1 map to t. For example, the activity measurement (a_{t}) for time t = 1 represents the summary of accelerometer readings between t = 1 and t = 2. We define the vector x_{t} = [x_{1, t}, x_{2, t}]^{′} as, respectively, the observed easting and northing GPS coordinates at time t. From this time series, we calculated apparent step lengths (l_{t}) as stepwise Euclidean distance (e.g. \( {l}_t=\sqrt[2]{{\left({x}_{1,t+1}{x}_{1,t}\right)}^2+{\left({x}_{2,t+1}{x}_{2,t}\right)}^2} \)). We use the term “apparent step length” to indicate that the observed step lengths do not necessarily indicate actual movement but result in part from GPS measurement error. If at least one GPS location was missing, we recorded an NA for l_{t}. This resulted in a paired set of bivariate time series for each collar deployment (Fig. 2); the first consisting of the observed easting and northings including missing locations, and the second consisting of a complete time series of activity measurement paired with the time series of observed and missing apparent step lengths.
Hidden Markov models
In the first stage of the analysis, we fit HMMs to the time series of activity measurements and apparent step lengths of each collar deployment to identify periods of resting and periods of movement. Numerous resources cover application of HMMs to animal movement (e.g. [6]). Here, we briefly review the model structure assuming a timehomogenous process and adapting the notation of LeosBarajas and Michelot [12]. For clarity, in this section and the next, we define the models generally for a single collar deployment. We fit the models independently to each collar deployment, yielding 21 sets of estimates for all parameters defined below.
In an HMM, the observed data vector at time t, y_{t}, depends on the underlying discrete behavioral state of the animal, s_{t}. Each element of y_{t} is conditionally independent and drawn from probability distributions that depend on the behavioral state. That is, P(y_{t}s_{t} = i)~f_{i} (y_{t}), where f_{i} is the probability distribution function of y_{t} when the animal is in behavioral state i. Note that f_{i} is a multivariate joint probability distribution if y_{t} is multivariate. If one assumes the contemporaneous conditional independence of the elements of y_{t}, then f_{i} can be expressed as the product of univariate distributions. HMMs assume that s_{t} can take on a finite number of distinct states (s_{t} = 1, …, N) and that transitions between behavioral states over time occur via a firstorder Markov chain. This means that the behavioral state at time t depends only on the behavioral state at time t – 1; that is, P(s_{t}s_{t − 1}, …, s_{1}) = P(s_{t}s_{t − 1}). The probability of transitioning from behavioral state i at time t – 1 to behavioral state j at time t is γ_{i, j} = P(s_{t} = js_{t − 1} = i), which is an element of the N x N transition probability matrix, Γ. To complete specification of the model requires the initial probability distribution, δ, which for a timehomogenous process can be the stationary distribution defined such that δ = δΓ. The marginal likelihood of the observed data can then be written as:
where P(y_{t}) = diag(f_{1}(y_{t}), …f_{N}(y_{t})) and 1 is a lengthN vector of 1’s.
In terrestrial applications, animal movement HMMs have typically been applied to animal location data by decomposing GPS coordinates and modelling the step lengths and turning angles between consecutive locations (e.g. [7]). We reasoned that fishers would be almost entirely stationary during times of resting but recognized that our location data were subject to GPS measurement error, resulting in nonzero apparent step lengths and spurious turning angles [43]. Because GPS measurement error affects turning angles to a greater degree than step lengths [44], we chose to use only the accelerometerderived activity measurements and apparent step lengths. We expected that resting would be associated with lower activity measurements and shorter apparent step lengths. Because a_{t} was constrained between 0 and 1, we modeled activity as Beta distributed, where 0 < μ_{i} < 1 is the mean and ϕ_{i} > 0 is the precision of a Beta distribution given behavioral state i under the parameterization of Ferrari and CribariNeto [45]. We used a lognormal distribution for apparent step lengths, where m_{i} is the log mean and σ_{i} > 0 is the log standard deviation of a lognormal distribution given behavioral state i. Assuming independence between activity and step length within each behavioral state, the conditional joint likelihood when both variables are observed is: f_{i}(y_{t} = {a_{t}, l_{t}}) = f_{i}(a_{t})f_{i}(l_{t}) and when only activity measurements are observed is: f_{i}(y_{t} = {a_{t}, NA}) = f_{i}(a_{t}). When x_{t} = x_{t + 1,}l_{t} = 0, which is inadmissible under the log normal distribution. In these instances, we set l_{t} = 1 for the HMM model, which affected 193 of 57,788 observed l_{t} measurements. Because we were mainly interested in separating resting from movement, we fit a twostate model under the assumption that all other types of behaviors that involved movement would be captured by the remaining behavioral state. We used a single timehomogenous parameter, ψ_{i} as the probability of remaining in state i and let γ_{i, i} = ψ_{i} and γ_{i, j} = 1 − ψ_{i}. We refer to ψ_{i} as the “state persistence probability.”
We identified resting periods by estimating the underlying behavioral state sequence given the complete observed data and the estimated parameters, Θ. We used the forwardfiltering, backwardsampling (FFBS) algorithm to globally decode the posterior marginal probability of the behavioral state sequence [46]. While the Viterbi algorithm has commonly been used to decode the most likely state sequence from HMMs [6], we chose to use the FFBS over the Viterbi because the Viterbi calculates the jointly most probable state sequence, whereas the FFBS calculates the posterior state assignment probability for each observation. We wanted to be conservative in assigning observations to the resting state, so we chose the FFBS because it gave a wider estimate of uncertainty of the true state for each timestep. We assumed that the behavioral state associated with lower activity corresponded to resting, which we labeled: s_{t} = 1. For each collar deployment, we used the posterior distribution of decoded behavioral state sequences to assign the inferred behavioral state, \( {s}_t^{\ast }=0\ \mathrm{if}\ P\left({s}_t=1\ \ {\boldsymbol{y}}_{1,}\dots, {\boldsymbol{y}}_T,\boldsymbol{\varTheta} \right)>0.95;\mathrm{and}\ 1\ \mathrm{otherwise}. \)
State space models
Assuming the inferred behavioral states as given, we used a linear Gaussian SSM to estimate fisher resting locations and movement between resting locations while accounting for GPS measurement error. The linear Gaussian SSM assumes two conditional relationships: first, that the time series of observed locations, x_{t} depends on a sequence of unobserved continuous latent locations α_{t}; and second, that α_{t + 1} depends on α_{t}. Like HMMs, SSMs have been widely applied to animal movement data (e.g. [19]) to take advantage of the KF [17] to efficiently estimate parameters related to movement and measurement error and to estimate (“smooth”) the time series of latent locations from imperfectly observed coordinate data. Here, we briefly review the mathematical details as they pertain to our implementation.
The linear Gaussian SSM consists of two related equations, the observation equation and the system equation:
where x_{t} are the observation locations, α_{t} are the latent locations, Z_{t} is the design matrix that maps α_{t} to x_{t}, ε_{t} is the meanzero normallydistributed observation error vector with covariance matrix H_{t}, U_{t} is the transition matrix mapping α_{t} to α_{t + 1}, η_{t} is the meanzero normallydistributed process error vector with covariance matrix Q_{t}, and R_{t} is the system selection matrix.
We parameterized our SSMs to alternate between periods of movement and resting based on the inferred behavioral state (\( {s}_t^{\ast } \)) arising from the HMM. We assumed that the fisher was stationary during periods of resting (\( {s}_t^{\ast }=0 \)), meaning that α_{t + 1} = α_{t}. During periods of movement (\( {s}_t^{\ast }=1 \)), we assumed that fisher moved according to a firstdifference correlated random walk (DCRW) [21]. To account for the DCRW we set α_{t} = [α_{1, t}, α_{2, t}, α_{1, t − 1}, α_{2, t − 1}]^{′} where α_{1, t} and α_{2, t} are, respectively, the latent easting and northing coordinates of the fisher at time t. For all t, we set the design matrix \( {\boldsymbol{Z}}_{\boldsymbol{t}}=\left[\begin{array}{cccc}1& 0& 0& 0\\ {}0& 1& 0& 0\end{array}\right] \) to reduce the fourelement α_{t} vector to the twoelement x_{t}. We assumed a constant GPS error over time and independent in each direction, that is H_{t} = \( \left[\begin{array}{cc}h& 0\\ {}0& h\end{array}\right] \) where h is the observation error standard deviation. Together, the transition matrix and system selection matrices encoded either nonmovement or a DCRW with correlation coefficient, ρ, as:
Finally, we assumed a twoelement system disturbance vector, η_{t} = [η_{1, t}, η_{2, t}], representing movement in the easting and northing directions, respectively, with covariance matrix Q_{t} = \( \left[\begin{array}{cc}q& 0\\ {}0& q\end{array}\right] \) where q is the process error (movement distance) standard deviation. Thus, when the fisher is inferred to be moving (\( {s}_t^{\ast }=1 \)) and using α_{t} = [α_{1, t}, α_{2, t}, α_{1, t − 1}, α_{2, t − 1}]^{′} as the frame of reference, these relationships yield:
which states that fisher location at time t + 1 is equal to the location at time t plus some contribution of its velocity from the previous timestep, as represented by the firstdifference, α_{1, t} − α_{1, t − 1}, plus some additional movement deviating from the velocity and direction of the last time step. When we inferred the fisher to be resting (\( {s}_t^{\ast }=0 \)), we have simply α_{t + 1} = [α_{1, t}, α_{2, t}, α_{1, t}, α_{2, t}]^{′}, which states that the fisher did not move in the previous timestep.
Model fitting and inference
We fit both models in the Stan probabilistic programming language [47]. While several standalone R packages have been developed for animal movement models (e.g., [12]) and Bayesian models have been published using BUGStype language (e.g., [20]), we chose to use Stan because of its speed and customizability [48]. We completed all data formatting and model fitting in R version 4.0.1 [49] using the rstan package [50]. We adapted the Stan code of LeosBarajas and Michelot [13] for the HMMs and that of Arnold [51] for the SSMs.
We first fit independent HMM models to the time series of activity measurements and apparent step lengths for each collar deployment. We used identical semiinformative priors for each time series: logit(μ_{1})~N(−3, 0.5); logit(μ_{2})~N(0, 1); log(ϕ_{1})~N(1, 1); log(ϕ_{2})~N(0, 1); m_{1}~N(3, 1); m_{2}~N(5, 1); σ_{1}~ ∣ N(3, 1) ; σ_{2}~ ∣ N(3, 1) ; ψ_{1}~Unif(0, 1); ψ_{2}~Unif(0, 1). For each model, we ran four Hamiltonian Monte Carlo (HMC) chains with a warmup of 1000 iterations and sampled for 1000 iterations. We used semiinformative priors with the intention to avoid identifiability issues common to mixture models [52]. We selected priors based on preliminary graphical analysis of apparent step lengths and activity measurements. While our data structure of multiple animals with one or more collar deployments lends itself to a hierarchical model wherein the movement parameters for each time series are drawn from a common distribution, we chose to fit the model independently because our main focus for inference was on the location of each rest event rather than on variation of higherlevel parameters across animals. Given the size of our data, we did not expect a hierarchical model to improve state assignment estimates [53].
Because general behavioral information on fishers is limited, we were interested in describing diurnal patterns of resting. We used the R package suncalc to compute the positions of the sun for our study area on days when GPS collars were deployed [54]. We defined eight categories of daylight: (1) “night after solar nadir” from the solar nadir to the beginning of astronomical twilight; (2) “morning twilight” from the beginning of astronomical twilight to dawn (the start of civil twilight); (3) “morning” from dawn to the end of the golden hour (approximately 1 h after the end of sunrise); (4) “day before solar noon” from the end of morning to solar noon; (5) “day after solar noon” from solar noon to the start of the evening golden hour; (6) “evening” from the end of the day through dusk (evening civil twilight); (7) “evening twilight” from dusk to the end of astronomical twilight; and (8) “night before solar nadir” from the end of twilight to the solar nadir. For each category, we calculated the proportion of each HMM time series classified as resting out of all time series observations.
We next fit independent SSM models to time series of GPS fix attempts and the inferred state assignment from the HMM. For the SSM models, we used identical weaklyinformative priors that were uniform over the support of the parameters. We again ran four HMC chains with a warmup of 1000 iterations and sampled for 1000 iterations for each model. After fitting the SSM, we estimated the true location at each timestep. We used the meancorrection simulation smoother to draw samples from the posterior of true locations [17]. For each rest event, defined as a consecutive set of timestamps for which \( {s}_t^{\ast }=0 \), we used the ks package in R to estimate the 50, 80, and 95% credible ellipses [55] for the location associated with each rest event. We were also interested in identifying potential reuse of rest sites by one or multiple fishers because reuse of landscape features as resting sites or use by multiple fishers has not been welldocumented and such features may be of greater conservation value than resting sites that are used only once. Out of all rest events identified across the 21 deployments, we selected those with at least one successful GPS fix, because preliminary analysis of the results showed that rest events with no successful GPS fixes had very large credible ellipses. We identified spatially unique rest locations by spatially joining overlapping 95% credible ellipses. We defined reuse of a rest location if rest events for two different fishers or rest events of the same fisher were separated by at least 24 h and had overlapping 95% credible ellipses.
We conducted a partial field validation of our method using the integrated VHF transmitter on the GPS collars which we used to opportunistically locate rest and den locations during field data collection. We acoustically determined whether collared fishers were inactive as indicated by a consistent VHF signal heard through a passive telemetry receiver (R1000, Communication Specialists, Orange, California, USA). If the fisher was inactive for > 2 min, we attempted to locate the resting or denning fisher by “walking in” and decreasing the gain, or power, on the receiver as we approached the structure. The field observer rated their confidence of the location between 1 and 5 (1 = fisher visually seen in structure, 2 = signal consistent with one structure, 3 = fisher within 2–3 adjacent structures, 4 = visual of animal leaving but structure not identified, 5 = within 100 m of animal but structure not identified). We compared the VHF identified rest sites to our model output as a partial validation of our approach. For this comparison, we restricted our sample to VHFlocations with a confidence of 3 or less.
Results
Hidden Markov models
The HMM models fit independently to 21 collar deployments converged successfully to similar estimates that clearly discriminated between two latent behavioral states. The number of effective posterior samples for the main parameters of the models ranged from 777 to 7126 (mean 4142) and each parameter had \( \hat{R}<1.01 \), indicating successful convergence. Although the models were fit independently to each collar deployment time series, parameter estimates were qualitatively similar across all 21 models (Fig. 3), which suggests that the models describe a consistent suite of behaviors across the population of 9 fishers. Similarly for each collar deployment, the HMM models showed strong separation in the posterior probability of being in the presumed resting (s_{t} = 1) or moving states (s_{t} = 2), as only 3% of all observations had 0.05< P(s_{t} = 1) < 0.95 (range 1 to 7%; Table 2).
The resting state was strongly associated with small values of a_{t}. Across 21 collar deployments, the posterior mean of the Beta distribution mean for the resting state (μ_{1}) ranged from 0.01 to 0.06 (average 0.03) and precision (ϕ_{1}) ranged from 9.5 to 164 (average 50). These parameters resulted in a set of Beta distributions for the resting state with a mode near 0 and with little probability density for activity values greater than 0.25 (Fig. 4). For the moving state, the posterior mean of μ_{2} ranged from 0.52 to 0.75 (average 0.65). However, precisions (ϕ_{2}) were less for the moving state, ranging from 1.9 to 11.3 (average 5) and resulting in more diffuse distributions (Fig. 4a).
Smaller apparent step lengths, l_{t}, were also associated with the resting state for all models, and greater values of l_{t} were associated with the moving state. Across collar deployments, the posterior mean of median apparent step length \( \left({e}^{m_i}\right) \) ranged from 7.5 to 20 m (average 13 m) for the resting state and from 59 to 383 m (149) for the moving state. Log standard deviations (σ) diverged less between states, with posterior means ranging from 0.8 to 1.0 (average 0.9) for the resting state and 0.8 to 1.7 (average 1.4) for the moving state. This resulted in more diffuse apparent step length distributions for the moving state than the resting state (Fig. 4b). However, apparent step lengths for observations assigned to moving state exhibited negative skewness which was not fully captured by the lognormal distribution (Supplemental Figure 1).
Both a_{t} and l_{t} displayed bimodality when aggregated across time. Because this bimodality was apparent for timeaggregated individual time series, we display data aggregated across all 21 collar deployments in Fig. 4c and d. This bimodality and the correlation between a_{t} and l_{t} (Fig. 4e) contributed to the model’s ability to separate between states. For example, the most common value of a_{t} across all 21 collar deployments was the minimum 0.005, which was recorded for 25,411 of 81,024 observations. Of these, only 65 had a P(s_{t} = 1) < 0.95. While we display the timeaggregated data (Fig. 4) to demonstrate the correspondence between observed bimodality and estimated statedependent probability distribution, it is important to note that this view ignores the time dependence that affects any given observation’s state probability.
The duration of rest events varied within a narrow range across collar deployments and individual fishers. Because we used a timehomogenous model, the duration of resting events is a direct consequence of the probability of remaining in the resting state, ψ_{1}. The posterior mean of this parameter across collar deployments ranged from 0.89 to 0.96 (average 0.93) for collars set to 15 min increments which corresponds to expected rest event durations of 2.3–6.3 h. For collars set to 30 min increments, the posterior mean of ψ_{1} ranged from 0.88 to 0.92 (average 0.91) which corresponds to expected rest event durations of 4.1–6.2 h (Fig. 3). The probability of remaining in the moving state, ψ_{2}, was more variable across deployments but generally lower than ψ_{1} (average posterior mean of 0.88 and 0.79 for 15 and 30 min increment deployments, respectively). As a result, we found that fishers spent more time in the resting state than in the moving state (Table 2) except for summer collar deployments for fishers F01T and F03T. We identified 4068 unique rest events across the 21 collar deployments. The duration of individual rest events was variable, but most lasted between 2 and 8 h (Table 2, Supplemental Figure 2).
The proportion of time resting during each part of the day was variable across fishers and seasons, but in general fishers were least active during the night and most active during the morning (Fig. 5). The proportion of time resting was highest during the night after solar nadir, with a median proportion across collar deployments of 0.78 (range 0.52 to 0.94). For the morning, the median proportion of time resting across collar deployments was 0.36 (range 0.13 to 0.67). The median proportion of time resting during the day after solar noon and during the evening were second and third highest across collar deployments, each at 0.68, but had the highest variability across deployments as reflected in the wide range of proportions (0.11 to 0.95 for day after solar noon, and 0.13 to 0.91 for evening). This may reflect seasonal changes in fisher behavior, as summer deployments of 2016 were on the low end of this range and fall and winter deployments were on the high end. Both night categories had the narrowest range of values in the proportion of time resting (0.42 to 0.84 for night before the solar nadir), which suggests a consistent tendency for our observed fishers to rest at night.
State space models
As with HMM models, the SSM models fit independently to 21 collars deployment converged successfully to similar estimates. The number of effective posterior samples for the main parameters of the models ranged from 2917 to 4961 (mean 3998) and each parameter had \( \hat{R}<1.01 \). Parameter estimates displayed broadly similar pattern across all 21 models with observation error a magnitude of order less than process error and positive autocorrelation (Fig. 6).
The precision of any given rest event location estimate depended on observation error standard deviation and the number of observed GPS fixes. Estimates of observation error standard deviation were on the same order as field tests of GPS precision [56, 57], with an average of posterior means of 20 m (range 7 to 36 m), similar to the average of 21 m estimated precision during field trials [56]. The SSM was able to estimate a location for rest events in the absence of any observed GPS fix with the area 95% credible ellipses for such events (n = 850) ranging from 11.8 ha to 2709 ha (Fig. 7). Rest events with the greatest precision were those with at least one successful fix (n = 3218) with a median area of 95% credible ellipses of 0.06 ha (range 19 m^{2} to 2.5 ha). While credible ellipses for location estimates of rest events without observed GPS fixes were too large to identify meaningful overlap with other rest events, in all instances, the SSMs constrained the possible location of each fisher within bounds of our study area (Fig. 8, top panel).
We found evidence of reuse of rest sites by individuals fishers over time and by multiple fishers (Fig. 8). Of rest events with at least one successful GPS fix, we identified 1176 geographically distinct areas by spatially joining overlapping 95% credible ellipse. Of these, we found 148 locations were revisited by a fisher after at least 24 h since last use with a median elapsed time between revisits within a single collar deployment of 7.6 days. We found 57 locations where the 95% credible ellipse of the rest event for one fisher overlapped with the 95% credible ellipse for another fisher. Of these, we found two instances of the 95% credible ellipses rest events for three fishers overlapping spatially.
During periods of movement, the scale of fisher movement was an order of magnitude greater than the GPS observation error. The posterior mean standard deviation of process error, which describes variance in the step length between successive fisher locations during periods of movement, averaged 165 m (range 115 m to 233 m) for the 15 collar deployments set at 15 min intervals and 428 m (range 256 m to 618 m) for collar deployments set at 30 min intervals. When fishers were moving, there was a tendency for them to continue moving in the same direction as the previous timestep, as evidenced by the positive correlation in movement direction (ρ average posterior mean 0.437).
The simulation smoother of the SSM also estimated the distance moved at each timestep when the fisher was in the moving state. The average posterior mean cumulative distance moved per day (Supplemental Figure 3) ranged from 3.6 km/day (F03T Fall 2015) to 12.8 km/day (M09T Winter 2017) with an overall average of 6.7 km/day (3.5 SD) for females and 7.7 km/day (6.8 SD) for males. The total distance moved per day was interspersed with resting events and so daily movement distances where made up of one or more periods of movement between resting. For these movement events, the average posterior mean cumulative distance travelled between rest events ranged from 0.6 km (F03T Fall15) to 4.2 km (M09T Winter17). For the two female fishers with multiple deployments across seasons (F01T and F03T), we found a similar pattern of greater daily cumulative distance travelled during summer deployments than during fall and winter deployments. This corresponds with a lesser proportion of time resting during summer deployments for these fishers (Table 2), particularly during the day (Fig. 5). However, F02T showed similar distances travelled between fall and summer deployments (Supplemental Figure 3), which suggests that the seasonal patterns observed in F01T and F03T may be due to variation among individuals rather than a general pattern of behavior by female fishers. The greatest daily distance travelled was for fisher M03T, who traversed an estimated 47 km (90% credible interval: 46–48 km) on March 25, 2017. Similar long daily distances were estimated for fisher M09T (35 km on March 11, 2018) and M08T (32 km on March 22, 2017). In contrast, the longest daily distance observed for a female fisher was 17.6 km (90% credible interval: 17.1–18.1 km) for F01T. We suspect that long distance movement of these males reflects mating behavior. During mating season, male fishers may travel to many females during their short period of estrus [58]. The timing of the longdistance movement we observed is consistent with seasonal mating as the average natal den initiation date we observed for this fisher population in a separate study was March 30 [56] and mating takes place approximately 10 days after birth [58].
Field validation
We identified 79 suspected resting events using VHF telemetry across 21 deployments. Of these 12 (14%) did not correspond temporally with a resting event identified by the HMM models, but for 5 of these 12 the HMM estimated a P(s_{t} = 1) ≥ 0.6. For those 5 events, the location identified via VHF telemetry was within the 95% credible ellipse of the HMM/SSM identified rest event nearest in time to VHF record. Of those which could be matched temporally to an HMM/SSM rest event, 8 (9%) were matched to rest events with no successful GPS fixes and thus large credible ellipses. For all but one of these, the VHF rest location was within the 95% credible ellipse (median area: 25 ha). 22 (25%) of VHF identified rest events were located within the 95% credible ellipse of the temporally corresponding HMM/SSM rest event (median area: 265 m^{2}). For another 27 (32%) the VHF identified location was within 25 m (median distance from edge: 6 m) of the nearest edge of the corresponding 95% credible ellipse (median area: 141 m^{2}). Finally, 10 (12%) of VHF identified locations were greater than 25 m from the nearest edge (median distance from edge: 52 m) of the corresponding 95% credible ellipse with a maximum distance of 433 m (median). Examples from each of these categories are depicted in the Supplemental Figure 4.
Discussion
Our twostage approach to modelling fisher GPS location and accelerometer data allowed us to overcome the dual problems of missing and imperfectly observed location data and unobserved behavioral states to obtain credible estimates of fisher resting sites. The data we analyzed represent the largest published GPS dataset for fishers in North America and allowed us to make inferences across individuals and within individuals across seasons. Despite fitting models independently to each dataset, we found remarkable agreement among the parameter estimates for both the HMMs and SSMs. This suggests that our models captured general patterns of fisher behavior and provided specific estimates of the location of resting sites for these individuals. Thus, we were able to both improve the general state of science on fisher ecology and produce directly actionable information to improve management of the study population.
The models described here allows use of the complete set of collected data to support analysis of fisher use of their landscape. A central tenet of ecology is that animals make movement choices to maximize their fitness [59]. We found that fisher spend a substantial amount of time resting and thus we infer that rest sites are important to their life history and constitute spatially explicit unique habitat components. The use of multiple, welldistributed rest sites suggests that regular access to a rest site is also important. We surmise that rest events have important influence on fisher fitness and therefore the physical location (e.g. slope, aspect), and the surrounding vegetation structure and composition, of the rest structure contribute to the fitness. The precise estimate of rest sites provided by our models gives a complete accounting for likely rest sites from which abiotic and biotic features can be summarized and evaluated. This could lead to a better description of rest site that could facilitate management decisions intended to retain rest site function or create potential rest sites (e.g. selective slash pile retention or creation). The data here should be considered to represent a minimum number and spacing of rest sites because of the relatively short duration of monitoring and the limited number of reused rest sites. We suspect that with a greater monitoring period additional rest sites could be identified.
By using accelerometerderived activity measurements matched to the sampling interval of GPS positions, we were able to successfully decode latent behavioral states despite GPS measurement error. HMMs and the closely related hidden semiMarkov models [60] applied to location data alone are valid only when measurement error is negligible and GPS position fix rate is high. Measurement error in excess of approximately 10% of the scale of movement can severely affect parameter estimates for movement models [16]. GPS measurement error and low fix rate affects both step lengths and turning angles metrics, but turning angles are more affected [43, 44]. Because accelerometers require only battery power and data storage, they are unaffected by the error associated with satellite and telemetry positioning methods, which rely on transmission of signals between the monitoring device and satellites or telemetry receivers. Accelerometers can collect large amounts of data that can be summarized into multiple statistics appropriate for application in HMMs [13]. The activity measure we used for this analysis was ideal for our purposes because it required only limited battery power and data storage, matched directly to sampling interval of GPS positions resulting in a complete time series with no missing observations, and could be readily modelled with the flexible Beta distribution. While we understood that apparent step lengths were affected by measurement error, the bimodality observed in both metrics and correlation between the two metrics suggest that apparent step lengths provided meaningful information to the HMM likelihood, although, we note that the lognormal distribution did not adequately capture the negative skew of steplengths for observations assigned to the moving statement. The joint bimodality in the observed data probably contributed to the strong separation in the posterior state assignment probability, where greater than 95% of observations had a state assignment probability in excess of 0.95.
Although the HMMs allowed us to estimate when each fisher was resting, the SSMs allowed us to quantify uncertainty about the location of fisher rest events to achieve our main goal of identifying potential resting sites for the study population. Fishers use forest features such as snags and cavities for resting sites which have been described as forest features of high conservation value [38], but relative importance of such structures is unknown. Camera traps and direct observation through VHF telemetry can confirm use of such features by fishers, but sampling of these sites is usually opportunistic and time consuming. For instance, 9 years of VHF and camera monitoring of 127 fishers led to a similar number of rest occasions (1040) and locations (933) and slightly lower estimates of reuse rates, 9.3% for females and 10.4% for males [61]. Although deducing the number of rest locations was not the goal of either study, our methods provided a similar amount of data to address this question, with a complete resting and movement chronology during the period of observation. Our estimate that 15% of rest locations were reused either by a single fisher or by multiple fishers may imply that these locations may be more important than those used only once. Overlapping, credible ellipses for rest sites suggest potential reuse of specific landscape features that may be more important for fisher ecology than features used only once. However, the overall number of rest locations was much greater than those reused, which suggests that fishers will use many features on the landscape for resting.
Our VHF telemetry field validation showed reasonable agreement with the GPS identified rest events, although far fewer rest events were identified using VHF telemetry. While only about half (43) of the VHFidentified rest locations (85) were inside or within 10 m of the edge of the temporally corresponding rest event 95% credible ellipse, it is not clear that disagreement between the two methods can be attributed to error in one method alone. For example, during our analysis we found a field data entry error on the time of observation that after correcting resulted in a match between the methods. For the two VHF rest events with the greatest distance from the edge of the credible ellipse for temporally corresponding HMM/SSM rest event, the HMM/SSM identified a resting location corresponding to the VHF location but at a different time, but we could not confirm this resulted from a data entry error. On the other hand, we found that the areas of rest site credible ellipses corresponding to VHF rest events which were nearby but not within the credible ellipse tended to be smaller than those for which the VHF rest event was within the ellipse. This could be due to violations of the assumptions of the SSM resulting in location estimates that were too precise. Adjustments to the SSM model, for example by allowing for bias or correlation in GPS measurement error, or nonGaussian measurement error would have results in more diffuse credible ellipses which would have resulted in greater agreement between the two methods. We propose that direct observation from crewbased telemetry, cameras, and GPS data would best provide an understanding of which structures fishers perceive to be most important for survival and recruitment as the SSMs produced a catalogue of resting event location estimates that can be used to prioritize confirmatory field surveys. Further field validation of the rest sites identified by this study is planned.
Our SSM formulation allowed us to estimate the location of the fisher for each timestep of collar deployment, even when a GPS fix was not acquired. A broadrange of models have been developed to deal with the problem of apparent temporal irregularity in GPS telemetry as a result of missing GPS fixes by parameterizing animal movement as a continuous process that is observed at specific timepoints [62]. For example, continuous time methods have been developed to account for measurement error of observed locations [19], to account for stateswitching [63] and have been applied to fishers [27]. Recently, Michelot and Blackwell [64] introduced extend continuous time models to allow for switching between behavioral states between any two irregularly spaced observations. These approaches, particularly that of Michelot and Blackwell [64], are flexible and powerful due to their ability to accommodate correlation across multiple scales and accommodate irregularly observed data. However, because we had an activity measurement and thus a HMMbased estimate of the behavioral state for each scheduled GPS fix, we chose to parameterize the SSMs in discrete time and use the simulation smoother to estimate the location of the fisher for each timestep even if no GPS fix was acquired. This resulted in a complete timeseries of both behavioral state and location estimates. The error for time steps with no GPS fix was greater than those with a successful fix, but because we assumed that fishers remained stationary for rest events, rest events with at least one successful fix had markedly greater location estimates than those with none. ‘.
While the relatively simple paired HMMs and SSMs we fit to the data allowed us to achieve our research and management goals, the data structure and results suggest several possible model expansions. We collected location and activity data from nine individual fishers, four females and five males, with collar deployments representing different seasons. These data could be fit to a hierarchical (mixed) model that assumes that individual parameters–HMM state persistence probabilities, HMM beta and lognormal mean and variance terms, or SSM process and observation error standard deviations–arise from common probability distribution across individuals and deployments. Such a model would allow for partial pooling of parameter estimates and quantify variation within and among individuals [65], but at the cost of increased model complexity and computation time. While adopting a hierarchical model allows inference across individuals, it may not improve estimates of state assignment when individual track data are large [53]. Additionally, by modelling all collar deployments in the same model, it would be possible to look for seasonal differences in state persistence. While we did not observe obvious difference in state persistence across seasons, that was not a goal of this study and so we cannot rule out such differences. A second potential avenue for model expansion is inclusion of exogenous variables and temporal heterogeneity HMMs or SSMs. For example, we identified a nocturnal pattern in the proportion of time fishers spent in the resting state. Nocturnal rest structures and microsites were more specific for Pacific martens, dominated by cavities and chambers [66], which were likely more beneficial for thermoregulation than similar structures used during the day. Nonetheless, finding such structures with human observers is nearly exclusively completed during the day. Understanding influence and timing of fishers’ preferential resting periods could be more directly modelled by allowing state persistence probabilities in the HMMs to vary based on time of day [8]. Additionally, the observed lack of fit of apparent steps lengths to the lognormal distribution in the moving state may be indicative of an additional behavioral state or of temporal heterogeneity in movement that could be better captured by a threestate and/or timeinhomogeneous HMM. Because we were not interested in inferring anything about step lengths in the HMM step due to the confounding of measurement error and because our focus on this analysis was on identifying resting sites, we were less concerned about this lack of fit than if our inferential goal was on characterizing fisher movement behavior. For the SSMs, one area of expansion is to allow measurement error to vary based on recorded GPS quality measurements, such as horizontal dilution of precision. Similarly, we could explore whether fishers move greater distances during the day by allowing process error to vary based on time of day or weather. Conceivably, spatiallyvarying variables (e.g., local forest structure, distance from roads) can also affect fisher movement and behavior. However, including these variables is more complicated because of GPS measurement error and missing fixes. One way to address this problem involves discretizing the continuous spatial state into a spatial grid [18]. However, grid cell size would need to be carefully chosen to ensure that spatial covariates (e.g., forest stand type and structure) are approximately homogenous within a cell and are at scales meaningful to forest management and fisher ecology and to ensure reasonable computation time. Finally, reuse of rest sites suggests an alternative HMM parameterization involving “centers of attraction” where fishers may switch between a moving state and a “resting” state centered around a particular resting site [67]. This approach has the advantage of estimating the location of only a single center of attraction for each reused rest site, as opposed to our approach of identifying reuse that relied on the overlap of unconstrained location estimates. However, because the number of centers of attraction would need to be specified a priori in such a model, our approach remains useful as a first step to inform the construction of such a model. Further, fishers moved an average of 7.4 km daily, suggesting the capacity to move among most rest structures within a short period. If their territory boundaries are patrolled regularly, similar to every 4.6 days on average for Pacific martens (Martes caurina) [68], then territorial boundaries may also be coded as a movement attraction boundary.
Conclusions
Here, we demonstrated that a relatively simple twostage approach can overcome the twin challenges of latent behavioral states and measurement error in the context of management of sensitive species. By leveraging accelerometer data to resolve latent states in a set of HMMs, we were able to use HMM output in a set of SSMs to quantify location error to produce a set of rest event location estimates. The models provided both basic information on fisher ecology such as rest event duration, diurnal patterns of resting, cumulative distance moved between rest events, and specific information on locations of rest events that can be used to inform conservation decisions. Our relatively simple method produced rapid results and can be readily applied to additional collar deployments. Given the flexibility of both HMMs and SSMs, either stage of the model can also be expanded to include hierarchical structure or appropriate covariates for additional ecological insights.
Availability of data and materials
The dataset supporting conclusion in this article is included with the article and its additional file. Complete R and Stan code for reproducing the statistical analyses in this manuscript is provided in Additional file 2. Location data have been centered to remove reference to physical locations, but original data is available via email.
Change history
23 April 2021
A Correction to this paper has been published: https://doi.org/10.1186/s4046202100262w
Abbreviations
 HMM:

Hidden Markov model
 SSM:

State space model
 DCRW:

First difference correlated random walk
 FFBS:

Forwardfiltering, backwardsampling
 KF:

Kalman filter
 GPS:

Global Positioning System
 O&C Lands:

Oregon and California Railroad Revested Lands
 HMC:

Hamiltonian Monte Carlo
 BLM:

Bureau of Land Management
References
 1.
McClintock BT, London JM, Cameron MF, Boveng PL. Bridging the gaps in animal movement: hidden behaviors and ecological relationships revealed by integrated data streams. Ecosphere. 2017;8(3):e01751. https://doi.org/10.1002/ecs2.1751.
 2.
Hooten MB, Johnson DS, Mcclintock BT, Morales JM. Animal Movement: Statistical Models for Telemetry Data. 1st ed. Boca Raton: CRC Press; 2017. [cited 2020 Aug 11]. https://www.taylorfrancis.com/books/9781466582156
 3.
Jonsen ID, Basson M, Bestley S, Bravington MV, Patterson TA, Pedersen MW, et al. Statespace models for biologgers: a methodological road map. Deep Sea Res Part II Top Stud Oceanogr. 2013;88–89:34–46.
 4.
Patterson TA, Parton A, Langrock R, Blackwell PG, Thomas L, King R. Statistical modelling of individual animal movement: an overview of key methods and a discussion of practical challenges. ArXiv160307511 QBio Stat. 2017; [cited 2020 Aug 11]; http://arxiv.org/abs/1603.07511.
 5.
Patterson TA, Basson M, Bravington MV, Gunn JS. Classifying movement behaviour in relation to environmental conditions using hidden Markov models. J Anim Ecol. 2009;78(6):1113–23. https://doi.org/10.1111/j.13652656.2009.01583.x.
 6.
Zucchini W, MacDonald IL, Langrock R. Hidden markov models for time series: an introduction using R. 2nd ed. Boca Raton: CRC Press, Taylor & Francis Group; 2016.
 7.
Morales JM, Haydon DT, Frair J, Holsinger KE, Fryxell JM. Extracting more out of relocation data: building movement models as mixtures of random walks. Ecology. 2004;85(9):2436–45. https://doi.org/10.1890/030269.
 8.
Li M, Bolker BM. Incorporating periodic variability in hidden Markov models for animal movement. Mov Ecol. 2017;5(1):1. https://doi.org/10.1186/s4046201600936.
 9.
McKellar AE, Langrock R, Walters JR, Kesler DC. Using mixed hidden Markov models to examine behavioral states in a cooperatively breeding bird. Behav Ecol. 2015;26(1):148–57. https://doi.org/10.1093/beheco/aru171.
 10.
Whoriskey K, AugerMéthé M, Albertsen CM, Whoriskey FG, Binder TR, Krueger CC, et al. A hidden Markov movement model for rapidly identifying behavioral states from animal tracks. Ecol Evol. 2017;7(7):2112–21. https://doi.org/10.1002/ece3.2795.
 11.
LeosBarajas V, Michelot T. An Introduction to Animal Movement Modeling with Hidden Markov Models using Stan for Bayesian Inference. ArXiv180610639 QBio Stat. 2018; [cited 2020 Aug 11]; http://arxiv.org/abs/1806.10639.
 12.
Michelot T, Langrock R, Patterson TA. moveHMM: an R package for the statistical modelling of animal movement data using hidden Markov models. Methods Ecol Evol. 2016;7(11):1308–15. https://doi.org/10.1111/2041210X.12578 McInerny G, editor.
 13.
LeosBarajas V, Photopoulou T, Langrock R, Patterson TA, Watanabe YY, Murgatroyd M, et al. Analysis of animal accelerometer data using hidden Markov models. O’Hara RB, editor. Methods Ecol Evol. 2017;8(2):161–73. https://doi.org/10.1111/2041210X.12657.
 14.
Pedersen MW, Righton D, Thygesen UH, Andersen KH, Madsen H. Geolocation of North Sea cod (Gadus morhua) using hidden Markov models and behavioural switching. Can J Fish Aquat Sci. 2008;65(11):2367–77. https://doi.org/10.1139/F08144.
 15.
Phillips JS, Patterson TA, Leroy B, Pilling GM, Nicol SJ. Objective classification of latent behavioral states in biologging data using multivariatenormal hidden Markov models. Ecol Appl. 2015;25(5):1244–58. https://doi.org/10.1890/140862.1.
 16.
Bradshaw CJA, Sims DW, Hays GC. Measurement error causes scaledependent threshold erosion of biological signals in animal movement data. Ecol Appl. 2007;17(2):628–38. https://doi.org/10.1890/060964.
 17.
Durbin J, Koopman SJ. Time series analysis by state space methods. 2nd ed. Oxford: Oxford University Press; 2012. https://doi.org/10.1093/acprof:oso/9780199641178.001.0001.
 18.
Pedersen MW, Patterson TA, Thygesen UH, Madsen H. Estimating animal behavior and residency from movement data. Oikos. 2011;120(9):1281–90. https://doi.org/10.1111/j.16000706.2011.19044.x.
 19.
Johnson DS, London JM, Lea MA, Durban JW. Continuoustime correlated random walk model for animal telemetry data. Ecology. 2008;89(5):1208–15. https://doi.org/10.1890/071032.1.
 20.
Jonsen ID, Flemming JM, Myers RA. Robust state–space modeling of animal movement data. Ecology. 2005;86(11):2874–80. https://doi.org/10.1890/041852.
 21.
Ghahramani Z, Hinton GE. Variational learning for switching statespace models. Neural Comput. 2000;12(4):831–64. https://doi.org/10.1162/089976600300015619.
 22.
McClintock BT. Incorporating telemetry error into hidden Markov models of animal movement using multiple imputation. J Agric Biol Environ Stat. 2017;22(3):249–69. https://doi.org/10.1007/s1325301702856.
 23.
Zielinski WJ, Truex RL, Schmidt GA, Schlexer FV, Schmidt KN, Barrett RH. Resting habitat selection by fishers in California. J Wildlife Manag. 2004;68(3):475–92. https://doi.org/10.2193/0022541X(2004)068[0475:RHSBFI]2.0.CO;2.
 24.
Aubry KB, Raley CM, Buskirk SW, Zielinski WJ, Schwartz MK, Golightly RT, et al. Metaanalyses of habitat selection by fishers at resting sites in the Pacific coastal region. J Wildlife Manag. 2013;77(5):965–74. https://doi.org/10.1002/jwmg.563.
 25.
Stewart FEC, Darlington S, Volpe JP, McAdie M, Fisher JT. Corridors best facilitate functional connectivity across a protected area network. Sci Rep. 2019;9(1):10852. https://doi.org/10.1038/s4159801947067x.
 26.
Stewart FEC, Fisher JT, Burton AC, Volpe JP. Species occurrence data reflect the magnitude of animal movements better than the proximity of animal space use. Ecosphere. 2018;9(2):e02112. https://doi.org/10.1002/ecs2.2112.
 27.
Blackwell PG, Niu M, Lambert MS, LaPoint SD. Exact Bayesian inference for animal movement in continuous time. O’Hara RB, editor. Methods Ecol Evol. 2016;7(2):184–95. https://doi.org/10.1111/2041210X.12460.
 28.
Nams VO. Combining animal movements and behavioural data to detect behavioural states. Ecol Lett. 2014;17(10):1228–37. https://doi.org/10.1111/ele.12328 Moorcroft P, editor.
 29.
Laliberte AS, Ripple WJ. Range contractions of north American carnivores and ungulates. BioSci. 2004;54(2):123–38. https://doi.org/10.1641/00063568(2004)054[0123:RCONAC]2.0.CO;2.
 30.
Zielinski WJ, Truex RL, Schlexer FV, Campbell LA, Carroll C. Historical and contemporary distributions of carnivores in forests of the Sierra Nevada, California, USA. J Biogeogr. 2005;32(8):1385–407. https://doi.org/10.1111/j.13652699.2005.01234.x.
 31.
Weir RD, Phinney M, Lofroth EC. Big, sick, and rotting: why tree size, damage, and decay are important to fisher reproductive habitat. For Ecol Manag. 2012;265:230–40. https://doi.org/10.1016/j.foreco.2011.10.043.
 32.
Matthews SM, Green DS, Higley JM, Rennie KM, Kelsey CM, Green RE. Reproductive den selection and its consequences for fisher neonates, a cavityobligate mustelid. J Mammal. 2019;100(4):1305–16. https://doi.org/10.1093/jmammal/gyz069.
 33.
Purchell KL, Mazzoni AK, Mori SR, Boroski BB. Resting structures and resting habitat of fishers in the southern sierra Navada, California. Forest Ecol Manag. 2009;258(12):2696–706. https://doi.org/10.1016/j.foreco.2009.09.041.
 34.
Zielinski WJ, Gray AN. 2018. Using routinely collected regional forest inventory data to conclude that resting habitat for the fisher (Pekania pennanti) in California is stable over∼ 20 years. For Ecol Manag. 2018;409:899–908. https://doi.org/10.1016/j.foreco.2017.12.025.
 35.
United States Fish and Wildlife Service. Threatened Species Status for the West Coast District Population Segment of Fisher. FWSR8ES20180105. Fed Reg. 2019;50(CFR 17):60278–305.
 36.
United States Fish and Wildlife Service. Endangered and Threatened Wildlife and Plants; Endangered Species Status for Southern Sierra Nevada Distinct Population Segment of Fisher. Docket No. FWSR8ES20180105, FF09E21000 FXES11110900000 201. Fed Reg. 2020;85:29532–89.
 37.
Oregon Department of Fish and Wildlife. Oregon Conservation Strategy. Salem: Oregon Department of Fish and Wildlife; 2016. Viewed online (March 2020): https://www.oregonconservationstrategy.org/
 38.
Aubry KB, Raley CM, Cunningham PG. Selection of rest structures and microsites by fishers in Oregon. J Wildlife Manag. 2018;82(6):1273–84. https://doi.org/10.1002/jwmg.21479.
 39.
Aubry KB, Lewis JC. Extirpation and reintroduction of fishers (Martes pennanti) in Oregon: implication for their conservation in the Pacific states. Biol Conserv. 2003;114(1):79–90. https://doi.org/10.1016/S00063207(03)00003X.
 40.
Raine RM. Winter habitat use and responses to snow cover of fisher ( Martes pennanti ) and marten ( Martes americana ) in southeastern Manitoba. Can J Zool. 1983;61(1):25–34. https://doi.org/10.1139/z83002.
 41.
Barry BR. Distribution, habitat association, and conservation status of Pacific fisher (Pekania pennanti) in Oregon: MS thesis. Oregon State University; 2018. https://ir.library.oregonstate.edu/concern/graduate_thesis_or_dissertations/2f75rf103.
 42.
Mortenson JA, Moriarty KM. Ketamine and midazolam anesthesia in Pacific martens (Martes caurina). J Wildlife Dis. 2015;51(1):250–4. https://doi.org/10.7589/201402031.
 43.
Hurford A. GPS measurement error gives rise to spurious 180° turning angles and strong directional biases in animal movement data. PLoS ONE. 2009;4:e5632 Rands S, editor.
 44.
Jerde CL, Visscher DR. GPS measurement error influences on movement model parameterization. Ecol Appl. 2005;15(3):806–10. https://doi.org/10.1890/040895.
 45.
Ferrari S, CribariNeto F. Beta regression for modelling rates and proportions. J Appl Stat. 2004;31(7):799–815. https://doi.org/10.1080/0266476042000214501.
 46.
FrühwirthSchnatter S. Data augmentation and dynamic linear models. J Time Ser Anal. 1994;15(2):183–202. https://doi.org/10.1111/j.14679892.1994.tb00184.x.
 47.
Carpenter B, Gelman A, Hoffman MD, Lee D, Goodrich B, Betancourt M, et al. Stan: A Probabilistic Programming Language. J Stat Softw. 2017;76 [cited 2017 Feb 2]. http://www.jstatsoft.org/v76/i01/.
 48.
Monnahan CC, Thorson JT, Branch TA. Faster estimation of Bayesian models in ecology using Hamiltonian Monte Carlo. Methods Ecol Evol. 2017;8(3):339–48. https://doi.org/10.1111/2041210X.12681 O’Hara RB, editor.
 49.
R Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2020. https://www.rproject.org
 50.
Stan Development Team. RStan: the R Interface to Stan. 2020 [cited 2020 Aug 18]. mcstan.org.
 51.
Arnold JB. State Space Models in Stan. 2016 [cited 2020 Aug 20]. https://jrnold.github.io/ssmodelsinstan/index.html.
 52.
Betancourt M. Identifying Bayesian Mixture Models. 2017 [cited 2020 Aug 18]. https://mcstan.org/users/documentation/casestudies/identifying_mixture_models.html.
 53.
McClintock BT. Worth the effort? A practical examination of random effects in hidden Markov models for animal telemetry data. Ecology. 2020. https://doi.org/10.1101/2020.07.10.196410.
 54.
Thieurmel B, Elmarhraoui A. suncalc: compute sun position, sunlight phases, moon position and lunar phase. 2019 [cited 2020 Aug 18]. https://CRAN.Rproject.org/package=suncalc.
 55.
Chacón JE, Duong T. Multivariate Kernel Smoothing and its Applications. 1st ed: Chapman and Hall/CRC; 2018. [cited 2020 Aug 13]. https://www.taylorfrancis.com/books/9780429939143
 56.
Moriarty KM, Kelsey CM, Matthews SM. Assessing den, rest site, and movement characteristics by Pacific fisher (Pekania pennanti) in the southern Oregon cascades: final report: USDA Forest Service Pacific Northwest Research Station; 2019.
 57.
Moriarty KM, Epps CW. Retained satellite information influences performance of GPS devices in a forested ecosystem. Wildlife Soc Bull. 2015;39(2):349–57. https://doi.org/10.1002/wsb.524.
 58.
Powell R. The fisher. Life history, ecology, and behavior. Minneapolis: University of Minnesota Press; 1982.
 59.
Nathan R, Getz WM, Revilla E, Holyoak M, Kadmon R, Saltz D, et al. A movement ecology paradigm for unifying organismal movement research. Proc Nat Acad Sci. 2008;105(49):19052–9. https://doi.org/10.1073/pnas.0800375105.
 60.
Langrock R, King R, Matthiopoulos J, Thomas L, Fortin D, Morales JM. Flexible and practical modeling of animal telemetry data: hidden Markov models and extensions. Ecology. 2012;93(11):2336–42. https://doi.org/10.1890/112241.1.
 61.
Green RE, Purcell KL, Thompson CM, Kelt DA, Wittmer HU. Microsites and structures used by fishers (Pekania pennanti) in the southern Sierra Nevada: a comparison of forest elements used for daily resting relative to reproduction. For Ecol Manag. 2019;440:131–46. https://doi.org/10.1016/j.foreco.2019.02.042.
 62.
Calabrese JM, Fleming CH, Gurarie E. Ctmm: an r package for analyzing animal relocation data as a continuoustime stochastic process. Freckleton R, editor. Methods Ecol Evol. 2016;7(9):1124–32. https://doi.org/10.1111/2041210X.12559.
 63.
Noonan MJ, Fleming CH, Akre TS, DrescherLehman J, Gurarie E, Harrison AL, et al. Scaleinsensitive estimation of speed and distance traveled from animal tracking data. Mov Ecol. 2019;7(1):35. https://doi.org/10.1186/s4046201901771.
 64.
Michelot T, Blackwell PG. Stateswitching continuoustime correlated random walks. Methods Ecol Evol. 2019;10(5):637–49. https://doi.org/10.1111/2041210X.13154.
 65.
Gelman A, Carling JB, Stern HS, Dunson DB, Vehtari A, Rubin D. Bayesian data analysis. Third. Boca Raton: CRC Press; 2014.
 66.
Tweedy PJ. Diel rest structure selection and multiscale analysis of Pacific marten resting habitat in Lassen National Forest, California. Corvallis: Oregon State University; 2018.
 67.
McClintock BT, King R, Thomas L, Matthiopoulos J, McConnell BJ, Morales JM. A general discretetime modeling framework for animal movement using multistate random walks. Ecol Monogr. 2012;82(3):335–49. https://doi.org/10.1890/110326.1.
 68.
Moriarty KM, Linnell MA, Chasco B, Epps CW, Zielinski WJ. 2017. Using highresolution shortterm location data to describe territoriality in Pacific martens. J Mammal. 2017;98(3):679–89. https://doi.org/10.1093/jmammal/gyx014.
Acknowledgments
Field work was under the direction of Drs. Katie Moriarty and Sean Matthews (Oregon State University). Survey efforts were funded by the USDOI Bureau of Land Management (BLM) Portland Office and Lakeview District. Additional funding and support were provided by the USDI Fish and Wildlife Service (USFWS) Klamath Falls Office; USDA Forest Service Fremont Winema National Forest and Region 6; NCASI; Hancock Timber Resource Group; and Weyerhaeuser Company. Thanks to Elizabeth Willy, Jenniffer Bakke, Sue Livingston, and Jessica Homyack for support and grant facilitation. Considerable aid with field logistics, vehicles, BLM staff, and equipment were provided by the Klamath Falls BLM Resource Area Lakeview District (Steve Hayner, Matt Broyles). Jeffery Stephens and the Medford BLM office helped during project initiation as well as with logistics. Thanks to team members Caylen Kelsey, Sean Anderson, Alyssa Roddy, Craig Heinrich, Isaac Kelsey, Todd Musser.
Author information
Affiliations
Contributions
KMM initiated this study, collected field data, and established research objectives. DJH developed statistical models and conducted all statistical analyses. DJH, KMM, BH, and RWP wrote the manuscript. All authors approved of final manuscript submission.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
We captured and processed fishers using methods approved by the USDA Forest Service’s Institute for Animal Care and Use Committee (Permit 2015–003) and Oregon Department of Fish and Wildlife Scientific Take permit (Permits 034–16, 027–17).
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.
The original version of this article was revised: the omitted file for Additional file 2 was added.
Supplementary Information
Additional file 1: Supplemental Figure 1.
Comparison of observed activity (left two panels) and observed apparent steplengths (rightpanel) to estimated statedependent beta and lognormal probability distribution, respectively. The probability distributions lines consist of 100 independent draws from the joint. Supplemental Figure 2. Duration of rest events (hours) for 21 GPS collar deployments on 9 individual fishers. A rest event is defined as continuous segments of time series for which P(s_{t} = 1) > 0.95. Supplemental Figure 3. Cumulative distance moved per day for 21 GPS collar deployments of 9 individual fishers. Points depict the posterior mean cumulative distance moved per day. Error bars on the points depict the 5^{th} and 95^{th} percentile of the posterior distribution. Boxplots are based off the posterior mean cumulative distance moved between each rest event and depict median (center line), first and third quartile (lower and upper hinge), with whiskers extending to 5^{th} and 95^{th} percentile of the posterior mean. Supplemental Figure 4a. Example comparisons from field validation of HMM/SSM identified rest events and contemporaneous radio telemetry identified rest events. Of 67 suspected rest events radio telemetry resting locations which could be matched in time to a HMM/SSM identified rest events, 8 were matched to rest events with no successful GPS fixes and large credible ellipses. Four of these are depicted if Figure 4a including the only example where the radio telemetry rest event was outside of the 95% credible ellipse for the HMM/SSM resting location (bottom right panel). The yellow target symbol is the radio telemetry identified resting location, and the white contours represent the 50%, 80%, and 95% credible ellipses for the HMM/SSM resting location. Supplemental Figure 4b. Example comparisons from field validation of HMM/SSM identified rest events and contemporaneous radio telemetry identified rest events. Of 67 suspected rest events radio telemetry resting locations which could be matched in time to a HMM/SSM identified rest events, 22 lay within the 95% credible ellipse estimated from the HMM/SSM. Four of these are depicted in Figure 4b. The yellow target symbol is the radio telemetry identified resting location, the red points are all GPS location fixes associated with the rest event, and the white contours represent the 50%, 80%, and 95% credible ellipses for the HMM/SSM resting location. Supplemental Figure 4c. Example comparisons from field validation of HMM/SSM identified rest events and contemporaneous radio telemetry identified rest events. Of 67 suspected rest events radio telemetry resting locations which could be matched in time to a HMM/SSM identified rest events, 27 were outside of the 95% credible ellipse but within 25m of the nearest edge. Four of these are depicted in Figure 4c. The yellow target symbol is the radio telemetry identified resting location, the red points are all GPS location fixes associated with the rest event, and the white contours represent the 50%, 80%, and 95% credible ellipses for the HMM/SSM resting location. Supplemental Figure 4d. Example comparisons from field validation of HMM/SSM identified rest events and contemporaneous radio telemetry identified rest events. Of 67 suspected rest events radio telemetry resting locations which could be matched in time to a HMM/SSM identified rest events, 10 were located further then 25 meters from the nearest edge of the corresponding HMM/SSM 95% credible ellipse. Four of these are depicted in Figure 4d including the example furthest from the closet edge where the radio telemetry rest event was identified 433 meters away from the 95% credible ellipse for the HMM/SSM resting location (bottom right panel). The yellow target symbol is the radio telemetry identified resting location, and the white contours represent the 50%, 80%, and 95% credible ellipses for the HMM/SSM resting location.
Additional file 2:
Complete R and Stan code for reproducing the statistical analyses in this manuscript.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Hance, D.J., Moriarty, K.M., Hollen, B.A. et al. Identifying resting locations of a small elusive forest carnivore using a twostage model accounting for GPS measurement error and hidden behavioral states. Mov Ecol 9, 17 (2021). https://doi.org/10.1186/s40462021002568
Received:
Accepted:
Published: