A hidden Markov movement model for rapidly identifying behavioral states from animal tracks

Abstract Electronic telemetry is frequently used to document animal movement through time. Methods that can identify underlying behaviors driving specific movement patterns can help us understand how and why animals use available space, thereby aiding conservation and management efforts. For aquatic animal tracking data with significant measurement error, a Bayesian state‐space model called the first‐Difference Correlated Random Walk with Switching (DCRWS) has often been used for this purpose. However, for aquatic animals, highly accurate tracking data are now becoming more common. We developed a new hidden Markov model (HMM) for identifying behavioral states from animal tracks with negligible error, called the hidden Markov movement model (HMMM). We implemented as the basis for the HMMM the process equation of the DCRWS, but we used the method of maximum likelihood and the R package TMB for rapid model fitting. The HMMM was compared to a modified version of the DCRWS for highly accurate tracks, the DCRWSNOME, and to a common HMM for animal tracks fitted with the R package moveHMM. We show that the HMMM is both accurate and suitable for multiple species by fitting it to real tracks from a grey seal, lake trout, and blue shark, as well as to simulated data. The HMMM is a fast and reliable tool for making meaningful inference from animal movement data that is ideally suited for ecologists who want to use the popular DCRWS implementation and have highly accurate tracking data. It additionally provides a groundwork for development of more complex modeling of animal movement with TMB. To facilitate its uptake, we make it available through the R package swim.


| 2113
WHORISKEY Et al. knowledge informs the management and conservation of both species and ecosystems. In aquatic environments, where direct observation of animal movement and behavioral states is often impossible, researchers are rapidly expanding the use of electronic telemetry for documenting animal movement through time (Hussey et al., 2015).
Satellite telemetry and acoustic positioning systems are common types of telemetry technology for estimating an aquatic animal's location in continuous space and yield time series of locations along an animal's path, usually referred to as tracks. Inferring behavioral states from animal tracks is possible by assuming that different types of movement, and therefore behavioral states, can be reflected by changes in characteristics of an animal's path. For example, while foraging can often be characterized by a tortuous track, a more directed path may suggest traveling between habitats (e.g., Jonsen, Mills Flemming, & Myers, 2005).
Hidden Markov models (HMMs) are a popular tool used to identify behavioral states from animal telemetry data with negligible error (e.g., Morales, Haydon, Frair, Holsinger, & Fryxell, 2004;Langrock, King, Matthiopoulos, Thomas, Fortin, & Morales, 2012). HMMs are a large class of models distinguished in the most general case by a set of observations that depend on an unobserved, underlying Markov process (Zucchini, MacDonald, & Langrock, 2016). In the context of animal movement, the latent Markov process is used to model the discrete behavioral states of interest, while the set of observations follow a movement process that can also be Markovian. The observations can consist of either location data (e.g., Jonsen et al., 2005) or metrics derived from the observed track, such as turning angles and step lengths (e.g., Morales et al., 2004;Langrock et al., 2012).
While current HMMs can be fitted rapidly using maximum-likelihood (ML) methods, with the exception of the formulation of Pedersen, Patterson, Thygesen and Madsen (2011), they are unable to account for measurement error associated with the technology used to obtain animal tracks. For those tracks measured with error, state-space models (SSMs) provide a more accurate and reliable method for identifying behavioral states, but are typically fitted using comparatively slow Bayesian methods such as Markov chain Monte Carlo (MCMC) sampling because of large numbers of random effects (e.g., Jonsen et al., 2005;McClintock, King, Thomas, Matthiopoulos, McConnell, & Morales, 2012;Jonsen, 2016). Therefore, the ideal tool for identifying behavioral states from animal tracks should incorporate features of both HMM and SSM implementations, such that measurement error can be accounted for within a ML framework that keeps the computational burden of estimation relatively small.
Here, we introduce a new HMM for estimating behavioral states from highly accurate animal tracks which is similar to the DCRWS, but does not account for measurement error. We directly implement the process equation of the DCRWS as the basis for our HMM, but we adjust the model for fitting within a ML framework to allow for rapid estimation. Model fitting and parameter estimation are performed using the R package TMB (Kristensen, Nielsen, Berg, & Skaug, 2016), which has previously shown great promise for analyzing animal tracking data (Albertsen, Whoriskey, Yurkowski, Nielsen, & Mills Flemming, 2015;Auger-Méthé et al. 2017). We make this model, which we entitle the hidden Markov movement model (HMMM), available through the R package swim (see supplementary material). To demonstrate the accuracy and applicability of the HMMM, we apply it to simulated animal tracks and to real tracks from multiple aquatic species.
We additionally compare our HMMM results to those obtained using its Bayesian counterpart and to results from the moveHMM package (Michelot, Langrock, & Patterson, 2016). We assess the advantages and disadvantages of each approach by comparing their computational efficiency, accuracy, and sequences of behavioral states.

