Linking resource selection and step selection models for habitat preferences in animals

The research community working on species-habitat associations in animals is currently facing a paralysing methodological conundrum, because its two dominant analytical approaches have been shown to reach divergent conclusions. Models fitted from the viewpoint of an individual (step selection functions), once scaled up, do not agree with models fitted from a population viewpoint (resource selection functions). We explain this fundamental incompatibility, and propose a solution by introducing to the animal movement field a novel use for the well-known family of Markov chain Monte Carlo (MCMC) algorithms. By design, the step selection rules of MCMC lead to a steady-state distribution that coincides with a given underlying function: the target distribution. We therefore propose an analogy between the movements of an animal and the movements of an MCMC sampler, to guarantee convergence of the step selection rules to the parameters underlying the population's utilisation distribution. We introduce a rejection-free MCMC algorithm, the local Gibbs sampler, that better resembles real animal movement, and discuss the wide range of biological assumptions that it can accommodate. We illustrate our method with simulations on a known utilisation distribution, and show theoretically and empirically that locations simulated from the local Gibbs sampler arise from the correct resource selection function.


INTRODUCTION
Understanding how animals use a landscape in response to its habitat composition is a crucial question in pure and applied ecology. Such insights are achievable only by confronting species-habitat association models with usage data, collected either via transect surveys or via biologging methods. Statistical inference, to link these data to environmental variables, can be approached from a population perspective, using resource selection functions (RSF; Manly et al. 2002). Alternatively, if individually referenced data (i.e., telemetry) are available, the question can be addressed from the viewpoint of the single animal, via step selection functions (SSF; Thurfjell et al. 2014). The population/individual dichotomy between these two approaches is not always clear-cut, because RSFs can be applied to the utilization distribution of single animals, and SSFs can combine joint insights from multiple individuals. Nevertheless, the two methods roughly fall at opposite ends of the Eulerian-Lagrangian spectrum outlined by Turchin (1998). Therefore, researchers in this area have tended to think of the habitat preference parameters obtained via SSFs as the microscopic rules of movement, while the corresponding parameters of an RSF are implicitly thought of as the macroscopic patterns obtained in the long term. Hence, SSF models are increasingly concerned with the geometry of movement trajectories (e.g., step lengths and turning angles in different behavioural states in Squires et al. [2013]), while RSF predictions often make a pseudo-equilibrium assumption (Guisan and Thuiller 2005), which is a biological term reminiscent of the mathematical idea of steady-state distributions. But herein lies a fundamental problem for this entire field of statistical analysis. A correctly formulated framework of movement must work across scales, such that, when the microscopic rules of individual movement are scaled up in space and time, they give rise to the expected macroscopic distribution of a population. However, there is now both analytical (Barnett and  and numerical (Signer et al. 2017) evidence that the distribution constructed from the coefficients of a SSF does not match the spatial predictions of the RSF fitted to the same data. Here, we explain how this discrepancy arises and propose a solution.
A RSF w(c) is proportional to the probability of a unit of habitat c being used (Boyce and McDonald 1999).
Depending on the type of usage data available, RSFs are derived in two steps. First, a model is fitted to the response and explanatory data. For example, a point process model (Aarts et al. 2012) or a use-availability logistic regression McDonald 1999, Aarts et al. 2008) can be used for telemetry data, and a log-linear regression can be used on count data from regular grids or line transects. Second, irrespective of the type of response data and model fitting method, the linear predictor of the resulting statistical model is transformed via a non-negative function (Manly et al. 2002: Chapter 2), of which the most common is the exponential where c is a vector of m covariate values, and b 1 , b 2 , . . . , b m are the associated regression coefficients. The RSF can be used to model the utilization distribution pðxÞ, i.e., the distribution of the animal's space use where the functions c 1 , c 2 , . . ., c m associate a spatial location x to the corresponding covariate values, and X is the study region. The utilization distribution is normalized by the denominator in Eq. 2 to ensure that it defines a valid probability distribution for x, hence the lack of an intercept in the linear predictor. Although they can encompass a wider range of environmental conditions, the covariates are often called resources in this context. In the following, we use "covariates" and "resources" interchangeably. Resource selection function approaches are commonly used to estimate the apparent effect of a spatial covariate on a species. The resource selection coefficients b k characterize this effect for each of the m covariates (b k > 0: preference; b k < 0: avoidance; b k = 0: indifference; see Avgar et al. [2017] for a discussion of the interpretation of the b k in terms of selection strength). However, recent work has shown that these interpretations are highly sensitive to the context in which the organisms are being studied, in particular, the availability of all habitat types to the animals (Beyer et al. 2010, Matthiopoulos et al. 2011, Paton and Matthiopoulos 2016. Thus, in this framework, the definition of habitat availability, determined by assumptions of spatial accessibility (Matthiopoulos 2003), is important in deducing preference from observed usage. For example, when using RSFs to analyze a time series of positions from a ranging animal, it may not be plausible to assume that all locations in the home range are accessible by the animal at every step (Northrup et al. 2013). Resource selection function approaches are often forced to treat such non-independence as a statistical nuisance (Aarts et al. 2008, Fieberg et al. 2010, Johnson et al. 2013), but step selection approaches treat it as an asset.
In step selection analyzes, the likelihood p(y|x) of a potential displacement by the animal to a location y over a given time interval (typically, the sampling interval) is modeled in terms of the habitat composition in the neighborhood of the animal's current position x where /ðÁjxÞ is defined over a spatial domain X, and, for any location x, c(x) = (c 1 (x), c 2 (x), . . ., c m (x)). The function /ðÁjxÞ is called the resource-independent movement kernel around x (Rhodes et al. 2005, Forester et al. 2009), and it describes the density of endpoints for a step starting in x, in the absence of resource selection. To link the movement to environmental covariates, w is modeled using the same log-linear link as the RSF, given in Eq. 1. In this context, the term "step selection function" is most often used for w (e.g., by Fortin et al. 2005, Thurfjell et al. 2014; however, note that it is sometimes used for the whole numerator in the right-hand side of Eq. 3 (see Forester et al. 2009). In the following, we call w the SSF.
The choice of the function / characterizes accessibility, and hence determines availability, in a step selection model; it corresponds to the distribution of feasible steps over one time interval, with origin x, when the resources do not affect the movement. It can, for example, be a uniform distribution on a disc around the current location x (e.g., Arthur et al. 1996), or obtained from the empirical distributions of movement metrics (e.g., step lengths and turning angles in Fortin et al. [2005]).
Step selection functions are most often fitted using conditional logistic regression on matched use-availability data, where each observed step x t ! x tþ1 is matched to a set of random steps generated from /ðÁjx t Þ (Thurfjell et al. 2014). Duchesne et al. (2015) showed that a step selection model defines a movement model equivalent to a biased correlated random walk. Biased correlated random walks are routinely used in ecology as a flexible basis for models of individual movement (Turchin 1998, Codling et al. 2008. Avgar et al. (2016) extended the step selection approach to allow simultaneous inference on habitat selection and on the movement process, making it a very attractive framework to estimate habitat preference from movement data (Prokopenko et al. 2017, Scrafford et al. 2018. Step selection models have been used to analyze the impact of landscape features on animal space use (e.g., Coulon et al. 2008, Roever et al. 2010, as well as animal interactions (Potts et al. 2014b).
Although the RSF and SSF are typically described with the same notation, and used for the same purpose of estimating habitat preference, it can be shown that their steady-state predictions do not generally coincide. For a known utilization distribution, Signer et al. (2017) showed empirically that the normalized SSF ("naive" estimate) differed from the utilization distribution. In particular, the difference was greater when / was narrow compared to the scale of habitat features. Similarly, Barnett and Moorcroft (2008) showed that, for the step selection model defined in Eq. 3, the steady-state distribution of the animal's location (i.e., its utilization distribution) is given by That is, the steady-state distribution of the model is generally not proportional to the SSF w, and that discrepancy crucially depends on the choice of the resource-independent movement kernel /. An example of this is their earlier result ) that under one specific set of assumptions, the steady-state distribution is approximately proportional to the square of the SSF.
Although it may seem disconcerting that the two approaches lead to different estimates of w, the cause of this apparent paradox is partly due to the notational misuse of the same symbol for what are, in effect, different objects. The SSF captures local aspects of the animal's movement, because it only considers a neighborhood of the current location of the animal (determined by /) and only becomes a better approximation of the RSF when the scale of / increases . The parameters of the two objects coincide in the limiting case of unconstrained mobility, i.e., when the availability assumed by both methods is global. However, in every other case, the two methods are different. Schl€ agel and Lewis (2016) also noted that, unlike RSF models, standard SSFs are scale dependent, in that their habitat selection estimates depend on the time scale of the observations (although see Hooten et al. [2014] for a SSF approach with a user-defined scale of selection).
Several approaches have been suggested to approximate the steady-state distribution of SSF movement models. In particular, Avgar et al. (2016) and Signer et al. (2017) showed that simulations from a fitted SSF could be used to obtain estimates of the underlying utilization distribution. Similarly, Potts et al. (2014a) described a numerical method to compute the utilization distribution given in Eq. 4, as it generally has no closed form expression. Those approaches are useful to predict space use from SSFs, but they do not allow the steadystate distribution of locations to be modeled in a simple parametric form, as in Eq. 2. One important consequence is that, because the utilization distribution of SSF models is not modeled by a RSF, joint inference from telemetry data and survey data into habitat selection and space use has not been possible with existing approaches.
Rather than seeking an equivalence of the parameters estimated by RSF and SSF methods, a better question to ask is: under what assumptions do the parameters estimated by a SSF lead to movement that scales to the distribution yielded by the parameters of a RSF model? In A model of step selection using a movement-MCMC analogy, we reconcile resource selection and step selection conceptually, with a new step selection model for which the longterm distribution of locations is guaranteed to be proportional to the RSF. Our method uses an analogy between the movement of an animal in geographical space and the movement of a Markov chain Monte Carlo (MCMC) sampler in its parameter space. In The local Gibbs sampler, we make these concepts applicable in practice, by developing a family of MCMC algorithms with considerable potential for encompassing realistic movement assumptions. In Simulations, we illustrate our method using simulations on a known utilization distribution. We verify that the distribution of simulated locations corresponds to the correct RSF, and we present a proof-of-concept analysis to demonstrate the potential of the method for estimating resource selection coefficients and parameters of the movement process from telemetry data.

A MODEL OF STEP SELECTION USING A MOVEMENT-MCMC ANALOGY
Markov chain Monte Carlo methods are a general framework to sample from a probability distribution, termed the target distribution (Gilks et al. 1995). This approach is mostly used for Bayesian inference, to sample from the (posterior) distribution of a set of unknown parameters (Gelman et al. 2014: Chapter 11). It includes a very wide class of algorithms, among them the widely used Metropolis-Hastings and Gibbs samplers. An MCMC algorithm describes the steps to generate a sequence of points x 1 , x 2 , x 3 . . ., whose long-term distribution is the target distribution. Each MCMC algorithm is defined by its transition kernel p(x t+1 |x t ), which determines (for any t = 1, 2, . . .) how the point x t+1 should be sampled, given x t . For example, in a Metropolis-Hastings algorithm, the transition kernel is a combination of the proposal distribution and the acceptance probability In general, given some easily satisfied technical conditions, a sufficient condition for pðx tþ1 jx t Þ to define a valid MCMC algorithm for the target distribution p (i.e., to ensure that the distribution of samples will converge to p) is the detailed balance condition 8x; y; pðyÞpðxjyÞ ¼ pðxÞpðyjxÞ: That is, if the process is in equilibrium with distribution p, then the rates of moves in each direction between any x and y balance out.
We propose an analogy between an animal's observed movement in n-dimensional geographical space, and the movement of an MCMC sampler in a n-dimensional parameter space, for which the target distribution is the utilization distribution. That is, we consider that a tracked animal "samples" spatial locations in the short term from some movement model and, in the long run, from its utilization distribution, in the same way that an MCMC algorithm samples points in the short term from some transition kernel and in the long term from its target distribution. An MCMC algorithm then defines a movement model, for which the steady-state distribution is known. The dynamics of the movement process x t are described by the transition kernel of the algorithm such that, at each time point t = 1,2,. . ., the next location x t+1 is sampled from p(x t+1 |x t ). By the properties of MCMC samplers, the steady-state distribution for x t is p. The utilization distribution can be modeled with the RSF, as defined in Eq. 2, to link the target distribution of the movement model to the distribution of resources.
An MCMC algorithm, if viewed as a movement model, can then be used to analyze animal tracking data, in the following steps. Although we focus on step 1 in this paper, we illustrate steps 2 and 3 with a simulated example in Local Gibbs estimation.
1) Choose an MCMC algorithm, to be used as a model of animal movement and habitat selection. We suggest one such algorithm in The local Gibbs sampler. 2) Write the likelihood of the model. Under an MCMC movement model, the likelihood of an observed step from x t to x t+1 is a function of the resource selection coefficients and of the other parameters of the sampler, given by the transition kernel p(x t+1 |x t ). 3) Use maximum likelihood estimation, or other likelihood-based methods, to estimate the resource selection and movement parameters.
In this framework, the choice of the MCMC algorithm determines the movement model. For example, with a Metropolis-Hastings model, different proposal distributions might capture different features of the animal's movement. The parameters of the algorithm, which are usually regarded as tuning parameters, are here parameters of the movement process. For example, the variance of the proposal distribution can be thought of as a measure of the animal's speed. It is important to make a distinction between these parameters of movement, and the parameters of the target distribution (i.e., the resource selection parameters). Two different samplers might have the same target distribution, but the rate at which it is approached by the MCMC samples will depend on the choice of algorithm. Indeed, part of the success of MCMC in its Bayesian context is the flexibility in choosing the transition kernel for a given target distribution. The suitability of an MCMC sampler is usually assessed by the speed of convergence of the simulated samples to the target distribution. However, for our application, we want an algorithm corresponding to a realistic model of movement, in addition to having the correct target distribution. It could happen that an MCMC algorithm that describes animal movement very realistically has a slow rate of convergence to the target distribution. This would merely mean that the animal, when observed at the time step of the observations, does not sample efficiently from its utilization distribution. In such a case, inference about the utilization distribution would be limited regardless of the modelling framework that is used.
In rejection-based MCMC algorithms such as Metropolis-Hastings, a relocation is proposed at each time step, and is accepted with some probability. If the proposed step is not accepted, the process remains in the same location. Although it can happen that a tagged animal is immobile over several time steps (in particular if temporal resolution is high), many telemetry data sets do not include such "rejections." Classic MCMC algorithms might thus seem to be an unnatural choice to analyze those data, because the animal will almost always change position in the process of sampling a new candidate location. To circumvent this problem, we design a new rejection-free MCMC algorithm in The local Gibbs sampler.

