An introduction to statistical models used to characterize species-habitat associations with animal movement data

Abstract

Understanding species-habitat associations is fundamental to ecological sciences and for species conservation. Consequently, various statistical approaches have been designed to infer species-habitat associations. Due to their conceptual and mathematical differences, these methods can yield contrasting results. In this paper, we describe and compare commonly used statistical models that relate animal movement data to environmental data. Specifically, we examined selection functions which include resource selection function (RSF) and step-selection function (SSF), as well as hidden Markov models (HMMs) and related methods such as state-space models. We demonstrate differences in assumptions while highlighting advantages and limitations of each method. Additionally, we provide guidance on selecting the most appropriate statistical method based on the scale of the data and intended inference. To illustrate the varying ecological insights derived from each statistical model, we apply them to the movement track of a single ringed seal (Pusa hispida) in a case study. Through our case study, we demonstrate that each model yields varying ecological insights. For example, while the selection coefficient values from RSFs appear to show a stronger positive relationship with prey diversity than those of the SSFs, when we accounted for the autocorrelation in the data none of these relationships with prey diversity were statistically significant. Furthermore, the HMM reveals variable associations with prey diversity across different behaviors, for example, a positive relationship between prey diversity and a slow-movement behaviour. Notably, the three models identified different “important” areas. This case study highlights the critical significance of selecting the appropriate model as an essential step in the process of identifying species-habitat relationships and specific areas of importance. Our comprehensive review provides the foundational information required for making informed decisions when choosing the most suitable statistical methods to address specific questions, such as identifying expansive corridors or protected zones, understanding movement patterns, or studying behaviours. In addition, this study informs researchers with the necessary tools to apply these methods effectively.

Background

Understanding the relationships between animals’ space use and their physical environment is a fundamental aspect of ecological research and is essential for species conservation [1–3]. Movement data from biologging devices provide valuable information on animal space use [4, 5], and statistical models are often used to relate the movement data (in terms of occurrence, movement, or behaviour) to indicators of resources (e.g., vegetation [6]), proxies of energy (e.g., terrain ruggedness [7]), or perceived predation risk (e.g., canopy cover, [8]). These models can provide insight into the environmental conditions associated with key ecological concepts such as home range, habitat selection, movement corridors, behaviour, and critical habitat [5, 9]. Identifying and designating critical habitat is one of the main ways governments protect species from rapid and widespread habitat degradation and climate changes (e.g., under the Environment Protection and Biodiversity Conservation (EPBC) Act in Australia, Species at Risk Act (SARA) in Canada, Wildlife and Countryside Act in the United Kingdom, Endangered Species Act (ESA) in the United States). Thus, it is essential that statistical models that use movement data to identify critical habitat are chosen, implemented, and interpreted properly.

Various models aim to link animal movement data to environmental covariates, but each model is appropriate for specific research questions [10, 11]. Many of these models can be easily implemented using R packages (e.g., amt, [12]; momentuHMM, [13]; see [14] for a relevant review of packages), and their use is ubiquitous in animal movement literature (e.g., [15, 16]). However, choosing an appropriate method is a complex and challenging aspect of movement modeling [14]. For example, different models address ecological questions at different scales, from large-scale questions on important areas for species (e.g., first- or second-order selection at the species or home range scale, respectively), to smaller-scale movement- or behaviour-specific questions (e.g., third- or fourth-order selection at the habitat usage or actual food intake scale, respectively; [17–19]). Some methods are designed to understand habitat preference via selection functions (e.g., resource selection function [RSF], step-selection function [SSF], or integrated-SSF) and others are focused on understanding how habitat relates to discrete behavioural states (e.g., hidden Markov model [HMM]). As such, consideration of scale and the behaviour of interest when choosing a model is imperative for meaningful interpretation of results.

Here, we describe and compare three mainstream, and readily implemented, models for linking animals’ movement data to environmental covariates: RSF, SSF, and HMM. While these models answer different ecological questions and require different resolutions of data, they are all commonly applied to characterize relationships between species and their environments using movement data. Our objective is to help clarify the intended use and mathematical underpinnings of each model to help ecologists properly choose between them. For example, RSFs and SSFs are typically used to address similar questions on habitat selection, yet SSFs generally require relatively high-frequency data compared to RSFs. SSFs and HMMs generally both require movement data at a finer temporal resolution, yet SSFs are commonly used to address questions regarding movement and habitat selection whereas HMMs are used to link specific animal behaviours to environmental covariates. Reviews comparing selection functions (e.g., [20, 21]) or selection functions with inhomogeneous point process models (IPPs, e.g., [22]) are available, but here, our review is unique as we compare selection functions and HMMs. HMMs are a fundamentally different, but increasingly popular, model to relate animal movement data to environmental covariates, and thus, a review of RSFs, SSFs, and HMMs is warranted.

Description of statistical models

Resource selection functions (RSF)

A resource selection function (RSF) is a widely used function that relates habitat characteristics to the relative probability of use by an animal [23–26]. RSFs can assume that the environmental conditions selected by an animal provide desirable resources, and that the probability of habitat selection (Table 1) is a function of predictor variables representing a resource distribution on a landscape [27]. Although RSFs are now often referred to as habitat selection functions as they can include covariates other than resources (e.g., predator probability; [28]).

Table 1 Terminology used in this paper

When applied to movement data, RSFs often compare observed animal locations, which represent a sample of the habitat used by the animal, to a sample of locations randomly selected within an animal’s home range (e.g., the minimum convex polygon (MCP), or extent, of the observed locations, Fig. 1). Observed locations are often referred to as used locations, while the random background locations are often referred to as available locations [26]. RSFs are widely applied in part due to their ease of use (e.g., using the amt R package, [12]) and that they can provide broad scale information on species-habitat relationships [21].

Fig. 1

A observed (i.e., “used”) animal locations, B the home range (here via minimum convex polygon), and C available locations randomly sampled from within the minimum convex polygon

The RSF, $w(\mathbf{x})$, encodes the relative strength with which an animal selects a given set of habitat covariate values $\mathbf{x}$ over all habitat covariate values available to them (e.g., one of these covariates x may be the density of food at a location [25, 29, 33, 34]). The RSF, $w(\mathbf{x})$, is typically defined and interpreted in the exponential form:

$$w\left( {\mathbf{x}} \right) = {\text{exp}}\left( { \beta_{1} x_{1} + \beta_{2} x_{2} + \cdot \cdot \cdot + \beta_{k} x_{k} } \right) ,$$

(1)

where $\mathbf{x}=\{{x}_{1},\dots , {x}_{k}\}$ denotes the values of k predictor habitat variables and ${\beta }_{1}$,…, ${\beta }_{k}$ are the associated selection coefficients [23]. In practice, the coefficients of the RSF (${\beta }_{1}$,…, ${\beta }_{k}$) are commonly estimated using a logistic regression. Namely, suppose that a total of $n$ resource units are monitored, and the response $\mathbf{y} =\{{y}_{1},...,{y}_{n}\}$ is a set of binary random variables representing observed (${y}_{i}=1)$ and available (${y}_{i}=0$) resource units. Then, the probability that resource unit $i$ is used given its environmental covariates ${\mathbf{x}}_{{\varvec{i}}}=\{{x}_{1,i},\dots , {x}_{k,i}\}$ is modeled as:

$$Pr(y_{i} = 1|{\mathbf{x}}_{i} ) = \frac{{{\text{exp}}\left( {\beta_{1} x_{1,i} + \beta_{2} x_{2,i} + \cdot \cdot \cdot + \beta_{k} x_{k,i} } \right)}}{{1 + {\text{exp}}\left( {\beta_{1} x_{1,i} + \beta_{2} x_{2,i} + \cdot \cdot \cdot + \beta_{k} x_{k,i} } \right)}} .$$

(2)

The RSF from Eqs. 1–2 is written in environmental space, but RSFs can also be written in geographic space as an inhomogeneous Poisson point process (IPP). RSFs written as IPPs describe the density of points (e.g., animal locations) across all physical locations available to the animal (e.g., s may be the easting and northing coordinate of an animal in Eq. 3) [29, 35, 36]. An important difference between the two is that RSFs written in environmental space generally model directly the relationships with environmental covariates, $w\left(\mathbf{x}\right),$ while IPPs account for the relationship with environmental covariates through modelling the density of points [26, 37]. Specifically, the IPP models the density of points as the exponential of a linear function of spatial predictors—interpreted as the intensity function, λ(s):

$${\uplambda }\left( s \right){ } = {\text{exp}}\left( {\beta_{0} + \beta_{1} x_{1} \left( s \right) + { }\beta_{2} x_{2} \left( s \right){ } + { }.{ }.{ }.{ } + { }\beta_{k} x_{k} \left( s \right)} \right) ,$$

(3)

where s is the location in geographical space, ${x}_{1}$(s), …, ${x}_{k}$(s) are k predictor habitat variables associated with location s, ${\beta }_{0}$ is an intercept term, and${\beta }_{1}$, …, ${\beta }_{k}$, are the selection coefficients similar to Eq. 1. The intensity function $\lambda (s)$ is proportional to $w(\mathbf{x})$ from Eq. 1 and the selection coefficients can be estimated with a logistic regression. As the number of available points sampled for the logistic regression grows to infinity, the selection coefficients (${\beta }_{1}$, …,${\beta }_{k}$) from Eq. 2 (which are estimated using logistic regression) converge to the selection coefficients (${\beta }_{1}$, …,${\beta }_{k}$) from Eq. 3 and the IPP model [37]. Note that if the model intercept ${\beta }_{0}$ is estimated it must be interpreted with care as it is often not biologically meaningful (i.e., it is related to the density of observations, but for movement data that may reflect frequency of observations; [26, 34]).

While the RSF $w(\mathbf{x})$ represents selection, it more formally represents the ratio between the use (i.e., utilization) distribution ${F}^{U}(\mathbf{x})$ (the frequency distribution of habitat covariates used by an animal) and available distribution ${F}^{A}(\mathbf{x})$ (the frequency distribution of habitat covariates assumed to be available) (i.e., $w(\mathbf{x})\propto {F}^{U}(\mathbf{x})/{F}^{A}(\mathbf{x})$). This relationship can be seen from the following definition of the use distribution, ${F}^{U}(\mathbf{x})$:

$$F^{U} \left( {\mathbf{x}} \right) = \frac{{w\left( {\mathbf{x}} \right)F^{A} \left( {\mathbf{x}} \right)}}{{\mathop \smallint \nolimits_{{{\mathbf{x^{\prime}}}\epsilon {\Omega }}}^{{}} w\left( {{\mathbf{x^{\prime}}}} \right)F^{A} \left( {{\mathbf{x^{\prime}}}} \right)d{\mathbf{x^{\prime}}}}} ,$$

(4)

where in the denominator we integrate over the entire domain of the environment, $\Omega$, to ensure that ${F}^{U}(\mathbf{x})$ is a valid probability distribution that integrates to one. Equation 4 shows that the use distribution is a function of both the available habitat (as defined by ${F}^{A}(\mathbf{x})$) and the selection function, $w(\mathbf{x})$ [25, 29, 38, 39]. As movement datasets only contain information on used habitat variables, available resource units, which should not erroneously be treated as true absences, must be sampled from ${F}^{A}$. A similar definition of the use distribution can be made in geographical space:

$$F^{U} \left( s \right) = \frac{{w\left( {{\mathbf{x}}\left( s \right)} \right) F^{A} \left( s \right)}}{{\mathop \smallint \nolimits_{s^{\prime}\epsilon G}^{{}} w\left( {{\mathbf{x}}\left( {s^{\prime}} \right)} \right)F^{A} \left( {s^{\prime}} \right)ds^{\prime}}} ,$$

(5)

where the selection function $w(\mathbf{x}(s))$ is now a function of $\mathbf{x}(s)$, the habitat covariates x at location $s$, and the available distribution ${F}^{A}(s)$ is a function of location and is generally assumed to be constant (i.e., ${F}^{A}(s) = \frac{1}{\left|G\right|}$; [26]). In the denominator, we integrate over the geographical area available to the animal $G$ [38].