| The DCRWS movement process
The DCRWS is a SSM that estimates the true locations, behavioral states, and parameters of a movement process from an Argos satellite system track (Jonsen et al., 2005). Given the true location x t at time t, the process equation of the DCRWS is a correlated random walk on the first differences of the true locations, d t = x t − x t−1 : The stochastic term in the movement process is a bivariate Gaussian (N 2 ) with mean 0 and covariance matrix Σ, where σ lat and σ lon are the standard deviations in the latitude and longitude axes, respectively, and ρ is the correlation between the two axes. Like Breed et al. , and these parameter values are dependent on the behavioral state for distinguishing between multiple behavioral states at each location. Typically of interest are two states: the first is directed movement characterized by traveling in the same direction (θ ≈ 0) and at a similar, high speed (γ > 0.5), and the second is tortuous movement characterized by frequent course reversals (θ ≈ π) at dissimilar, slower speeds (γ < 0.5). Parameter sets for each state are identified with the appropriate subscript, either 1 or 2.

| The HMMM
Our HMMM uses the movement process described by (1), but instead of using a Bayesian framework like Jonsen et al. (2005), we employ a ML framework and fit the process equation as a HMM via TMB, which requires that the likelihood function be coded in the C++ programming language. The probability distribution of the movement process is conditional on the assumed behavioral states b t and is given by The likelihood of the HMMM is that of a HMM (Zucchini et al., 2016): Assuming two behavioral states, the 2×1 vector δ contains the initial probabilities of being in each state. A is a 2×2 transition probability matrix containing the switching probabilities α i,j that describe the probability of switching from state i at time t−1 to state j at time t.
Because the rows of A sum to 1, we need only estimate two switching probabilities instead of four; we choose to estimate α 1,1 and α 2,1 . P is a 2×2 diagonal matrix with diagonal entries equal to (f(d t |b t−1 = 1), f(d t |b t−1 = 2)), that is, the probabilities of being at the observed locations given each behavioral state as described by the movement process. 1 is a 2×1 vector of ones. We estimate the parameters of the movement process directly from the likelihood within TMB and then use the Viterbi algorithm to estimate the unobserved behavioral states (Zucchini et al., 2016).

