How to account for behavioral states in step-selection analysis: a model comparison

Hidden Markov step-selection model

We use {x_0,1, x_0,2, …, x_0,T} to denote the sequence of two-dimensional animal locations observed at regular time intervals, which forms the observed movement track. We here use the subscript “0” to distinguish these used locations from control locations (introduced below). Conditional on the previous location x_0,t−1, a step from the current location x_0,t to the next location x_0,t+1, is characterized by its step length l_0,t+1, i.e., the straight-line distance between the two consecutive locations, and its turning angle α_0,t+1, i.e., the directional change. The corresponding covariate vector Z_0,t+1 stores the feature values of the step, and we use Z to denote the collection of covariate values for all possible locations in the given spatial domain D in which the animal moves.

In the hidden Markov step-selection model (Nicosia et al., 2017; Klappstein, Thomas & Michelot, 2023; Prima et al., 2022), we assume the observed steps to be driven by an underlying hidden state sequence {S₁, S₂, …, S_T} with N discrete states. Thus, each state variable S_t at time t can take one of N state values (S_t ∈ {1, …, N}). These states serve as proxies for the unobserved behavioral modes of the animal that influence its movement and habitat selection (illustrated in Fig. 1). We assume the state sequence to be a homogeneous N-state Markov chain, characterized by its transition probabilities γ_ij = Pr(S_t = j∣S_t−1 = i) to switch from state i to state j, summarized in the N × N transition probability matrix Γ, and the initial state distribution δ which contains the probabilities to start in a certain state.

Each state i (i = 1, …, N) is associated to a state-dependent density f_i generating the next location. Its functional form is similar to the basic step-selection model (Forester, Im & Rathouz, 2009), but with movement and habitat selection parameters being state-dependent. Thus, conditional on locations x_0,t−1 and x_0,t, and covariates Z, the current state S_t = i determines the following distribution for a step to location x_0,t+1: (1)${\displaystyle f_i\left({\mathbf{x}}_{0,t+1}\mid {\mathbf{x}}_{0,t},{\mathbf{x}}_{0,t-1},\mathbf{Z};{\theta }_i,{\beta }_i\right)=\frac{\stackrel{\text{movement kernel}}{\overbrace{\varphi \left({\mathbf{x}}_{0,t+1}\mid {\mathbf{x}}_{0,t},{\mathbf{x}}_{0,t-1};{\theta }_i\right)}}\cdot \stackrel{\text{selection function}}{\overbrace{\omega \left({\mathbf{Z}}_{0,t+1};{\beta }_i\right)}}}{\underset{\text{normalizing constant}}{\underbrace{{\int }_{\tilde{\mathbf{x}}\in D}\varphi \left(\tilde{\mathbf{x}}\mid {\mathbf{x}}_{0,t},{\mathbf{x}}_{0,t-1};{\theta }_i\right)\cdot \omega \left(\tilde{\mathbf{Z}};{\beta }_i\right)d\tilde{\mathbf{x}}}}},}$ where $\tilde{\mathbf{Z}}$ in the integral denotes the covariate vector of location $\tilde{\mathbf{x}}$. The density consists of three components: (i) The movement kernel ϕ(⋅) describes the space use in a homogeneous landscape and is usually defined in terms of step length l_0,t+1 (e.g., gamma distribution) and turning angle α_0,t+1 (e.g., von Mises distribution). The corresponding state-dependent parameters for state i are summarized in the movement parameter vector θ_i; (ii) The movement kernel is weighted by the selection function ω(⋅) which indicates a possible selection for or against the covariates in Z_0,t+1. It is usually assumed to be a log-linear function of the (state-dependent) selection coefficient vector β_i,

${\displaystyle \omega \left({\mathbf{Z}}_{0,t+1};{\beta }_i\right)=exp\left({\mathbf{Z}}_{0,t+1}^{\top }{\beta }_i\right),}$ where a positive selection coefficient indicates preference for, and a negative coefficient avoidance of a corresponding covariate; (iii) The integral in the denominator ensures that f_i integrates to one. Usually, it is analytically intractable and must therefore be approximated, for example, using numerical integration methods. We provide an example of state-dependent step-selection densities in a 2-state scenario in Fig. S1.

