Flexible hidden Markov models for behaviour-dependent habitat selection

Background There is strong incentive to model behaviour-dependent habitat selection, as this can help delineate critical habitats for important life processes and reduce bias in model parameters. For this purpose, a two-stage modelling approach is often taken: (i) classify behaviours with a hidden Markov model (HMM), and (ii) fit a step selection function (SSF) to each subset of data. However, this approach does not properly account for the uncertainty in behavioural classification, nor does it allow states to depend on habitat selection. An alternative approach is to estimate both state switching and habitat selection in a single, integrated model called an HMM-SSF. Methods We build on this recent methodological work to make the HMM-SSF approach more efficient and general. We focus on writing the model as an HMM where the observation process is defined by an SSF, such that well-known inferential methods for HMMs can be used directly for parameter estimation and state classification. We extend the model to include covariates on the HMM transition probabilities, allowing for inferences into the temporal and individual-specific drivers of state switching. We demonstrate the method through an illustrative example of plains zebra (Equus quagga), including state estimation, and simulations to estimate a utilisation distribution. Results In the zebra analysis, we identified two behavioural states, with clearly distinct patterns of movement and habitat selection (“encamped” and “exploratory”). In particular, although the zebra tended to prefer areas higher in grassland across both behavioural states, this selection was much stronger in the fast, directed exploratory state. We also found a clear diel cycle in behaviour, which indicated that zebras were more likely to be exploring in the morning and encamped in the evening. Conclusions This method can be used to analyse behaviour-specific habitat selection in a wide range of species and systems. A large suite of statistical extensions and tools developed for HMMs and SSFs can be applied directly to this integrated model, making it a very versatile framework to jointly learn about animal behaviour, habitat selection, and space use. Supplementary Information The online version contains supplementary material available at 10.1186/s40462-023-00392-3.

Here, we derive the SSF covariates needed to model step lengths with a gamma distribution.
Following Rhodes et al. (2005) and Forester et al. (2009), the two-dimensional spatial distribution of the endpoint y of a step (given that it started at x) is related to the distribution of the step length L = ||y − x|| as follows This defines the relationship between the two-dimensional distribution p(y | x) and the onedimensional distribution of the distance L between the start point x and end point y. Therefore, the distribution of y for gamma-distributed step lengths with shape a and scale b is obtained by replacing p(L) by the gamma probability density function, ( 2) where Γ is the gamma function. To model this in the SSF, we write the movement kernel as the exponential of a linear predictor, such that ϕ(y | x) = exp(c(x, y) · β). Therefore, to identify the terms of linear predictor, we can take the natural log of the movement kernel, Note that the last term, log(b a Γ(a)2π), is a constant (i.e., not dependent on the step length) and can be ignored. Finally, this leaves two terms that need to be included in the linear predictor of the SSF: step length L and its logarithm log(L). The equation gives us the relationship between the gamma distribution parameters and the selection coefficients (here, referred to as β 1 and β 2 for L and log(L), respectively): β 1 = −1/b and β 2 = a − 2. In turn, we can derive mean µ and standard deviation σ of the step length distribution, in terms of the selection parameters. For the gamma distribution, it is known that from which we can derive the following Here, we describe the approach for for the gamma distribution, but this could be applied to other distributions from the exponential family suitable to model step lengths (Forester et al., 2009;Avgar et al., 2016). For example, the exponential distribution is a special case of the gamma distribution, where a = 1.

A.2 Turning angle distribution
Similarly to the procedure used above for step length, we can derive covariates required to model a given distribution of turning angle. Following Avgar et al. (2016) and Nicosia et al. (2017), we focus our attention on the von Mises distribution, with probability density function where θ ∈ (−π, π] is the turning angle, µ ∈ (−π, π] is the mean parameter, κ > 0 is the concentration parameter, and I 0 is the modified Bessel function of the first kind of order 0. Intuitively, κ is inversely related to the variance of the distribution, and large κ corresponds to a high peak around µ. We first focus on the case where µ = 0, corresponding to a tendency to persist in direction, which is the most prevalent in the context of animal movement (but we discuss the case µ = π below). Assuming µ = 0, we take the log of the density function to identify terms that should be included in the SSF, The term log[2πI 0 (κ)] does not depend on turning angle, and we can therefore omit it. Finally, we get the result that κ cos(θ) should be added to the SSF linear predictor to model turning angle with a von Mises distribution with mean 0. In practice, this means that the cosine of turning angle is included as a covariate in the model, and the corresponding selection parameter is equal to the concentration parameter of the distribution.
The above assumed that the selection parameter β θ associated with cos(θ) was strictly positive (because κ > 0 and β θ = κ). We can relax this assumption by noticing that This implies that, if β θ is negative, we can interpret −β θ as the concentration parameter of a von Mises distribution centred on π. A turning angle distribution centred on π corresponds to frequent reversals in direction, which is commonly observed as an artifact of discrete-time observations when animals are inactive.
In summary, including cos(θ) as a covariate in the SSF makes it possible to model turning angle with a von Mises distribution, either centred on 0 (if the corresponding selection coefficient β θ is positive) or centred on π (if β θ is negative). This covers the vast majority of animal movement scenarios, as other values of the mean turning angle are virtually never used.