| Data analysis and simulation study
To evaluate the performance of the HMMM, we compared it with two other approaches for estimating behavioral states from animal tracks with negligible measurement error. The first was the switching movement process described by (1), fitted using a Bayesian framework and Markov chain Monte Carlo (MCMC) sampling via rjags (Plummer, 2015). This first model is the DCRWS of Jonsen et al. (2005) modified such that it does not include a measurement equation, and therefore differs from the HMMM solely in implementation (i.e., Bayesian vs. ML inference). Hereafter, we refer to it as the DCRWS NOME . Although the DCRWS NOME has not been fitted before, implementations of the DCRWS for tracking data with minimal errors do exist (Jonsen, 2016), and the DCRWS NOME is the most direct implementation of the DCRWS when no measurement error is assumed. To fit the DCRWS NOME , we used a burn-in period of 40,000 samples and then sampled 20,000 from the posterior distribution but only kept every 20th sample (thinning). We fitted and compared two MCMC chains to each track to check for convergence. All prior distributions were specified as in the R package bsam (Jonsen et al., 2005) that fits the original DCRWS model, with the exception of those for the error covariance matrix Σ.
Instead, by setting ρ = 0, which we believe is more appropriate, we were able to specify separate vague uniform priors on σ lon and σ lat as opposed to using the original Wishart prior on the entire matrix (Jonsen et al., 2005). Parameters and behavioral states were estimated as the posterior medians of the samples from the two chains combined. We additionally fitted a HMM to the turning angles (rad) and step lengths (km) of the animal tracks with the R package moveHMM (Michelot et al., 2016), using a von Mises (mean μ and concentration parameter c) and Weibull (shape λ and scale parameter k) distribution, respectively. Behavioral states were again identified via the Viterbi algorithm, using functions from moveHMM.
We fitted these three models to three animal tracks: (1)  Because geolocation data can be error-prone, the track was processed with the Wildlife Computers (2015) GPE3 software (a SSM) to improve positioning accuracy by estimating the true blue shark locations, as suggested by the manufacturers. Although a SSM for geolocation data would have been the ideal approach for analyzing the blue shark data, our approach of fitting HMMs to true location SSM estimates has been previously adopted (e.g., Eckert et al., 2008). Because the data were collected in continuous-time but all three models assume underlying discrete-time Markov processes, we had to approximate the locations in discrete-time and then assume that these were known.
We linearly interpolated the data sets over time using a 6-hr, 15-min, and 12-hr time step for the grey seal, lake trout, and blue shark data, respectively, yielding data sets with 1,227, 2,187, and 393 locations.
Different time steps were required based on the different temporal resolutions of the tracks.
Additionally, using the parameter estimates from the grey seal HMMM and DCRWS NOME fits (Table 1), we conducted a simulation study to formally compare the accuracy of the HMMM and DCRWS NOME and compare their results with those obtained using moveHMM. We simulated 50 tracks corresponding to the HMMM from a known parameter set with turning angles θ 1 = 0 and θ 2 = π; autocorrelation γ 1 = 0.8 and γ 2 = 0.05; process error standard deviations σ lon = 0.07 and σ lat = 0.05; and switching probabilities α 1,1 = 0.89 and α 2,1 = 0.20. These two behavioral states are classically interpreted as transiting (θ 1 , γ 1 ) and foraging (θ 2 , γ 2 ). We then fitted the HMMM, moveHMM, and the DCRWS NOME to each simulated track and calculated the parameter estimates and interval measures of uncertainty for these estimates. For the HMMM and moveHMM, this consisted of the 95% confidence interval based on the standard error estimates. For the DCRWS NOME , we determined the 95% credible interval as the 2.5% and 97.5% quantiles of the posterior samples. We found the behavioral state error rate, that is, the proportion of states that were incorrectly identified, for each model, and additionally calculated the root mean squared error (RMSE) for each parameter estimate Θ from the HMMM and DCRWS NOME fits as We were unable to calculate the RMSE for the moveHMM fits because the data were simulated according to the HMMM movement process and the moveHMM implementation does not involve the same parameters.