RSFs’ assumptions of independence and large-scale availability make them particularly well-suited for analyzing coarse animal movement data to understand the broad-scale spatial ecology of animals. For example, RSFs have been used to identify broad conservation corridors [39] and high wildlife density areas to improve rabies vaccination programs [35]. Additionally, RSFs have been used to understand how habitat selection is affected by seasons [36] anthropogenic developments [40], and the presence of predators [28], and can be used to understand the cumulative effects of such factors [41].

When applied to movement data, simple RSFs are most suitable for data collected at a coarse temporal resolution. Autocorrelation in movement data can introduce spatio-temporal correlation that cannot be explained by habitat-related covariates. If not accounted for, this autocorrelation can result in underestimation of standard errors, potentially inflating Type I error rates [29, 42]. This issue is particularly common when using high-resolution datasets [43] and is expected to become more significant with the advancement of technology that enables increasingly fine-resolution movement data [44]. Various approaches have been proposed to mitigate this problem, including data thinning [45], weighting points [46], correcting the standard errors [47], and adding an autocorrelation term (e.g., intrinsic conditional autoregressive (ICAR) model, [48]).

Step selection functions (SSF)

A step selection function (SSF) builds on RSFs and estimates an animals’ resource selection at each observed sequential location or step (i.e., the linear segment between two consecutive locations; Fig. 2) based on both habitat covariates and movement constraints [49–51], partially providing a solution to the autocorrelation issues associated with RSFs (but see [45, 52]). Unlike RSFs, which assume the entire home range is available over the study period, SSFs assume that an animal can only access a limited area near its previous location within the time scale of a step. There are many versions of SSFs, including many versions that focus simply on the influence of environmental covariates on selection. In our review, we focus on a commonly applied version often referred to as integrated-SSFs [53] because it allows for environmental covariates to influence animal movement itself in addition to selection. These SSFs allow the movement characteristics (including step length and turning angle) and environmental covariates to be combined in a single linear predictor, which allows for joint estimation, and potential interactions, between movement and environment, thereby relaxing the assumption that observed movement attributes (i.e., velocities and directional persistence) are independent of resource selection [53, 54]. That is, these SSFs simultaneously estimate the habitat selection and movement kernel. The movement kernel describes how an animal would move in the absence of habitat selection and thus what is available to the animal. In these SSFs, where environmental covariates influence both animal movement and selection, availability is estimated simultaneously with habitat selection through the movement kernel, rather than being assumed a priori [50, 53, 55].

Fig. 2

Example of an animal movement track with linear observed steps between locations. An observed turning angle is relative to the trajectory from the preceding step. At each observed location, random control locations (m = 3; and their associated step lengths and turning angles) are generated to estimate availability

Step selection functions model habitat selection and movement explicitly, by defining use distribution as the product of the habitat selection w and a selection-free movement kernel $\phi$ [53, 54]. The movement kernel is essentially equivalent to having the availability distribution of an RSF, ${F}^{A}(\mathbf{x})$, vary temporally, and it represents how the animal would move in the absence of habitat selection [53, 54]. Thus, the probability density of an animal using location ${s}_{t}$ at time $t$ given its previous locations ${\mathbf{s}}_{1:t-1}=\{{s}_{1},{s}_{2},\dots , {s}_{t-1}\}$ is modeled as:

$$F^{U} {(}s_{t} {|}{\mathbf{s}}_{1:t - 1} ) = \frac{{w\left( {s_{t - 1} ,s_{t} } \right)\phi \left( {s_{t} , {\mathbf{s}}_{1:t - 1} } \right) }}{{\mathop \smallint \nolimits_{{s_{{}}^{,} \epsilon G}}^{{}} w\left( {s_{t - 1} ,s_{{}}^{,} } \right)\phi \left( {s_{{}}^{,} , {\mathbf{s}}_{1:t - 1} } \right)ds_{{}}^{,} }} ,$$

(6)

where $w({s}_{t-1},{s}_{t})$ is the habitat selection function and $\phi ({s}_{t} , {\mathbf{s}}_{1:t-1})$ is the selection-free movement kernel (although this terminology is inconsistently used, see [55]). As in Eq. 5, the denominator integrates over the entire geographic domain $G$ to normalize the use distribution. Step selection functions typically require observed locations to be at a constant time interval, and the parameters of the SSF and movement kernel are dependent on time index $t$ [53, 54]. Similar to RSF coefficient estimation with a logistic regression, SSFs coefficients are often estimated by fitting a conditional logistic regression. SSFs coefficients can also be estimated using Monte Carlo integration techniques [55].

Fitting a conditional logistic regression to estimate the SSF coefficients generally entails drawing $m$ control points at each time t. We denote the vector with the locations of observed and control points as ${\mathbf{s}}_{t,0:m}=\{{s}_{t,0},{s}_{t,1},\dots , {s}_{t,m}\}$, where ${s}_{t,0}$ is the location (i.e., the coordinates) of the observed point at time t and ${s}_{t,i}$ for $i= 1, \dots , m$ is the location of the i^th control point. Further, ${y}_{t,i}$ is a binary variable that equals 1 when $i=0$ (observed) and 0 when $i= 1,\dots ,m$ (control). Control points (also known as random, sample, or integration points) are often generated using observed step length and turning angle distributions (Fig. 2; [55]). While such control points can be misinterpreted as solely representing available movement options, they are often used as a computationally efficient tool to estimate availability at high resolution close to the animal’s last location [55]. The conditional probability that the point at location ${s}_{t,i}$ is the observed location given the environmental covariates at that location ${\mathbf{x}}_{t,i}$ and the animal’s previous locations ${\mathbf{s}}_{1:t-\text{1,0}}$ [53, 55] is:

$$Pr\left( {y_{t,i} = 1{|}{\mathbf{X}}_{t} ,{\mathbf{s}}_{t,0:m} ,{\mathbf{s}}_{1:t - 1,0} } \right) = \frac{{\exp \left( {\eta_{t,i} } \right)}}{{\mathop \sum \nolimits_{j = 0}^{m} \exp \left( {\eta_{t,j} } \right) }} ,$$

(7)