B Simulation study
We conducted simulations to verify that our implementation method is able to return the correct parameter values. The basic structure and objective of the simulation study was to: i) simulate a movement track from the HMM-SSF (i.e., with known parameter values), and ii) fit the HMM-SSF with the implementation method described in Section 2.4 and check how well the parameters were estimated.

B.1 Methods
For all simulations, we followed the same algorithm to produce each movement track: 1. Simulate a state sequence {S 1 , S 2 , ..., S T }, with S 1 sample with probabilities given by δ and S 2 , ..., S T based on Γ 2. Generate starting location y 1 randomly within the study area Ω 3. To generate each y t in {y 2 , y 3 , ..., y T } (a) Generate many possible endpoints {z 1 , z 2 , ..., z J } uniformly on a disc with radius r, centred on y t and evaluate/interpolate relevant covariates.
(b) Select y t+1 from the possible endpoints {z 1 , z 2 , ..., z J } for j ∈ 1, 2, ..., J with probabilities given by the state-dependent SSF (for the known state S t = k), i.e., In practice, we generated J = 10, 000 endpoints within a disc of radius r = 5. J was chosen to be very large to reduce the risk of bias arising from the simulation itself. The initial location y 1 was generated from the middle 50% of Ω. Simulations were run for 100 iterations, each for T = 750 observations.
We assumed a two-state model meant to represent a generally ideal dataset (i.e., with distinct states and moderate effect sizes), while still being realistic. We assumed that the step lengths followed a gamma distribution (i.e., we included step length and its natural log as covariates) with state-specific means µ (1) = 0.5 and µ (2) = 2 and standard deviations σ (1) = 0.3 and σ (2) = 1. We assumed that turning angles followed a von Mises distribution with a mean of 0 (i.e., we included the cosine of the turning angle as a covariate) and state-specific angular concentration parameters κ (1) = 0.25 and κ (2) = 5. These parameters can be translated to the corresponding distribution and β parameters following Appendix A. We also considered selection for two additional habitat covariates, which were simulated using a moving average window with spatial autocorrelation parameter ρ (Avgar et al., 2016;Klappstein et al., 2022). Both covariate rasters were simulated with dimensions of 2000 × 2000 with a resolution of 1. The general raster simulation process was to: i) generate a random value for each raster cell from [U ∼ (0, 1)], and ii) determine the covariate value based on a circular moving average window with radius ρ (i.e., a higher ρ indicates more spatial autocorrelation). We defined ρ 1 = 5 for covariate 1 and ρ 2 = 25 for covariate 2 (Figure 1). In state 1, there was selection for covariate 1 (β = 3) and selection against covariate 2 (β = −1). This selection pattern was reversed in state 2, with selection against covariate 1 (β = −1) and selection for covariate 2 (β = 3).
We additionally included covariates on the transition probabilities (in addition to the SSF covariates, previously described). We simulated data at a 1-hour time resolution (i.e., 24 steps in a day). Following Towner et al. (2016), we included time of day as a cyclic covariate, such that the linear predictor in Equation 8 of the main text becomes where τ t is the time of day at time t. We set our simulation parameters for the two-state model as, to represent a higher probability of switching into state 1 during the day and a higher probability

B.2 Results
Generally, movement and habitat selection parameters were estimated well, with the exception of the estimate for covariate 2 in state 1 (Figures 2 and 3). This is most likely because there is high spatial autocorrelation in covariate 2 and short step lengths in state 1. This seems to cause bias in the direction of the (positive) selection in state 2 because, when the animal switches into state 1, it will tend to already be in an area of high values of covariate 2, and will not be able to move to low values before switching again. The general pattern of transition probabilities in relation to time of day was also captured, but with higher variance on the transition from state 2 to state 1 ( Figure 4). States were decoded correctly 98.2% of the time. (1-2) and from state 2 to 1 (2-1). Each line is one of the 100 simulation iteration estimates and the black dashed line is the true relationship.