| Identifying behavioral states
We applied the HMMM, the DCRWS NOME , and moveHMM to a grey seal, lake trout, and blue shark track estimated with negligible measurement error. All three models performed similarly and identified two clearly distinct behavioral states for the grey seal and lake trout tracks (Figures 1 and 2). For both animals, HMMM and DCRWS NOME parameter estimates were similar except for the tortuous turning angle θ 2 , which was estimated at a similar distance from the number π but in the opposite direction (i.e., while one was estimated turning slightly to the left, the other was estimated turning to the right; Tables 1 and   2). This, along with the relatively large confidence intervals for θ 2 , is not unusual because together with small γ 2 , it suggests that the animal is exhibiting tortuous movement, in which case the mean turning angle does not have as much influence because the animal is more equally likely to travel in any direction. Switching probabilities (α 1,1 and α 2,1 ) were similar among all three models for the seal track.
moveHMM estimated switching probabilities for the lake trout track different from the HMMM and DCRWS NOME , although the estimated parameters of the three models led to similar decoded behavioral state sequences. For the seal track, the DCRWS NOME took 6.4 hr to fit, moveHMM took 0.9 s, and the HMMM took 0.06 s. For the lake trout data, the DCRWS NOME took 9.4 hr to fit, moveHMM took 1.8 s, and the HMMM took 0.16 s.
All three models identified two states from the blue shark track, although half of the switching probabilities estimated by moveHMM differed greatly from those estimated by the HMMM and DCRWS NOME (Table 3), and this led to different state sequences (Figure 3). Specifically, all three models estimated a high probability of remaining in state 1, α 1,1 , but moveHMM estimated a low probability of switching from state 2 to state 1, α 2,1 , while the DCRWS NOME and HMMM estimated a high α 2,1 . The switching probabilities of the HMMM and DCRWS NOME therefore led to state sequences containing long stretches of state 1 interspersed with short (length 1 or 2) stretches of state 2. By contrast, moveHMM estimated a state sequence with longer stretches of both behavioral states. While the DCRWS NOME took 1.7 hr to fit to the blue shark track, moveHMM took 1.2 s, and the HMMM took 0.02 s.