where ${\mathbf{x}}_{t,i}$ is row $i$ of the matrix ${\mathbf{X}}_{t}$ and ${\eta }_{t,i}$ is a linear predictor that arises from the fact that both the habitat selection function and the movement kernel of Eq. 6 can take an exponential form. The linear predictor typically includes a vector of k environmental covariate values ${\mathbf{x}}_{t,i}$ at location ${s}_{t,i}$ and movement variables defined by previous and current locations, such as the natural log of step length and cosine of turning angle as covariates (e.g., to model step length with a gamma distribution and turning angle with a von Mises distribution [43, 55]). To model a movement kernel that can change as a function of environmental covariates (e.g., slower walking speed in deep snow), ${\eta }_{t,i}$ can include interactions between movement variables and the environmental covariates, including, for example, environmental covariates at the previous observed location ${s}_{t-\text{1,0}}$. For example, the linear predictor can be expressed as:

$$\eta_{t,i} = \beta_{1} x_{1,t,i} + \beta_{l} l_{t,i} + \beta_{{{\text{ln}}}} {\text{ln}}\left( {l_{t,i} } \right) + \beta_{\theta } \cos \left( {\theta_{t,i} } \right) + \beta_{2} x_{2,t - 1,0} \ln \left( {l_{t,i} } \right) ,$$

(8)

where ${x}_{1,t,i}$ represents an environmental covariate at location ${s}_{t,i}$ that affects selection (e.g., food density), ${x}_{2,t-\text{1,0}}$ an environmental covariate at location ${s}_{t-\text{1,0}}$ that affects movement itself (e.g., snow depth), ${l}_{t,i}$ the step length between observed location ${s}_{t-\text{1,0}}$ and potential location ${s}_{t,i}$, ${\theta }_{t,i}$ the turning angle between observed locations ${\mathbf{s}}_{t-2:t-\text{1,0}}$ and potential location ${s}_{t,i}$, and the $\beta$s are coefficients to estimate. This formulation allows SSFs to integrate both movement dynamics and environmental influences into step-selection probability estimates.

Step selection functions are well-suited for fine-scale movement datasets and have been used for understanding how environmental covariates affect animal movement characteristics [56]. For example, SSFs have been used to show how roads increase step lengths [57]. SSFs have also been used to understand how predators or social networks influence movement [49, 58]. Further, SSFs have been used to determine how variation in temperature affects movement and use of thermal cover [15]. Additionally, SSFs have been used to understand how slopes associated with thermal uplift affect soaring bird movement patterns, including potential stopover movements to access carrion [59].

Hidden Markov models (HMM)

In the context of movement ecology, hidden Markov models (HMMs) are principally used to classify movement into a finite number of discrete behavioural states (e.g., searching, resting, traveling, [60]), and differ conceptually from RSFs or SSFs as they generally do not make inference on habitat selection. However, HMMs can be used to investigate relationships between environmental covariates and animal movement and behaviour [61, 62].

HMMs assume there is a hidden state process unfolding over time from which we obtain observation data (Fig. 3). That is, it is assumed that at any time, animals exhibit one of $N$ discrete and temporally autocorrelated states Z_t  ∈ {1, 2, …, N}, which represent hidden (unobserved) behaviours. These states are related to an observation time series (Y₁, …, Y_T), where T is the total number of time steps. In particular, given the hidden state process, Z₁, …, Z_T, each observation, Y_t, is assumed to depend only on its corresponding hidden state, Z_t (Fig. 3). In contrast to the binary (observed/available) values used in the observation time series of RSFs and SSFs, the observation used in animal movement HMMs, generally comprise time series of step lengths and turning angles [63, 64] and can include additional data streams such as dive depth or acceleration [65].

Fig. 3

The structure of a hidden Markov model, where the hidden states (Z_t) depend on the previous state, and the observed data (Y_t) depends on the hidden state. In the case of modeling animal movement, the observed process is composed of empirical data such as step lengths and turning angles, and the hidden state may represent behaviours such as foraging, resting, and traveling

For an N-state HMM (e.g., N = 3 states, representing foraging, resting, and traveling), the state transition probability matrix (${\varvec{\Gamma}}$) is an N x N matrix, where entry (i,j) is denoted as ${\gamma }_{i,j}$ and represents the probability of transitioning from state i to state j. Namely, we have ${\gamma }_{i,j}= Pr({Z}_{t+1}=j | {Z}_{t}= i)$, and each row of ${\varvec{\Gamma}}$ must sum to 1 (i.e.,${\sum }_{j=1}^{N}{\gamma }_{i,j}=1$). All observations are assumed to be independent of one another after conditioning on the underlying hidden state process: $Pr({y}_{t}|{Z}_{1}, ... , {Z}_{t}; {y}_{1}, ... , {y}_{t-1},{y}_{t+1}, ... , {y}_{T}) = Pr ({y}_{t}| {Z}_{t})$. We define an N x N diagonal matrix, $\mathbf{P}({y}_{t}) = diag({p}_{1}({y}_{t}), ... , {p}_{N}({y}_{t}))$, where ${p}_{i}({y}_{t})$ represents the probability density of observation ${y}_{t}$ given state $i$ (i.e., ${p}_{i}({y}_{t}) = Pr({Y}_{t}={y}_{t}|{Z}_{t}=i$)). When modeling animal movement, step length is often modeled with distributions such as gamma, exponential, or Weibull, while turning angle is modeled with circular distribution such as the von Mises or wrapped Cauchy [64].

Taken together, the HMM likelihood function (L) is:

$$L = {\mathbf{\delta P}}(y_{1} ){\mathbf{\Gamma P}}\left( {y_{2} } \right){\mathbf{\Gamma P}}\left( {y_{3} } \right), ...,\user2{ }{\mathbf{\Gamma P}}\left( {y_{T} } \right)1^{\prime} ,$$

(9)

where ${\varvec{\updelta}}$ is the initial state distribution of the Markov chain, and $1^{\prime}$ is a row vector of ones. The parameters can be estimated using numerical methods to maximize the likelihood (Eq. 9) or with the expectation–maximization (EM) algorithm [62].

Environmental covariates can be incorporated into two different parts of an HMM. First, to examine how covariates affect the probability to transition between states (e.g., how vegetation density may increase the probability to transition into foraging for an herbivore), they can be used to alter the transition probabilities, ${\gamma }_{i,j}$. Specifically, we can make transition probabilities, now denoted as ${\gamma }_{t,i,j}$, change at each time t as a function of a vector of environmental covariates ${\mathbf{x}}_{t}$, e.g., via a multinomial logit link:

$$\gamma_{t,i,j} = \frac{{{\text{exp}}\left( {\eta_{t,i,j} } \right)}}{{\mathop \sum \nolimits_{l = 1}^{N} {\text{exp}}\left( {\eta_{t,i,l} } \right)}} ,$$

(10)

and

$$\eta_{t,i,j} = \beta_{0,i,j} + \beta_{1,i,j} x_{1,t} + ... + \beta_{k,i,j} x_{k,t} ,$$

(11)

where ${\beta }_{0,i,j}$ is an intercept term and ${\beta }_{1,i,j},\dots ,{\beta }_{k,i,j}$ are the behaviour-specific regression coefficients for each $k$ environmental covariates with values ${{\mathbf{x}}_{t}=x}_{1,t},\dots ,{x}_{k,t}$ at time $t$ [63, 66, 67]. To ensure model identifiability, ${\beta }_{0,i,j},{\beta }_{1,i,j},\dots ,{\beta }_{k,i,j}$ are fixed to 0 for one element in each row of the state transition probability matrix, typically, the diagonal when $i=j$.

Second, environmental covariates (typically abiotic) can be used to alter the observation probability density function, $p_{i} ( y_{t} )$, which can be interpreted as examining how they affect the movement of the animal within state $i$ (e.g., how snow depth affects the speed of a walking animal). To do so, the observation probability density functions will change through time, now denoted $p_{t,i} (y_{t})$, and their parameters will depend on environmental covariates. For example, environmental covariates can affect step length mean $\mu_{t,i}^{\left( l \right)}$ for each state i at each time step t as

$$ln\left( {\mu_{t,i}^{\left( l \right)} } \right) = \beta_{0,i}^{{\left( {\mu^{\left( l \right)} } \right)}} + \beta_{1,i}^{{\left( {\mu^{\left( l \right)} } \right)}} x_{1,t} + ... + \beta_{k,i}^{{\left( {\mu^{\left( l \right)} } \right)}} x_{k,t} .$$

(12)

Similarly, environmental covariates can affect other observation parameters such as the standard deviation of the step length distribution and turning angle concentration (see Online Appendix 2.3). Additionally, movement can be modelled with an orientation bias (i.e., biased random walks) or further, with bias toward any angle relative to stimulus (e.g., olfactory search perpendicular to wind; [68]).

Hidden Markov models are specific types of state-space models [69], where the states are discrete (i.e., there is a finite number of hidden states). Other state-space models have been developed that describe animal movement behaviour as a continuum of move-persistence ranging from 0 to 1. In these models, the state represents the move persistence where values closer to 0 are indicative of less move persistence (i.e., area-restricted search), and values closer to 1 are indicative of more persistent movement (i.e., traveling; [70]). However, these state-space models assume that behaviour exists on a linear continuum of only two behaviours (e.g., area-restricted search to traveling; e.g., [71]), and it is unclear where other behaviours (e.g., resting) would fall along that continuum [72]. Thus, state-space models with discrete state space such as HMMs may be easier to interpret, as each state has defined respective turning angle and step length distributions [72].

Hidden Markov models are well-suited for high-frequency movement datasets and behaviour-specific research questions since they provide an approach to identify habitat features relevant to separate behaviours [67]. For example, HMMs have been used to quantify how foraging probability relates to canopy cover [73], how travel speed relates to snow depth [16], and how animal orientation relates to wind direction [68]. Further, HMMs have been used to determine how variation in human influence on the landscape affects movement mode [74].

Choosing a model

By design, RSFs are well-suited to address large-scale questions on important areas for species (e.g., second-order selection; at the home range scale), whereas SSFs and HMMs are typically well-suited to address smaller-scale movement- or behaviour-specific questions (Box 1, e.g., third-order selection; at the habitat usage scale; [17–19, 75]). Since RSFs often generate the availability sample within the home range (MCP) of the observed data, and SSFs utilize a movement kernel that constrains the model to the scale of a step [53], the SSF prediction coefficients provide insight at a smaller scale than in the large-scale (home range) case of RSFs (although see [76–78] for guidance on scaling up from SSFs). In contrast, HMMs classify movement into discrete behavioural states, and thus do not consider habitat covariates as a driver of habitat selection, rather, as a mechanism related to estimated behavioural states [61, 62].

Box 1. Which model to choose?

Decision tree for choosing an appropriate modeling approach for understanding animal relationships with environmental conditions using movement data, including if you are interested in how the environmental conditions influence animal space use or movement, and if the data are coarse- or fine-scale. For example, caribou movement in the winter could be related to snow depth. RSFs can be used for broad-scale inference relating caribou occupancy and snow depth, SSFs can be used to understand whether snow depth affects selection at finer scale, including if it slows down movement speed, and HMMs can identify the depths of snow associated with a behaviour (e.g., foraging)

Important considerations

Choosing an appropriate model can be difficult, and choosing an adequate model does not guarantee appropriate inference. For example, [79] used HMMs to identify white tailed deer (Odocoileus virginianus) male–female interaction events, but upon validation, found that the model incorrectly identified interactions, and its use would mislead researchers about when and where mating behaviours are occurring. Additionally, [80] used move-persistence models to identify relationships between foraging behaviour and prey abundance, and found mis-matches between foraging behaviour, prey biomass, and diving behaviour. These cautionary tales highlight the importance of understanding the chosen models’ design, abilities, and limitations (Box 2). Where possible, validation exercises (e.g., leave-one-individual-out cross-validation, [81]) and incorporation of additional data streams (e.g., dive data, [71]) is recommended.

Box 2 Model comparison chart