There are important relations between the hidden Markov step-selection model and two alternative movement models: (i) If all states share the same parameters, i.e.,  θ₁ = … = θ_N and β₁ = … = β_N, or if the number of states is set to one, i.e.,  N = 1, the model reduces to the basic step-selection model without state-switching (Forester, Im & Rathouz, 2009); (ii) If all selection coefficients are equal to zero, i.e.,  β₁ = … = β_N = 0, the model reduces to a basic movement HMM(Langrock et al., 2012, Patterson et al., 2017) with state-dependent step length and turning angle distributions as implied by the movement kernel ϕ(⋅) but without habitat selection. These relations are very convenient, for example in the context of model comparison and model selection, as it allows the use of standard tests or information criteria to select between these three candidate models.

We can simplify the step-selection density f_i by assuming step lengths follow a distribution from the exponential family (with support on non-negative real numbers, e.g., a gamma distribution) and turning angles follow either a uniform or von Mises distribution with fixed mean (Avgar et al., 2016; Nicosia et al., 2017). In this case, the product of the movement kernel ϕ(⋅) and the exponential selection function ω(⋅) is proportional to a single log-linear function of the corresponding model parameters and f_i reduces to:

(2)${\displaystyle f_i\left({\mathbf{x}}_{0,t+1}\mid {\mathbf{x}}_{0,t},{\mathbf{x}}_{0,t-1},\mathbf{Z};{\theta }_i,{\beta }_i\right)=\frac{exp\left({\mathbf{C}}_{0,t+1}^{\top }{\theta }_i+{\mathbf{Z}}_{0,t+1}^{\top }{\beta }_i-log\left(l_{0,t+1}\right)\right)}{{\int }_{\tilde{\mathbf{x}}\in D}exp\left({\tilde{\mathbf{C}}}^{\top }{\theta }_i+{\tilde{\mathbf{Z}}}^{\top }{\beta }_i-log\left(\tilde{l}\right)\right)d\tilde{\mathbf{x}}},}$ where $\tilde{\mathbf{C}}$, $\tilde{\mathbf{Z}}$ and $\tilde{l}$ depend on location $\tilde{\mathbf{x}}$. The vector C_0,t+1 can be interpreted as a movement covariate vector that contains different step length and turning angle terms. Its exact form depends on the chosen step length and turning angle distributions (Table S1, see also Avgar et al., 2016, Nicosia et al., 2017). For example, for gamma-distributed step length and von-Mises-distributed turning angles with either mean zero or π, we have C_0,t+1 = (log(l_0,t+1),  − l_0,t+1, cos(α_0,t+1))^⊤. The corresponding state-dependent movement coefficient vector is θ_i = (k_i − 1, r_i, κ_i)^⊤ with k_i and r_i being the shape and rate parameter of the gamma-distribution belonging to state i, respectively, and |κ_i| being the state-dependent concentration parameter of the von-Mises distribution(κ_i < 0 indicates a von Mises distribution with mean π; Nicosia et al., 2017). Thus, in this reduced representation of the step-selection density f_i, the parameterisation of the movement kernel might differ from the commonly used parameterisation of the corresponding step and angle distributions (e.g.,  k_i − 1 instead of k_i), but there is a direct relationship between the two (Tables S1 and S2). The negative log step length included in the exponential function is necessary to correctly represent the movement kernel in a Cartesian coordinate system.

The reduced form of f_i is very convenient. Justified by the law of large numbers, it allows for a joint parameter estimation of the state, movement and selection parameters based on a Markov-switching conditional logistic regression for case-control designs with M control, i.e., available, locations per observed, i.e., used, location (Nicosia et al., 2017, see also Avgar et al. 2016 for step-selection models without state-switching). This forms the basis for the standard HMM-iSSA.

Parameter estimation

The HMM-iSSA workflow is similar to the one of the standard iSSA. For each observed step, we choose M control steps, usually using a suitable parametric proposal distribution for step length and turning angle, respectively, and extract the corresponding habitat and movement covariate values. This builds the case-control data set. The model parameters are then estimated using a Markov-switching conditional logistic regression, i.e., a conditional logistic regression in which the regression coefficients depend on an underlying latent Markov chain. As Klappstein, Thomas & Michelot (2023), we use the forward algorithm, which is well-known especially in the context of HMMs (Zucchini, MacDonald & Langrock, 2016) to efficiently evaluate the corresponding likelihood. This allows for a numerical maximum likelihood estimation based on standard optimization procedures such as nlm in R (R Core Team, 2022). Afterwards, it is possible to decode the states, for example, using the Viterbi algorithm (Viterbi, 1967), which calculates the most likely sequence of states given the fitted model and the case-control data.