| Simulation Study
We simulated 50 tracks from the HMMM with a set of parameters representative of the grey seal track. The HMMM and DCRWS NOME provided accurate estimates of the model parameters (Figure 4), but the DCRWS NOME had a smaller average (over the parameters) RMSE (0.120 vs. 0.140; Table 4). The RMSEs for individual parameters were similar (within 0.01) between the two models with the exception of θ 2 , where the RMSE of the DCRWS NOME was smaller by 0.149 (Table 4).
The DCRWS NOME additionally had the smallest behavioral state error rate (0.175) which differed from the HMMM and moveHMM by approximately 1.5% (0.189) and 18.7% (0.362), respectively. Finally, the average time needed to fit the DCRWS NOME was 5.10 hr, while moveHMM took 1.2 s and the HMMM took 0.08 s. T A B L E 1 Parameter estimates from three models fitted to a grey seal track. The Lower and Upper columns are the lower and upper bounds of 95% uncertainty intervals around the estimates. These correspond to 95% confidence intervals for the HMMM and moveHMM, and 95% credible intervals for the DCRWS NOME . The only two parameters in common among all three models are the switching probabilities, α 1,1 and α 2,1 T A B L E 2 Parameter estimates from three models fitted to a lake trout track. The Lower and Upper columns are the lower and upper bounds of 95% uncertainty intervals around the estimates. These correspond to 95% confidence intervals for the HMMM and moveHMM, and 95% credible intervals for the DCRWS NOME . The only two parameters in common among all three models are the switching probabilities, α 1,1 and α 2,1 T A B L E 3 Parameter estimates from three models fitted to a blue shark track. The Lower and Upper columns are the lower and upper bounds of 95% uncertainty intervals around the estimates. These correspond to 95% confidence intervals for the HMMM and moveHMM, and 95% credible intervals for the DCRWS NOME . The only two parameters in common among all three models are the switching probabilities, α 1,1 and α 2,1 . Because this track had some step lengths equal to zero, the two parameters ζ 1 and ζ 2 were used to estimate zero-inflation for each behavior when using moveHMM With the HMMM, DCRWS NOME , and moveHMM, we identified two behavioral states within the lake trout track. The Drummond Island lake trout population spawn primarily at nighttime on rock rubble reefs in association with submerged drumlins (Riley et al., 2014;Binder et al., 2015). Lake trout show multiple behaviors characterized by tortuous movement, including spawning on the reefs. For example, lake trout (particularly males) often aggregate on the spawning reefs in the weeks leading up to spawning, a behavior known as staging (Muir, Blackie, Marsden, & Krueger, 2012). Because egg surveys have verified that no spawning occurs in some locations where our models identified tortuous behavior (T. Binder, unpublished observations), we believe the models are distinguishing, more generally, reef and nonreef behaviors. Being able to mathematically distinguish between reef and non-reef behaviors can allow for identification of key lake trout habitats for conservation such as spawning sites in places where direct observation is difficult. Furthermore, by building a dependence of the HMMM on one or more covariates, it may be possible to more acutely identify spawning behavior. For example, because the Drummond Island lake trout tend to spawn at night close to the substrate, time of day and lake trout depth (which is often recorded by positioning systems such as the VPS) may provide sufficient additional information for the HMMM to distinguish spawning behavior from other reef-associated behaviors. One possible way to achieve this is by allowing the switching probabilities of the HMMM to depend on these covariates in a linear fashion (as in, e.g., Bestley, Jonsen, Hindell, Guinet, & Charrassin, 2012;Michelot et al., 2016). Additionally, an extension to the HMMM which could estimate more than two behavioral states may be able to distinguish reef from spawning behavior. We chose to model only two states so that we could more directly compare results of the HMMM to our implementation of the original DCRWS (the DCRWS NOME ); however, the HMMM should be directly extendible.
When fitted to the blue shark track, moveHMM produced different state sequences than the HMMM and DCRWS NOME , as moveHMM estimated longer stretches of behavioral state 2 than either of the other models. This is likely because moveHMM models the distributions of the turning angles and step lengths calculated from an animal path, which is fundamentally different from the movement process of the HMMM and DCRWS NOME . Furthermore, McClintock et al. (2014) showed that the continuous-time analog to the movement process introduced by Jonsen et al. (2005)  Our simulation study results suggested that while the DCRWS NOME was slightly more accurate than the HMMM, the difference was marginal. The two models performed similarly while estimating model parameters with the exception of θ 2 , which the DCRWS NOME more accurately estimated. This result is likely explained by the rather informative priors on θ 1 and γ 1 when fitting the DCRWS NOME . The DCRWS NOME also more accurately estimated the behavioral states, an unsurprising result because while the DCRWS NOME directly estimates these random effects from the posterior likelihood, the HMMM uses a post hoc global decoding algorithm (the Viterbi algorithm) to identify the most likely sequence of states given the ML parameter estimates.
Predictably, moveHMM had the highest behavioral state error rate of the three approaches, likely because it was fitted to simulated data from a movement process not equivalent to its own. Finally, the HMMM was the fastest model to fit, with moveHMM and the DCRWS NOME taking on average 15 times and 229,500 times longer to fit than the HMMM, respectively. Quicker fits of the DCRWS NOME may be achieved by reducing burn-in and sampling sizes of the MCMC, but they would still take Our HMMM is a major advance in using TMB to solve animal movement problems. Highly accurate data are becoming more common in the marine realm, and the HMMM, as implemented through the R package swim, provides a fast and reliable tool for making meaningful inference from animal movement data. Fast methods for analyzing data will become more important as larger data sets are collected. The HMMM therefore additionally provides a baseline method for movement modeling in TMB that can be further developed for more specific and nontrivial animal movement problems such as determining relationships between movement and environmental covariates, or accounting for measurement error.

ACKNOWLEDGMENTS
The grey seal and blue shark tracks were obtained under the approval of animal care protocols 12-64 and 13-002, respectively, from the University Committee on Laboratory Animals, of Dalhousie This paper is contribution number 25 of the Great Lakes Acoustic Telemetry Observation System (GLATOS).
T A B L E 4 Parameter results from the simulation study (n = 50) comparing the HMMM to the DCRWS NOME . The Lower and Upper columns correspond to the 95% confidence and credible intervals for the HMMM and the DCRWS NOME , respectively.