Selection functions are subject to several assumptions. Selection functions assume that animals make resource-driven choices where at each location they choose the best habitat to move to next, but this may not always be the case (e.g., as in migration, or when relying on memory but see [84]). RSFs assume that locations are independent, and methods have been proposed to circumvent the violation of this assumption (e.g., data thinning [21], weighted approach [46]). Similarly, SSFs assume that the step length and turning angle are independent, yet autocorrelation may exist depending on the scale (e.g., animals may keep foraging in nearby locations; [63, 85]; but see examples to account for spatial and temporal structure: [86–88]). Thus, sampling intervals need to be examined because small time intervals can lead to misleadingly narrow confidence intervals [46, 53]. Additionally, SSFs assume independence between step length and turning angle, yet cross-correlation between the two may be unavoidable in some animal movement datasets [63]. Solutions to address this issue have been proposed [85, 89]. While many RSF and SSF models assume that individuals exhibit identical movement and selection, individual variation can be included in selection functions via random slopes in a mixed model (e.g., using the glmmTMB or mgcv package in R, [90, 91]) or by using a two-step approach of fitting individual models and averaging parameter estimates for inference at the population level (e.g., using the TwoStepClogit package in R, [44, 56]).

Choosing the number and the spatial extent of the availability sample are important computational and ecological considerations, respectively, and can affect the coefficient estimates and subsequently the ecological interpretation [29, 51, 55]. The number of available samples can be based on a constant (e.g., 10,000, [27, 33, 92]), point density (e.g., 1 available location per km², [93]), or by a ratio of observed:available locations (e.g., 1:20, [29]). Nevertheless, it is advisable to apply the selection function using a range of available sample sizes and choosing the sample size where the coefficient estimates stabilize [26, 29, 37, 92]. As approximation accuracy increases with the size of the available sample, it is generally better to include a larger availability sample, especially given quick computation time [29, 55, 91]. For RSFs, the spatial extent of the availability sample represents the home range of the animal, which is often represented as the MCP (e.g., absolute, buffered, or clipped) of the observed sample and can affect the scale of inference and thus should be chosen after consideration of the species’ ecology (see [94, 95]).

Movement HMMs rely on different assumptions that relax many assumptions of selection functions (RSF, SSF). Unlike selection functions, movement HMMs consider habitat covariates as related to estimated behavioural states, rather than as a driver of habitat selection. Nonetheless, HMMs acknowledge that telemetry data encompass distinct, unobserved behaviours and distinguish the state process from the observation process. HMMs also acknowledge that telemetry data are not independent, and explicitly model autocorrelation through the state process. Movement patterns and habitat associations are behaviour-specific, which can be modeled by the state-dependent distribution of the observation process and state transition probabilities.

These models (RSFs, SSFs, and HMMs) can be easily extended to incorporate inter- or intra-specific relationships. For example, interspecific relationships can be included by adding another species’ probability of occupancy as a covariate [28, 91, 96], or by using an individual-based model [97]. Additionally, more complex and dynamic intraspecific relationships can be included by considering individual animal’s movement relative to the group movement (see [98] for an HMM example).

Maps of probability of use are most easily generated using RSFs [24] and are particularly useful for conservation management and land-use planning [39, 99, 100]. Maps can be made using SSFs and HMMs [12, 13, 76–78], however for SSFs, require path simulation exercises with setting choices that are subjective (e.g., choosing the number of tracks and steps in the simulation, the burn-in size). Recent work has offered guidance on space use predictions from SSFs, including implementation in amt [12, 77, 78] (and see our step-by-step tutorial in Online Appendix 4).

Case study: ringed seals in Hudson Bay, Canada

Background and methods

We used movement data of one ringed seal (Pusa hispida) to conduct a comparative analysis of different methods to characterize their relationship with one covariate. We chose prey diversity as the covariate of interest for exploration (data from [83]), given that ringed seals are known as opportunistic predators (Online Appendix 2; [101]). Providing an accurate description of how ringed seals interact with their environment is important because of their ecological and cultural significance. Ringed seals are facing demographic decline in some areas of their range [102] underscoring the urgency to identify important areas of occupancy for this species. In this case study, we apply each of the models to the same dataset to demonstrate and compare the insight gleaned from each. The interpretation of this case study is limited, and focuses on illustrating how the inference will be affected by which model is applied.

We analyzed the estimated movement of a seal (Fig. 4A) equipped with an Argos satellite telemetry transmitter from 29 Oct 2012 to 17 Mar 2013 in Hudson Bay, Canada (see [80] for details). We used a correlated random walk state-space model (aniMotum R package, [103]) to filter and regularize the location data at a 24-h time step since Argos location data is observed irregularly in time and is prone to error [104]. This procedure resulted in one location per day from 29 Oct 2012 to 17 Mar 2013 (total of 140 observed locations). Then, we matched each seal location with prey diversity in the corresponding grid cell (Fig. 4D).

Fig. 4

The observed (used) and available locations for the A resource selection function (RSF) and B step-selection functions (SSF). Purple circles represent the observed locations and purple lines denote the path between observed points. Peach circles represent the available locations for the RSF and control locations for the SSF (see Fig. 2). The RSF availability sample was generated randomly within the minimum convex polygon of the track. C Decoded states from the hidden Markov model (HMM), including slow movement (green squares), moderate movement (blue circles), and fast movement (purple triangles). D Ringed seal estimated locations (yellow circles) and straight-line path between locations (yellow line) overlaid on estimated prey diversity

We used the prey diversity data as a covariate in two RSFs (one on the full dataset and one on a thinned dataset), two SSFs (one with and one without movement-habitat interactions), and an HMM (in the transition probabilities, but see Online Appendix 4 for how to apply covariates to the observation probabilities). For the RSF, we included ten available locations per observed location that were sampled from within the MCP of the observed locations. We used the amt package in R [12] to fit the RSFs and SSFs. Specifically, we estimated the RSF coefficients associated with the prey diversity (prey) covariate by fitting a logistic regression to the use-available observations $y = \left( {y_{1} , ... ,y_{n} } \right)$, where $y_{i} \in \left\{ {0,1} \right\}$. We used a logit link function to keep the probability $Pr(y_{i} = 1|x_{preydiv,i} )$ between 0 and 1:

$$Pr\left( {y_{i} = 1{|}x_{prey,i} } \right) = \frac{{{\text{exp}}\left( {\beta_{0} + \beta_{prey} x_{prey,i} } \right)}}{{1 + {\text{exp}}\left( {\beta_{0} + \beta_{prey} x_{prey,i} } \right)}} .$$