More precisely, for each step from location x_0,t to x_0,t+1 (t = 2, …, T − 1), we create a choice set ${\tilde{\mathbf{x}}}_{t+1}=\left\{{\mathbf{x}}_{0,t+1},{\mathbf{x}}_{1,t+1},\dots ,{\mathbf{x}}_{M,t+1}\right\}$ that includes the observed and the M control locations for the end point of the step. Usually, the control steps are randomly drawn from a suitable proposal distribution for step length and turning angle (Forester, Im & Rathouz, 2009). However, it is also possible to use a grid or a mesh (Arce Guillen et al., 2023). Here the devil is in the detail, as depending on the sampling procedure, the interpretation of the models’ movement coefficients might differ (see Section S2). The interpretation of the selection coefficients, however, remain unaffected if M is chosen large enough to provide stable results.

In this case-control setting, we model the state-dependent choice probability p_0,t,i of choosing the step to the observed location x_0,t+1 from the choice set ${\tilde{\mathbf{x}}}_{t+1}$ given the current state S_t = i, as:

(3)${\displaystyle p_{0,t,i}\left({\mathbf{x}}_{0,t+1}|{\tilde{\mathbf{x}}}_{t+1},\mathbf{C},\mathbf{Z};{\theta }_i,{\beta }_i\right)=\frac{exp\left({\mathbf{C}}_{0,t+1}^{\top }{\theta }_i+{\mathbf{Z}}_{0,t+1}^{\top }{\beta }_i\right)}{\sum _{m=0}^{M}exp\left({\mathbf{C}}_{m,t+1}^{\top }{\theta }_i+{\mathbf{Z}}_{m,t+1}^{\top }{\beta }_i\right)},}$ with C_m,t+1 and Z_m,t+1 being the movement and habitat covariate vectors belonging to location x_m,t+1 for m = 0, …, M. This case-control step-selection probability is closely related to direct numerical integration, which offers an alternative way to approximate the step-selection density f_i. Note that in contrast to Eq. (2), the term −log(l_⋅,t+1) does not appear in Eq. (3) because we change from Cartesian to polar coordinates when sampling control steps (i.e., step lengths and turning angles).

We derive the likelihood of the Markov-switching conditional logistic regression by plugging p_0,t,i into the HMM likelihood (Zucchini, MacDonald & Langrock, 2016), (4)${\displaystyle \mathcal{L}\left(\theta ,\beta ;{\tilde{\mathbf{x}}}_3,{\tilde{\mathbf{x}}}_4,\dots ,{\tilde{\mathbf{x}}}_T,\mathbf{C},\mathbf{Z}\right)={\delta }^{\top }\mathbf{P}\left({\tilde{\mathbf{x}}}_3\right)\mathbf{\Gamma }\mathbf{P}\left({\tilde{\mathbf{x}}}_4\right)\mathbf{\Gamma }\cdots \mathbf{\Gamma }\mathbf{P}\left({\tilde{\mathbf{x}}}_T\right)\mathbf{1},}$ where $\mathbf{P}\left({\tilde{\mathbf{x}}}_{t+1}\right)=\text{diag}\left(p_{0,t,1},\dots ,p_{0,t,N}\right)$ is a diagonal matrix including the state-dependent step-selection probabilities, Γ and δ are the transition probability matrix and the initial distribution of the underlying Markov chain, respectively, and 1 is an N-dimensional vector of ones. We can then estimate the model parameters using a numerical maximization of the log-likelihood (for details, see Zucchini, MacDonald & Langrock, 2016). In our implementation, we restrict the movement parameters to always remain in their natural parameter space, e.g., the shape and rate parameters of the gamma distribution are constraint to values greater than zero.

For initialization, the numerical maximization requires a set of starting values for the model parameters. To avoid ending up in a local maximum of the log-likelihood, it is necessary to test several sets of starting values, for example by randomly drawing initial values for each model parameter. We discuss this in more detail in Section S3.

Two-step approach

The TS-iSSA is based on the same idea as the HMM-iSSA. However, the TS-iSSA relies on a prior classification of the movement data into different movement states. Thus, in a first step, an N-state HMM with state-dependent step length and turning angle distributions as defined for the movement kernel is fitted to the data, e.g., using a gamma distribution for step length and a von Mises distribution for turning angles. Then, the Viterbi algorithm is used to assign each observed step to one of the N HMM movement states. Alternatively, local state decoding can be used. In the second step, state-specific (i)SSAs are applied to the N state-specific data sets using a case-control design and conditional logistic regression (e.g., Roever et al., 2014; Karelus et al., 2019). The control steps for the state-specific case-control data sets are usually sampled based on the respective state-dependent HMM step length and turning angle distributions.