THE LOCAL GIBBS SAMPLER
Standard Metropolis-Hastings samplers require a rejection step to ensure convergence to the target distribution. Viewing this as a movement model would imply the unlikely scenario of a return by the animal to its previous position, after having tested and rejected a relocation. Instead, it is more natural to think about tracking data as the outcome of a rejection-free sampler. Several such algorithms are possible; see Discussion. Here, we describe one such algorithm, which we call the local Gibbs sampler.
In the classic Gibbs sampler, each "step" involves updating just one of the n parameters, x j say, while keeping x 1 ,. . ., x jÀ1 , x j+1 , . . ., x n fixed; the values of j can be chosen systematically or randomly. Thus, each step is a move within a one-dimensional subspace of the parameter space, rather than over the whole space. It is used when the target distribution over each such one-dimensional space (the so-called "full conditional distribution") is mathematically tractable, so that when it is used as the transition kernel for that step, the acceptance probability is guaranteed to be 1.
The local Gibbs sampler uses the same idea of sampling from a restricted part of the target distribution: at each iteration t, the updated parameter x t+1 is sampled directly from the target distribution, truncated to some neighborhood of x t . The way in which this neighborhood is selected is crucial to ensuring that the algorithm samples from the required target distribution in the long run.
In explaining the details of the algorithm, we focus on the case of n = 2 dimensions, by far the most important case for ecological applications, though the algorithm works for any n with straightforward changes. For any point x, and r > 0, we define D r ðxÞ to be the disc of centre x and radius r.
The local Gibbs sampler for p is given by the following steps, and the notation is illustrated in Fig. 1. The track starts from a location x 1 , and moves to locations x t+1 over iterations t = 1,2,. . .. 3) Sample the next location x t+1 from D r ðcÞ according to the constrained pdfp.
The local Gibbs sampler has one parameter: the radius r > 0 of the relocation disc. Here, for simplicity, we only consider the case where r is fixed, but the algorithm would still work if r were generated independently at each iteration from a probability distribution.
Using the analogy introduced in A model of step selection using a movement-MCMC analogy between animal movement and MCMC sampling, the local Gibbs algorithm can be used as the basis for a model of animal movement and habitat selection, that we will call the local Gibbs model. It relies on the assumption that an animal "samples" locations from its utilization distribution based on the step selection rules described above.
Note that, at each time step, the overall relocation region of the local Gibbs model is symmetric around the animal's current location. The choice of the relocation disc D r ðcÞ, based on the selection of a point c in step 1 of the algorithm, might seem biologically unrealistic, because a moving animal would not relocate to a disc that is shifted at random from its current location. Nevertheless, because c is chosen uniformly from D r ðx t Þ, one should think of the relocation region once c has been integrated over, i.e., a disc of radius 2r around x t .
In the local Gibbs model, the parameter r determines the size of the area that is available to the animal over one time step. As in most step selection analyzes, the region of availability is a simplistic but useful model for a combination of the animal's mobility and perception.
Taking p to be the normalized RSF (Eq. 2), the local Gibbs algorithm defines a step selection (movement) model in which the distribution of the animal's space use is guaranteed to be proportional to the RSF. Indeed, it satisfies the detailed balance condition (Eq. 5), which can be shown as follows. Given r, we have Given c, y is sampled from D r ðcÞ with a density proportional to pðyÞ and, given x, c is sampled uniformly from D r ðxÞ, so as required. The local Gibbs model is superficially similar to the availability radius model of Rhodes et al. (2005), first introduced by Arthur et al. (1996). In that model, at each time step, the next location x t+1 is sampled from the RSF truncated and scaled on a disk centered on x t . That is, in step 1 of the algorithm described above, they take c = x t . This means that there is no mechanism in their approach to guarantee that the overall distribution of the sampled locations is the RSF. Specifically, the two sides of the detailed balance equation involve different normalization constants, and so their movement models do not have the normalized RSF as their equilibrium distributions. For this reason, the coefficients they estimate will differ from the resource selection coefficients estimated from a RSF approach.
We can derive the resource-independent movement kernel / LG ðyjxÞ of the local Gibbs model, to describe the distribution of steps on a flat target distribution. In the case where r is fixed / LG ðyjxÞ ¼ 1 ðpr 2 Þ 2 AðD r ðxÞ \ D r ðyÞÞ ifky À xk 2r 0 otherwise 1. Notation for the local Gibbs sampler in two dimensions. The point c is sampled uniformly from D r ðx t Þ, and the next location x t+1 is sampled from the resource selection function (RSF) truncated to D r ðcÞ.
where ky À xk is the distance between x and y, and AðD r ðxÞ \ D r ðyÞÞ is the area of the intersection of the discs of centers x and y, and of radius r. The point c is such that kc À xk \ r and kc À yk \ r, and so, in the absence of environmental effects, the relative probability of a step from x to y is proportional to AðD r ðxÞ \ D r ðyÞÞ. By construction, it is impossible to have a step between two points if the distance between them is larger than 2r, hence / LG ðyjxÞ ¼ 0 when ky À xk [ 2r. The detail of the derivation is given in Appendix S1. A graph of the density function / LG is shown in Fig. 2.
The transition kernel given in Eq. 6 and plotted in Fig. 2 describes the distribution of steps in the absence of habitat selection, in the case where the radius parameter r is fixed. A more flexible movement model can be obtained by taking r to be time varying, and drawn at each time step from a probability distribution (e.g., exponential or gamma distribution, to ensure r > 0).
It is important to note that the transition kernel of the local Gibbs algorithm cannot be written in the form given in Eq. 3, i.e., p(y|x) is in general not proportional to / LG ðyjxÞwðcðyÞÞ. For this reason, the local Gibbs model is not merely a special case of the step selection model described by Forester et al. (2009).

SIMULATIONS
The local Gibbs algorithm, described in The local Gibbs sampler, can be used to simulate tracks based on a known RSF. The truncation of the RSF to the disc D r ðcÞ requires the calculation of the normalizing constant C r ðcÞ. It is not generally possible to derive it analytically, but Monte Carlo sampling can be used to approximate it. In practice, to sample from the truncated target distributionp, n d points are generated uniformly in D r ðcÞ, and x tþ1 is sampled from those points, with probabilities proportional to their RSF values. Simulation using the local Gibbs algorithm is illustrated in Fig. 3.
Here, we illustrate the method described in A model of step selection using a movement-MCMC analogy, with the local Gibbs sampler. In Local Gibbs simulation, we show that our algorithm can produce movement tracks on a known utilization distribution and, in Local Gibbs estimation, we illustrate the use of the local Gibbs movement model for the estimation of resource selection and movement parameters from simulated data. The R code used for the simulations is available in the supplementary material, as Data S1.

Simulated resources
To mimic the type of environmental data of a real case study, we simulated two covariate distributions c 1 and c 2 as Gaussian random fields on square cells of size 1, using the R package gstat (Pebesma 2004). We restricted the study region to X ¼ ½À15; 15 Â ½À15; 15, to ensure that the target distribution is integrable. Plots of c 1 and c 2 are shown in Fig. 4A,B. The utilization distribution was defined by with b 1 = À1 and b 2 = 4 (i.e., avoidance for c 1 and preference for c 2 ). A plot of the RSF is shown in Fig. 4C.

Local Gibbs simulation
In this section, we demonstrate that the local Gibbs algorithm, described in The local Gibbs sampler, can be used to sample from a given probability distribution. We considered the utilization distribution p defined in Simulated resources. To analyze the behavior of the local Gibbs sampler at different spatial scales, we ran three simulations, with three different values for the radius r of the movement kernel: r = 0.5, r = 2, and r = 8. The value of r affects the range of perception of the animal and, indirectly, its speed. For each r, 5 9 10 5 locations were simulated with the local Gibbs algorithm, starting from the point x 1 = (0,0). (Given the length of the simulated tracks, the choice of the starting point has only a minor impact on the overall distribution of sampled locations.) For comparison, we also illustrate the results of Barnett and Moorcroft (2008), that the steady-state distribution of a standard SSF model (p in Eq. 4) differs from the normalized SSF. We sampled a movement track from a step selection model with uniform sampling, as defined by Forester et al. (2009), that we denote SSF unif . We simulated 5 9 10 5 locations from SSF unif , as follows. We started from x 1 = (0,0). Then, at each time step t = 1,2,. . ., we generated 100 proposed locations y 1 , y 2 ,. . ., y 100 uniformly from a disc of radius r = 3 centred on x t . The next location x t+1 was sampled from the proposed locations, with each point y i having a probability to be picked proportional to p(y i ). That is, we use p as the (normalized) SSF to simulate from the uniform sampling model. Here, we chose r = 3 because it gave rise to approximately the same mean step length as the local Gibbs sampler with r = 2 (i.e., comparable speed of spatial exploration).
The first 300 steps of each simulated track, and the density of all simulated points, are shown in Fig. 5. The density of points simulated from the local Gibbs sampler (right column, first three plots) displays the same patterns as the true RSF (Fig. 4C). By contrast, the density of the locations obtained in the SSF unif simulation (right column, last plot) fails to capture many features of the landscape, as the process spends a disproportionate amount of time in areas of high values of w(x).
To compare the empirical distribution of simulated points to the distribution p used in the simulations, we plotted the (normalized) count of locations simulated in each grid cell against the corresponding value of p. The comparison is presented in Fig. 6. Alignment with the identity line indicates similarity between the empirical distribution and p. For the three local Gibbs simulations, the points align well with the identity line, in particular in the experiments with r = 2 and r = 8, in which the speed of spatial exploration is higher than when r = 0.5. This confirms that the local Gibbs algorithm can sample movement trajectories on a given target distribution. It defines a movement model for which the long-term distribution of locations is known. However, the plot for the SSF unif simulation reveals a clearly nonlinear relationship between the density of simulated points and the normalized SSF. This confirms the results of Barnett and Moorcroft (2008), Avgar et al. (2016), andSigner et al. (2017): the coefficients of a SSF do not measure the underlying steady-state distribution. (Note that SSF models may be used to estimate space use, with simulations, as in Avgar et al. [2016], but the parameters of the SSF only measure local habitat selection.) We illustrated how the local Gibbs sampler can generate movement tracks that converge in distribution to the underlying RSF.

Local Gibbs estimation
The approach introduced in A model of step selection using a movement-MCMC analogy shows great promise for the estimation of movement and resource selection parameters from observed animal movement data. Considering the MCMC algorithm as a movement model, it is in principle straightforward to express the likelihood of observed steps, given the parameters of the sampler (e.g., radius r in the local Gibbs model) and of the RSF (b 1 , b 2 , . . .). In cases where the transition kernel of the chosen sampler, p(x t+1 |x t ), can be calculated, the Simulation using a local Gibbs sampler, with radius parameter r = 0.5 (first row), r = 2 (second row), and r = 8 (third row); and simulation using a step selection function with uniform sampling (r = 3, fourth row). The left column displays the first 300 simulated steps, and the background color represents the utilization distribution (i.e., the normalized resource selection function (RSF); the RSF is given in Fig. 4C). The right column shows the density of the 5 9 10 5 simulated locations, i.e., the normalized counts.
In this section, we wish to demonstrate its practical application, with the example of the local Gibbs model. We simulated a track of T = 3,000 locations from the local Gibbs sampler (described by the algorithm in The local Gibbs sampler), with r = 2, on the RSF defined in Simulated resources. Then, similarly to a real analysis, we used the local Gibbs model to recover estimates of the RSF (i.e., of b 1 and b 2 ) and of r, from the (simulated) movement data and covariate rasters.
The likelihood of an observed track under the local Gibbs model is obtained as the product of the likelihoods of the individual steps The details of the derivation are given in Appendix S1. This likelihood is a function of the movement parameter r, and of the coefficients b i of the RSF (which appear in the expression of p). Maximum likelihood techniques can then be used to obtain parameter estimates. We implemented the likelihood function of Eq. 7, and used the numerical optimizer nlminb in R to  6. Results of the simulations. In each plot, the distribution of simulated points (on the x-axis) is compared to the distribution p used in the simulations (on the y-axis). In the local Gibbs simulations, p denotes the (normalized) resource selection function (RSF) and, in the uniform step selection function (SSF unif ) simulation, p denotes the (normalized) SSF. Each dot represents the value associated with one grid cell. The closer the dots are to the identity line, the more similar the empirical distribution is to p. In the local Gibbs simulations, the empirical distributions are very similar to the RSF; the similarity increases with r, because a larger radius leads to faster spatial exploration. For the SSF unif model, there is a clear discrepancy between the empirical distribution and the SSF, as predicted by Barnett and Moorcroft (2008: Eq. 4). get maximum likelihood estimates of b 1 , b 2 , and r. The results are summarized in Table 1. Fig. 7 shows a plot of the estimated utilization value of each grid cell against its true utilization value. If we denoteŵ i the estimated value of the RSF in cell i, its estimated utilization valuep i is derived aŝ In Fig. 7, the alignment of the dots with the identity line indicates that the estimated utilization distribution captures the shape of the true utilization distribution well. In addition, the parameter r of the movement process was successfully estimated (Table 1).
This example demonstrates how the method can be used to estimate resource selection and movement parameters from tracking data. In real applications, unlike with simulated data, the true form of the movement process would not be known, and additional work would be needed to assess the fit. We discuss this further in Discussion.

DISCUSSION
We have presented a versatile class of models of animal movement, for which the steady-state distribution of locations is proportional to the same RSF that influences short-term movement. Our approach reconciles the resource selection and step selection approaches to the analysis of space use data. We anticipate that the resolution of this discrepancy between RSF and SSF models will have important implications for the study of individual movement and, also, species distributions. The central point of this paper is the idea that multiscale modelling of a dynamic system can be achieved using stochastic processes for which both the short-term transition density and the long-term stationary distribution are explicitly formulated (in particular, here, MCMC samplers). Although we have presented this method for the analysis of animal movement and resource selection, we expect that the underlying idea could have other ecological applications. For example, this problem is reminiscent of population genetics, where both the microscopic heritability laws and the macroscopic allele frequencies are of interest.
At the level of the individual, we have recognised a tendency in the current literature to embed increasingly realistic movement models in SSF analyzes. We hazard that the subtext of this trend is the intuitive notion that the habitat selection coefficients of SSF models that stay faithful to movement biology, will automatically correspond to the estimates of RSF models. As we have argued and demonstrated here, this is not necessarily the case, because SSF coefficients measure local habitat selection rather than long-term space use. Conversely, any given population distribution may be achievable by multiple movement models, just as, in the simplest of movement models, the same degree of population diffusivity can be achieved by an infinity of different movement rules, simply by trading off individual speed against path sinuosity. Although meticulous realism in movement turns out not to be a strict requirement for achieving agreement between the microscopic and macroscopic models of space use, our paper demonstrates how SSFs (through the application of statistical estimation and model selection) might in the future be used to learn about movement biology.
This manuscript serves as a proof of concept for the approach, but stops short of describing a complete workflow for the analysis of animal location data. In Local Gibbs estimation, using simulated data, we explained how the local Gibbs model can be used to estimate resource selection and movement parameters from a movement track. In a real data analysis, it would be necessary to investigate the goodness of fit. One possibility would be to simulate many locations from the fitted local Gibbs sampler, and compare the simulated and observed data in terms of some metrics of movement (e.g., distribution of step lengths). Discrepancies between features of the true and simulated data sets would point to possible model misspecifications. In addition, different models of individual movement, described by different MCMC algorithms but all guaranteed to scale up to the same long-term distribution, may be allowed to compete in a setting of statistical model selection, pointing to parsimonious explanations of the movement observations. In the estimation framework introduced in Local Gibbs estimation, likelihood-based model selection criteria, such as the Akaike information criterion (AIC), could be used to compare several candidate models. The likelihood derived from an MCMC movement model accounts for the serial correlation found in telemetry data. As such, it is a more defensible measure of likelihood than what might be obtained with other RSF approaches (Aarts et al. 2008, Fieberg et al. 2010). This modeling framework combines some of the advantages of process-based movement models and of distribution-based resource selection models. In addition to its advantages for individual-level inference, the prospect of reconciliation between RSF and SSF approaches will also benefit population-level results. In particular, the problem of formally combining the two major sources of space-use information, telemetry and transect data, has, in our experience, resisted several analytical attempts. The approach proposed here offers a solution to this problem of joint inference. For example, the steady-state distribution implied by an SSF fitted to telemetry data would be required to coincide with the utilization distribution generated by fitting a RSF to independently obtained transect data. As described in Local Gibbs estimation, the likelihood of a track (x 1 , . . ., x T ) under an MCMC movement model with transition kernel p(x t+1 |x t ) is L mov ¼ Q TÀ1 t¼1 pðx tþ1 jx t Þ and, in the same framework, the likelihood L ind of isolated survey locations {y 1 , . . ., y n } can be obtained using standard RSF methods (e.g., logistic regression or Poisson GLM). The two types of data can be combined by multiplying L mov and L ind , thus enhancing the effective sample size of the resulting estimates. Incorporating additional constraints, for example if the survey is confined to a subregion, is also straightforward.
Because it builds on the very wide and flexible class of MCMC samplers, various other movement rules could be considered. The slice sampler (Neal 2003) is an existing rejection-free sampler that shares some mathematical details with our local Gibbs sampler, and a "local" version may give some additional flexibility in movement modeling. Models of animal movement often incorporate directional persistence, such as the discrete-time and continuous-time correlated random walks (e.g., Jonsen et al. 2005, Johnson et al. 2008. Within the framework we described, this feature of movement could be modeled using non-reversible MCMC samplers, which often display this type of autocorrelation (e.g., Michel and S en ecal 2017). Such algorithms could be used for more realistic movement models.
Although we have focused on the case where the radius parameter r of the local Gibbs algorithm is taken to be constant, allowing r to be stochastic is straightforward, as mentioned above. The flexibility of the model depends in part on the choice of this distribution. More realistic features of animal movement, such as different distributions of step lengths, could thus be incorporated in the local Gibbs sampler by choosing a flexible parametric distribution for r (e.g., a gamma or Weibull distribution). A further refinement would be to combine this approach with the state-space modeling framework (Patterson et al. 2008), with the state of the process representing true location, thus incorporating measurement error on locations and giving some robustness against errors of measurement, classification, or registration in the habitat map.
The present paper therefore opens the way for future research in three vital directions: the exploration of the wealth of biological models that can be implemented with our MCMC analogues, the development of inferential methods for the integrated analysis of different data types, and the investigation into how population-level space use arises from individual rules of movement.