(13)

By default, the fit_rsf() function in the amt package calls glm()and thus fits logistic models with an intercept, however, in the context of an RSF based on use-available data, this intercept is not interpretable and is disregarded. Therefore, the RSF derived from this model excludes the intercept:

$$w\left( {x_{prey,i} } \right) = \beta_{prey} x_{prey,i} .$$

(14)

Since the full dataset was autocorrelated (see tutorial) and independent data is required for RSFs (see Box 1), we also fit an RSF on a thinned dataset comprising of every 10th location in the full dataset, using the same equations as the RSF on the full dataset.

We estimated the SSF coefficients by fitting a conditional logistic regression. At each time t, we sampled ten control locations $s_{t,i}$ for $i \in \left\{ {1, \ldots ,10} \right\}$ from a movement kernel $\phi$ fitted to the observed data and appended them to the observed location ${s}_{t,0}$. Then, the conditional probability that the point at location ${s}_{t,i}$ is the observed location given the environmental covariates at that location ${\mathbf{x}}_{t,i}$ and the animal’s previous locations ${\mathbf{s}}_{1:t-\text{1,0}}$ [53, 55] is:

$$Pr\left({y}_{t,i}=1\right|{{\mathbf{X}}_{t},\mathbf{s}}_{t,0:m},{\mathbf{s}}_{1:t-\text{1,0}}) =\frac{\text{exp}({\eta }_{t,i}) }{{\sum }_{j=0}^{10}\text{exp}({\eta }_{t,j}) } ,$$

(15)

where ${y}_{t,i}$ is a binary variable that takes the value of 1 when the i^th point is the observed location and 0 for the control locations and the linear predictor ${\eta }_{t,i}$ includes prey diversity, step length, turning angle as covariates (see Online Appendix 2):

$${\eta }_{t,i}={\beta }_{1}{x}_{prey, t,i}+{\beta }_{l}{l}_{t,i}+{\beta }_{\text{ln}}\text{ln}({l}_{t,i})+{\beta }_{\theta }cos({\theta }_{t,i}) .$$

(16)

We also fit a SSF where we allowed prey diversity to affect the movement kernel (specifically, the shape of the gamma distribution), which was achieved by including an interaction between prey diversity and the natural log of the step length:

$${\eta }_{t,i}={\beta }_{1}{x}_{prey, t,i}+{\beta }_{l}{l}_{t,i}+{\beta }_{\text{ln}}\text{ln}({l}_{t,i})+{\beta }_{\theta }cos({\theta }_{t,i})+ {\beta }_{2}{x}_{prey,t,i}\text{ln}({l}_{t,i}) .$$

(17)

Note that the fit_ssf() function in the amt package calls survival::clogit()and thus fits conditional logistic models without an intercept. Further, we chose to include the natural log (to model the gamma distribution) of step length and the cosine (to model the von Mises distribution) of turning angle as covariates due to their conventional use in the literature and model guides [12, 53, 105].

Finally, we fit a three-state HMM using the momentuHMM package in R [13], with prey diversity as a covariate on the transition probability:

$${\gamma }_{t,i,j}=\frac{\text{exp}({\beta }_{0,i,j}+{\beta }_{1,i,j}{x}_{prey,t })}{{\sum }_{l=1}^{N}\text{exp}({\beta }_{0,i,l}+{\beta }_{1,i,l}{x}_{prey,t })} ,$$

(18)

where $N$ is the number of behavioural states (set to $N = 3$), and $\beta_{0,i,j}$ and are the intercept and coefficient quantifying the effect of prey diversity on transition probability from state to, respectively. We fixed ${\beta }_{0,i,j} ={\beta }_{1,i,j}=0$ when $i=j$ (i.e., the diagonal of the transition probability matrix).