Simulation study

We used a simulation study with three state-switching scenarios to evaluate the performance of our HMM-iSSA approach and to demonstrate possible consequences of either ignoring the underlying latent states in the traditional iSSA or ignoring the uncertainty of prior state-decoding in the TS-iSSA. For each scenario, we generated movement data from a hidden Markov step-selection model with 2 states and state transition probabilities γ₁₁ = γ₂₂ = 0.9. A realization of a Gaussian random field with covariance σ² = 1 and range ϕ = 10, computed using the function grf from the R-package geoR (Ribeiro Jr et al., 2022), served as the habitat covariate Z (Fig. S2). For the movement kernel, we used state-dependent gamma and zero-mean von Mises distributions to model step length and turning angle, respectively (Fig. 2).

Table 1 summarizes the movement and selection parameters for each of the three simulation scenarios. Scenario 1 is chosen to represent a typical inactive-active scenario in which the first state (“inactive” state) is associated to small step length, less directive movement and no selection, while the second state (“active” state) corresponds to larger step length, more directed movement and attraction to the landscape feature Z. The second and the third scenarios cover the rather extreme cases in which either the selection or the movement parameters are shared across states: In Scenario 2 (“switching preferences”), the two states only differ in their selection patterns with avoidance of the feature in state 1, and attraction to the feature in state 2. In Scenario 3 (“HMM”), only the movement patterns differ across states while there is no selection for or against the landscape feature in either state. This corresponds to a basic movement HMM. To check the robustness of the HMM-iSSA, we additionally covered a fourth scenario without state-switching in which the data are generated based on a standard step-selection model (Section S5.2). Furthermore, to check the influence of the spatial variation in the habitat feature on the estimation results, we also considered a second landscape feature map with weaker spatial heterogeneity (Section S5.3).

In each of the 100 simulation runs per scenario, we simulated movement paths of length T = 1000 from the corresponding hidden Markov step-selection model and then applied 2-state HMM-iSSAs, 2-state TS-iSSAs and iSSAs to corresponding case-control data sets with randomly drawn control steps using a uniform distribution for turning angle and a proposal gamma distribution for step length, respectively (Section S2). To check whether the parameter estimates converge to stable values, we considered M = 20, M = 100 and M = 500 control locations per observed step. For model selection purposes, we further estimated the parameters of a 2-state HMM-iSSA with movement but without habitat covariates. This corresponds to a basic movement HMM, but fitted to the same case-control data as the candidate models. After parameter estimation, we computed AIC and BIC values (Burnham & Anderson, 2002) and standard p-values for the estimated selection coefficients. For TS-iSSAs and HMM-iSSAs we further computed the state missclassification rate, i.e., the percentage of states that were not correctly classified using the Viterbi algorithm. These metrics and the parameter estimates were used to evaluate the estimation and classification accuracy of the candidate models, and to assess the performance of standard model selection procedures.

Case study on bank vole interactions

To illustrate the use of HMM-iSSAs on empirical data, we applied them to existing movement data of synchronously tracked bank vole individuals (Myodes glareolus) as analyzed in Schlägel et al. (2019). The original sample included 10 group-level replicates with 4 animals each, combining animal personalities (not part of this study) and sex. All bank vole females were sexually mature, but not gravid or lactating. All bank vole males were mature. Individuals within a replicate were synchronously tracked in fenced quadratic outdoor enclosures of 2,500 m² for 3–5 days using collars with small radio telemetry transmitters (1.1 g, BD-2C, Holohil Systems Ltd., Canada) and a system of automatic receiving units (Sparrow systems, USA). For bank vole individuals tracked under natural conditions, the estimated home range sizes were on average 2, 029m² with a core area of 549m² (Schirmer et al., 2019). Thus, the size of the enclosures allowed the individuals to express their natural movement and space use. The original sample size allowed for potential technical failures of tracking equipment or escapes of single individuals from semi-natural enclosures, which indeed reduced the sample size to 8 populations with 2 males and 1–2 females each and a total of 28 animals. The tracking system produced 6-minute location data. Due to daily system maintenance, locations were missing for approximately one hour per day. Otherwise, movement paths were complete. Depending on the replicate, this resulted in 602–1,200 locations per individual split into 3–5 bursts of around 23 h each. After the experiment animals were captured from the enclosures, transmitter collars were removed, and animals were returned to the wild at their individual capture location. For more details on the bank vole capture and location data, see Schlägel et al. (2019). Details on animal care approvals and research permissions are provided in the Section 6.1.

To study interactions between the bank vole individuals, i.e., attraction, avoidance or neutral behavior towards each other, Schlägel et al. (2019) applied SSAs to each individual of each replicate, respectively, using occurrence estimates of the conspecifics as covariates. The occurrence estimate of an individual provides a map of the individual’s space use during a certain time window, indicating areas of higher and lower probability of occurrence during that time period. It is estimated from the discrete sample of observed locations through kriging (Fleming et al., 2016). To account for the movement of individuals, occurrence estimates were computed using a rolling time window (here 4 h).

The analysis focused on interactions between males and females: Males were expected to mainly show attraction towards females, while females could show any of the three behaviors depending on their reproductive state (Schlägel et al., 2019). The authors suggested that the relatively large number of non-significant selection coefficients, especially found for male interactions with females, might be caused by unobserved mixtures of different underlying behavioral modes. Bank voles are polyphasic with resting phases of approximately 3h and active phases of approximately 1h following each other (Mironov, 1990). We therefore applied 2-state HMM-iSSAs to the same data to investigate (i) if the state-switching model is capable to detect meaningful biological states, and (ii) if we find different significant selection coefficients using the state-switching approach.

For each individual, we used a 2-state HMM-iSSA with state-dependent gamma distributions for step length, and uniform distribution for turning angle, respectively. Occurrence estimates of each conspecific within the same replicate were used as covariates for the selection part of the model (Schlägel et al., 2019). As interactions were analyzed within replicates, i.e., only for animals tracked simultaneously in same enclosure, a nested approach controlled for potential confounds such as differences among enclosures. We did not include a resource covariate, as vegetation was sufficiently homogeneous within enclosures. Furthermore, we chose to use M = 500 control steps per used step, as preliminary analysis with increasing values for M provided stable results for this choice (Fig. S9). Thus, with M = 500 control steps, the corresponding selection covariate vector for individual k at time t and locations x_m,t, m = 1, …, 500, was given by Z_k,m,t = ({O_−k,m,t}), where {O_−k,m,t} denotes the set of occurrence estimates of the respective conspecifics withing the same replicate. The corresponding movement covariate vector was C_k,m,t = (log(l_k,m,t),  − l_k,m,t). Parameters were then estimated using 50 sets of random starting values. For model comparison, we further applied corresponding iSSAs (no state-switching) and HMM-iSSAs without selection covariates (i.e., HMM approximations; no selection) to the same case-control data set for each individual. For the 2-state TS-iSSAs (prior state-classification), control steps were sampled using the estimated state-dependent distributions of the movement HMMs fitted in the first step of the TS-iSSA.

Both the simulation and the case study were implemented in R (R Core Team, 2022). We used custom code for HMMs and HMM-iSSAs and the clogit-function from the R package survival (Therneau, 2020) for the iSSA and corresponding parts of the TS-iSSA. The R code is available in Supplementary Information 2.

Simulation study

Overall, the HMM-iSSA performed very well across all simulation scenarios and did not produce any evident bias even in the extremer Scenarios 2 (“state-switching preferences”) and 3 (“HMM”; Fig. 3, Tables S3 and S4). The TS-iSSA was able to detect two suitable states in both scenarios with state-dependent movement kernels (Scenarios 1 and 3), although there was a small but evident bias for some parameters, for example, for the selection coefficients in scenario 1 (0.18 in state 1, −0.13 in state 2 for M = 500), and for the shape parameter in scenario 3 (0.12 in state 2 for M = 500). For Scenario 1 (“active-inactive”), this is also reflected in the rather large percentage of significant selection coefficients across the simulation runs in state 1 (39–42% at a significance level of α = 0.05, Table 2), although the true coefficient is equal to zero. Thus, in contrast to the HMM-iSSA, the p-values of the TS-iSSA are not reliable in this active-inactive setting.

The iSSA is by its nature unable to distinguish between the underlying states and thus, did not recover the true underlying parameters in either scenario. Especially in Scenario 2 (“switching preferences”), the iSSA selection coefficients were estimated close to zero and the associated p-values would misleadingly indicate no selection for or against the landscape feature in 42% − 43% of the simulation runs (Table 2). Note that the TS-iSSA produced similar results to the iSSA in this scenario, since the inherent HMM classification was not able to distinguish between states that share the same movement kernel, and therefore all steps were classified to belong to the same state.