We predicted the probability of use relative to prey diversity (log of the relative selection strength [log-RSS] for selection functions; see [106]) and in space (steady-state utilization distributions [SSUDs]; [76]). Detailed methods are available in Online Appendix 2 and the data and a coding tutorial of the analysis are available in Online Appendix 3 and 4, respectively, and on Github (https://github.com/kflorko/movementstats_review).

Results and discussion

These models led to important differences in the magnitude, and sometimes direction, of the relationships with the same covariate, and thus resulted in different ecological inference. The estimated selection coefficient that represented the habitat selection relationship with prey diversity was significantly positive for the RSF on the full dataset, but not significantly different from zero for the RSF on the thinned dataset or for the SSFs (Table S1, Fig. 5A–B). Notably, the estimated coefficient for prey diversity from RSF on the thinned dataset, while not significant, was similar to that from the RSF on the full dataset, and both were much higher than that from the SSF (Table S1). The SSF with movement-habitat interactions showed no significant relationship between prey diversity and movement speed (Table S1), and as a result the relationship between selection and prey diversity was not significantly modified (Fig. 5B). The HMM characterized the states mostly on movement speed (i.e., step length, Fig. S1), and the states had contrasting relationships with prey diversity (Fig. 5C). The probability of being in the slow behaviour, which is often of higher importance for conservation because it is thought to be associated with foraging ([107], but see [80] for a discussion on the complexities of interpretation), had a positive relationship with prey diversity (Fig. 5C). The probability of being in the moderate speed behaviour decreased with prey diversity, while the probability of being in fast behaviour remained relatively constant with prey diversity, with a potential peak in low-medium prey diversity (Fig. 5C).

Fig. 5

Log relative selection strength (log-RSS) from the A RSF on the full dataset, RSF on the thinned dataset (every 10th location), and SSF without movement-habitat interactions and B SSF with movement-habitat interactions, where prey diversity affected step length; log-RSS was calculated for speeds at the 25th (slow), 50th (moderate), and 75th (fast) percentiles of step length, which we infer as “speed” since the locations are recorded at a fixed time interval (see [26, 106]). C Estimated state probabilities from the HMM as a function of prey diversity (on the transition probability). Shaded areas represent the standard error. Note the bowtie shape of the standard error in A–B is due to log-RSS calculating selection strength relative to a starting step, in this case, where prey diversity is set to its mean

If these results were being used for conservation, the HMM would suggest that high prey diversity is important for ringed seal space use and important behaviours, while the RSF on the thinned dataset and the SSFs would suggest that this covariate is not important, likely due to the conceptual attributes of each model. For example, the different selection coefficient estimates between the RSF and the SSFs may be attributed to their distinct definition of available habitat. The RSF assumes the entire study area (within the minimum convex polygon (MCP) in this case study) is accessible to the animal, whereas the SSF generates control points at each step and availability is estimated simultaneously with habitat selection through the movement kernel, rather than being assumed a priori (Fig. 4A–B). Further, it is expected that RSFs will have greater coefficients than SSFs due to scale [108, 109]. It is likely that the positive relationship with prey diversity for the RSF on the full dataset is significant due to autocorrelation in the dataset rather than the scale at which the model assesses habitat-environment relationships. Indeed, once the dataset was thinned, the estimate—while still positive and with a similar effect size—was not significant (Table S1, Fig. 5A), highlighting the importance of considering autocorrelation in the dataset (Box 1). The RSF on the full dataset and HMM slow movement state likely show similar results because slow movement is associated with greater residency (i.e., remaining in an area) and thus will result in more locations in an area.

Differences in coefficient estimates will affect maps of predicted probability of seal occurrence or behaviour (selection function or HMM, respectively; Fig. 6). Given the model’s positive coefficient estimate with prey diversity, both RSFs predicted that seals are most often found in a small portion of the study area where there is high prey diversity (Fig. 6A shows the RSF on the thinned dataset; the RSF of the full dataset renders a similar map as the coefficient is similar, Table S1). The SSFs did not have significant relationships with prey diversity, thus, the maps generated from the SSUDs showed that probability of use was more uniformly distributed in space (Fig. 6B–C). The HMM will produce one map per state (behaviour), and here with vastly contrasting spatial patterns. The slow movement state, which was positively related to prey diversity, showed a similar restricted spatial pattern to the RSF (Fig. 6D), whereas the moderate movement state, which had a negative relationship with prey diversity, showed the opposite pattern to the RSF prediction map (Fig. 6E). The fast movement state did not show much spatial variation (Fig. 6F), reflective of its non-significant relationship with prey diversity (Fig. 5F). The SSF-generated maps did not show any similarities with the other models’ maps (or relative to prey diversity), likely due to prey diversity not being relevant at the step-scale. Overall, we see that prey diversity may be relevant for the behavioural state transition at the whole track-level scale (HMM). We also note important considerations for the autocorrelation in the data, as once thinned, the RSF no longer estimated a significant relationship for prey diversity. This case study demonstrates how choosing the wrong model for the question at hand can mislead ecological insight and could misinform conservation and management decisions.

Fig. 6

Predicted density (A–C) and predicted state probability (D–F; herein prediction) from A resource selection function (RSF; thinned dataset), B step selection function (SSF) without movement-habitat interactions, C SSF with movement-habitat interactions, and hidden Markov model (HMM) D state 1: slow movement, E state 2: moderate movement, and F state 3: fast movement. Note that utilization distributions in B and C were created using simulated paths with 200,000 locations

In this case study, selection at the home range scale modeling approach (i.e., RSF) may be most relevant due to the coarse resolution of the prey data. To simplify comparisons across models, this case study was limited to one seal, one prey metric, and used the default settings for these models. Thus, to reflect realistic ecological dynamics, an analysis of a more comprehensive ringed seal dataset (e.g., more individuals, longer tracking, and more prey metrics) should be performed in the future, and further, the uncertainty associated with Argos data and in generating prediction maps should be accounted for.

Future directions

Selection functions and HMMs are widely applied models that provide essential insight from animal movement data, and several new and extended versions of these models offer deeper ecological understanding. For example, recent models have incorporated memory [110], diel variation [111], and habitat connectivity [112] to better understand animal movement. Further, integration of habitat availability into HMMs (termed HMM-SSF) has been developed to provide insight on behaviour-specific habitat selection (e.g., [113–115]). Additionally, new approaches incorporate additional data streams, such as those from spatial surveys [116]. These advancements extend conventional animal-habitat models for movement data and exemplify the recent exciting progress in the field.

Recent work has also addressed limitations of mainstream selection functions and HMMs by improving model fitting, checking, and simulation capabilities. For example, models have incorporated more complex spatiotemporal structures [117–119], non-linear and random effects [90], and cross-correlated movement kernels [85]. Notably, to support the growing complexity of these models, the recent development of the hmmTMB package enhances the implementation of HMMs by allowing for more complex models and inclusion of random effects [120]. Additionally, new methods and tools have been introduced to evaluate model goodness-of-fit [121], address challenges associated with incomplete data [122], and ease simulation for space use predictions from SSFs [77, 78].

There are a variety of machine learning approaches (e.g., MaxEnt, random forests) that we have not covered because our review is focused on statistical methods. However, it is important to note that MaxEnt is considered equivalent to an exponential RSF [92, 123], and generally may be better suited for occupancy data (i.e., a location can only be used once, as in species distribution data, [124]), but can also be applied to movement data [27]. Additionally, random forests may be used for modeling animal movement and generally relax the assumptions associated with parametric methods, but the results may be more challenging to interpret [125]. If prediction is the goal, machine learning approaches including MaxEnt and random forests, or a parametric approach using least absolute shrinkage and selection operator (LASSO) for model selection, are often recommended [38, 100].

Conclusions

The models reviewed in this paper for linking animals’ movement data to the environment are widely applicable and can be easily implemented but provide different ecological insights. Our review highlights the conceptual differences that will aid ecologists in choosing a model based on their goals. For example, RSFs may be well suited for identifying broad corridors or protected areas, SSFs for understanding movement patterns, and HMM for understanding behaviours. Particularly, we highlight that RSFs are the most broad-scale model and are useful for home-range studies, whereas SSFs and ultimately HMMs get increasingly more movement-specific. Also, we show how RSF results change when autocorrelation is removed. Further, the sampling—or exclusion as in the case of HMMs—of availability locations can influence the coefficient estimates, as can the inclusion of movement information in the model. All these models are invaluable for ecological research, yet careful consideration of the mathematical and conceptual underpinnings is necessary to address the ecological questions and applied conservation management that these results may inform.