For all three simulation scenarios, the number of available steps M only slightly affected the estimation results in this simulation exercise, especially the results for M = 100 and M = 500 are very similar (Fig. 3). Thus, the results seem to be stable. The HMM-iSSA with M = 500 available steps achieved the lowest missclassification rate (Table 3). As the data in Scenario 3 were simulated from an HMM, the HMM classification was equally accurate in this scenario. Overall, the HMM-iSSA clearly outperformed the other candidate models in its estimation and classification performance in all scenarios.

As the TS-iSSA involves an a-priori HMM classification, it does not provide a proper maximum likelihood value. It is therefore not possible to calculate corresponding AIC or BIC values for model selection. Thus, we only considered iSSAs without state-switching, HMMs without selection (fitted to the same case-control data sets) and HMM-iSSAs as candidate models to evaluate information-criteria based model selection in this modeling framework. For Scenario 1 and 2, AIC and BIC performed very well and selected the true underlying model in 100 and 98% of the simulation runs, respectively (Table 4). In Scenario 3 (“HMM”), the AIC tended to select the true HMM model in most of the cases but occasionally selected the more complex HMM-iSSA (9 − 14% of the runs), while the BIC again selected the correct model in all simulation runs.

Overall, the simulation runs with lower spatial variation in the landscape variable produced similar results (Section S5.3). However, the lower spatial variation reduces the influence of the habitat selection function on space use. Therefore, the variance in the estimates slightly increased, the HMM missclassification rate decreased in Scenario 1 and the HMM-iSSA missclassification rate increased in Scenario 2. In the supplementary simulation scenario without state-switching, the HMM-iSSA was able to recovered the true underlying values in state 1, but produced unusable estimates in state 2 (Section S5.2).

Case study on bank vole interactions

For most bank vole individuals, the HMM-iSSA approach could reasonably distinguish between two activity levels. State 1 was always associated to shorter step lengths compared to state 2 which could correspond to a rather inactive behavior (Fig. 4; mean of the estimated gamma distribution for step length ranging from 1.43 to 7.38 in state 1, and from 8.13 to 20.85 in state 2, respectively). According to the Viterbi decoded state sequences, the “less active” state 1 was occupied between 15.29% and 66.71% of the observed time period (Table S7), except for male 1 in replicate 4 which spent 96.43% of the time in state 1 according to its decoded state sequence. It is also the individual with the largest estimated mean step length in both states (7.38 in state 1, 20.85 in state 2). Thus, for this male, interpretation must be taken with care. The TS-iSSA provided mostly similar results for the movement kernel and state classification (Fig. S11).

For 21 of the 28 bank vole individuals, the HMM-iSSA results implied neutral behavior towards conspecifics in state 1 as all selection coefficients were non-significant (α = 0.05; Fig. 5 and Fig. S10). This matches well with the interpretation of a less active/inactive state. For two individuals, the results indicated avoidance behavior in state 1. In state 2 (“active state”), most bank voles showed attraction to at least one bank vole of opposite sex as implied by the positive and significant selection coefficients. However, for four males and four females, the coefficients for occurrence of individuals with opposite sex were non-significant in both states. These are mainly the individuals for which the iSSA also implied neutral behavior (Fig. 5). However, for 3 individuals, i.e., male 1 in replicate 7, female 1 in replicate 4, and female 1 in replicate 8, the HMM-iSSA indicated attraction towards another individual of opposite sex, while the iSSA indicated neutrality. The opposite is true for female 1 in replicate 7 for which only the iSSA indicated attraction. The selection coefficients for occurrence of individuals with same sex usually implied neutral behavior in state 1, and neutral or attraction behavior in state 2 (Fig. S10).

Overall, the results of the TS-iSSA are in line with the results of the HMM-iSSA (Figs. S10, S12), although the implications are slightly different for nine individuals. Regarding information-criteria based model selection, for most bank voles, AIC and BIC pointed to the hidden Markov step-selection model (Table S8). However, for 10 individuals, including half of the female individuals, BIC selected a simpler model, i.e., iSSA or HMM. The selection of HMMs mainly corresponded to cases with many non-significant HMM-iSSA selection coefficients. The iSSA was preferred by BIC for male 1 in replicate 4 and female 2 in replicate 8.