Differentiating the Lévy walk from a composite correlated random walk
Abstract
1. Understanding how to find targets with very limited information is a topic of interest in many disciplines. In ecology, such research has often focused on the development of two movement models: i) the Lévy walk and; ii) the composite correlated random walk and its associated arearestricted search behaviour. Although the processes underlying these models differ, they can produce similar movement patterns. Due to this similarity and because of their disparate formulation, current methods cannot reliably differentiate between these two models.
2. Here, we present a method that differentiates between the two models. It consists of likelihood functions, including one for a hidden Markov model, and associated statistical measures that assess the relative support for and absolute fit of each model.
3. Using a simulation study, we show that our method can differentiate between the two search models over a range of parameter values. Using the movement data of two polar bears (Ursus maritimus), we show that the method can be applied to complex, realworld movement paths.
4. By providing the means to differentiate between the two most prominent search models in the literature, and a framework that could be extended to include other models, we facilitate further research into the strategies animals use to find resources.
Word count: 6718
Short running head: Differentiating random search models
List of online material: Appendixes A B & C
Manuscript type: Standard article
Key words: Lévy flight, Arearestricted search, Areaconcentrated search, Animal movement, random search strategy, Hidden Markov model, Lévy foraging hypothesis
1 Introduction
Search strategies that allow targets to be found with very limited information are relevant to diverse fields of study (Bénichou et al., 2011). In particular, they have received much attention in the animal movement literature, where the two most prominent random search models are the Lévy walk and the composite correlated random walk (CCRW), with its associated arearestricted search behavior (Fauchald & Tveraa, 2003; Viswanathan et al., 2008; Dragon et al., 2012). The Lévy walk is a popular but controversial movement model that is defined as a random walk with a powerlaw distribution describing the step length frequency (Benhamou, 2007; Edwards et al., 2007; Humphries et al., 2012; Sims et al., 2012; Pyke, 2015). This distribution has a characteristic heavy tail that allows for arbitrarily long step lengths. Lévy walks are sometimes inaccurately referred to as Lévy flights in the movement literature (see Pyke, 2015). Arearestricted search (also known as areaconcentrated search) is the process whereby animals restrict their movement to the vicinity of recent captures, and is particularly useful in heterogeneous environments (Kareiva & Odell, 1987; Benhamou, 1992). Arearestricted search is one of two behaviors often modeled with CCRWs or similar composite random walks (Benhamou, 1992, 2007). Such twobehavior models typically consist of ‘extensive’ and ‘intensive’ phases, and are often used to identify foraging events or locate food patches from movement data (e.g., Jonsen et al., 2007; Dragon et al., 2012; Knell & Codling, 2012). Each behavior is related to a specific part of the landscape. The intensive search behavior is triggered by the encounter of a food item. This behavior is called arearestricted search because the animal uses low speed and large turning angles to remain within a food patch and thus increase the probability of detecting prey. The extensive search behavior is resumed after repeated unsuccessful searches. It uses fast and nearly straight movement to find the next food patch. Both the Lévy walk and CCRWs with arearestricted search have been claimed to be optimal under certain conditions (Benhamou, 1992; Viswanathan et al., 1999, but see James et al. 2011) and both have empirical support (e.g., Dragon et al., 2012; Humphries et al., 2012).
Although the processes underlying these two search models differ widely in their biological interpretation, their movement patterns are similar and difficult to differentiate. Many have argued that the CCRW could be confounded with the Lévy walk (Benhamou, 2007; Plank & James, 2008; Plank & Codling, 2009; Codling & Plank, 2011) and the disparate formulation of these models hinders their direct comparison (AugerMéthé et al., 2011). In response, new methods to identify the Lévy walk have been developed (Reynolds, 2012; Gautestad, 2012, 2013, but see AugerMéthé et al. 2014). However, these improved methods cannot be used to quantify the evidence for the CCRW. Quantifying the level of evidence for each model is important as it both reduces the potential for misidentification and allows for a more comprehensive analysis of search strategies. Recently, methods have been proposed that simultaneously fit the Lévy walk and models approximating the CCRW (Jansen et al., 2012; Plank et al., 2013). Although these methods represent significant improvements over previous approaches, they do not fully represent the CCRW as they lack turning angles and temporal correlation in behaviors. Turning angles are an essential part of movement and are crucial for distinguishing between the two movement behaviors found in the CCRW (Benhamou, 1992; Pyke, 2015). Temporal correlation in behaviors is an inherent characteristic of the CCRW because it is required to create the tortuous movement that allows the animal to remain in a food patch.
Here, we present a new method for differentiating between the movement patterns of the Lévy walk and CCRW. In the proposed method, the CCRW is represented by a hidden Markov model (HMM) that incorporates turning angles and behavioral persistence (similar to Langrock et al., 2012). For comparability, the Lévy walk and two null models are modified to incorporate turning angles. Likelihood functions for these models are created because they are essential for a set of statistical measures that assess both the relative and absolute support for each model. Using a simulation study, we show that our method can be used to successfully differentiate between the movement patterns of Lévy walks and CCRWs and to assess the relative and absolute fit of the models. We demonstrate the applicability of our method by applying it to the movement paths of two polar bears (Ursus maritimus).
2 Methods
2.1 Development of the proposed method
Our proposed method consists of likelihood functions representing each search model and statistical measures that use these likelihoods to assess the support for each model.
2.1.1 Likelihood functions
Our likelihood functions use the information from both movement measures of a step: . The step length, , is defined as the distance between the starting and ending locations of the step. The turning angle, , is defined as the angle of a step relative to the previous step direction. Only steps with a sufficient number of locations to measure both a step length (i.e., requires two locations) and a turning angle (i.e., requires three locations) are included. In addition, because we focus on the case where animals are moving and potentially searching (i.e., not performing behaviours such as resting) and because Lévy walks only model step lengths greater than 0 (see below), we exclude steps with identical start and end points. Excluding steps is possible because the models either assume that each measure of movement is independent and identically distributed or, in the case of the HMM, are built to handle missing steps. In this section, we present the development of the likelihood functions representing a CCRW, Lévy walk, and two null models. The four likelihoods differ mainly in the probability density functions (PDFs) chosen to describe the step length and turning angle frequencies.
A CCRW is a combination of two random walks, representing two behavioral modes. Similar to Plank & Codling (2009), we describe the tortuous movement of the intensive search (hereafter denoted with subscript i) with a Brownian walk (BW) and the directed movement of the extensive search (hereafter denoted with subscript e) with a correlated random walk (CRW). The BW and CRW are two common models that differ in their turning angle distribution. While an animal following a BW has no preferred turning direction, one following a CRW has a tendency to continue in the same direction as the previous step (Codling et al., 2008).
For each behavior, we define the turning angle frequency with one of two circular PDFs. To represent the intensive search as a BW, we use a circular uniform distribution, (Appendix A: Table A.1). For the extensive search, we chose the von Mises distribution. This distribution was used in recent studies comparing Lévy walks and CCRWs (Plank & Codling, 2009; Plank et al., 2013). The von Mises distribution has two parameters: , which is the location parameter and can be interpreted as the mean angle between steps; and , which is the scale parameter and can be interpreted as the size of the directional correlation. To represent the extensive search as a CRW, we set and estimate . This von Mises distribution is similar to a circular version of the Gaussian distribution centered at 0 (Forbes et al., 2011) and is represented as (Appendix A: Table A.1).
For each behavior, we model the step lengths with a slightly modified exponential distribution, (Appendix A: Table A.1). This modified exponential distribution is often used as an alternative to the Lévy walk (e.g., Codling & Plank, 2011; Edwards, 2011; Reynolds, 2012) and was used in previous attempts to compare multiphasic movement to the Lévy walk (Jansen et al., 2012). The exponential distribution defines the probability of a step length as exponentially decreasing with increasing size. The modified exponential distribution starts at the minimum step length, , rather than starting at 0. This modification is equivalent to applying the exponential distribution to the difference between the step length and the minimum step length, . The distribution often used to model Lévy walks requires the minimum step length, , to be greater than 0 (Edwards, 2011; Forbes et al., 2011; Edwards et al., 2012). As such, the data sets used in Lévy walk studies exclude step length of 0 (i.e., when the animal remains stationary). The modified exponential distribution can model data sets that exclude step lengths of 0 and thus makes our CCRW directly comparable to Lévy walk models.
Each exponential distribution has two parameters to estimate: the minimum step length, , and the rate parameter, . While the minimum step length, , is assumed to be the same for both behaviors, differs between behaviors: and . We can interpret as the inverse of the mean step length (Forbes et al., 2011), or more precisely as the inverse of the mean difference between step lengths and the minimum step length, . Thus a difference between and captures differences in the distances moved in each behavior. By combining the step length and turning angle distributions, we get the following observation PDFs associated with each behavior: {linenomath*}
(eqn 1) 
and {linenomath*}
(eqn 2) 
The observation PDFs describing the movement of each behavior are combined through what is referred as a mixing distribution. The choice of mixing distribution is an important difference between our model and the previous attempts to compare multiphasic movement to the Lévy walk (Jansen et al., 2012; Plank et al., 2013). Previous models combined the observation probabilities through an independent mixing distribution, where the probability of intensively searching and that of extensively searching are independent of previous probabilities and constant through time. Although these models provide good approximations to the movement of an animal that has two behaviors, they do not represent the temporal correlation in behaviors that a HMM can provide. Behavioral persistence is crucial when modeling the CCRW without including environmental variables as the trigger for behavioral switches. In our case, we implicitly represent the spatial correlation that a patchy landscape would create with firstorder temporal correlation in behavior. Thus, unlike models with an independent mixing distribution, the order of the observations is important in a HMM.
We used the methods of Zucchini & MacDonald (2009) to create a HMM from our observation probabilities. The mixing distribution is a firstorder Markovian process, which models the transition between the behavior of consecutive steps through the transition probability matrix: {linenomath*}
(eqn 3) 
where and are the probabilities of remaining in the intensive and extensive search behaviors, respectively, and and are the probabilities of switching from intensive to extensive and from extensive to intensive, respectively. Because the duration of each movement phase follows a geometric distribution, and can be interpreted as the mean number of steps the animal remains in the intensive and extensive search, respectively. Thus, an animal that remains on average more than two steps in the same search behavior will have and . As the probability of being in a behavior depends on the previous probabilities, we need to define the initial probability of being in each behavior: {linenomath*}
(eqn 4) 
where and are the probabilities of starting in the intensive and extensive search behaviors, respectively. The likelihood of the CCRW is: {linenomath*}
(eqn 5) 
where is a column vector of ones and is the observation probability matrix that incorporates the probability of being in each behavior as defined by eqn 1 and eqn 2: {linenomath*}
(eqn 6) 
The expanded formula of the likelihood can be found in Table 1. As mentioned above, we used the von Mises and exponential distributions in our HMM because they have been used in related research (Plank & Codling, 2009; Jansen et al., 2012; Plank et al., 2013). However, this approach can be generalized. The HMM framework is flexible, and other turning angle and step length distributions can be used to create CCRWs (e.g., wrapped Cauchy and Weibull distributions, see: Langrock et al., 2012).
To make the likelihood of the Lévy walk comparable to the CCRW, we used a PDF for the turning angle in addition to the PDF that is generally used to describe the step lengths of the Lévy walk (Table 1). Following others (e.g., Bartumeus et al., 2005; Plank et al., 2013; Pyke, 2015), we assume that the turning angle of the Lévy walk is uniform. Thus, we used the same circular uniform PDF, , as described above (Appendix A: Table A.1). Two step length PDFs can be used to describe the Lévy walk. One represents the pure Lévy walk, the other represents the truncated Lévy walk (TLW). Unlike the pure Lévy walk, the TLW places an upper bound on the size of possible step lengths, making it biologically plausible (Viswanathan et al., 2008). As a result, the TLW is often used as the Lévy walk model for animal movement (e.g., Sims et al., 2012). The step length PDF of the TLW is the truncated Pareto, (Appendix A: Table A.1). This distribution has three parameters to estimate: the shape parameter, , which increases the probability of long step length as it decreases, the minimum step length, , and the maximum step length, , which represents either the greatest step length an animal can make or the greatest distance between prey encounters. The likelihood for the TLW is: {linenomath*}
(eqn 7) 
While we focused on the TLW in the main text, we present analyses of the pure Lévy walk in Appendix A.
To verify that the complexity associated with the CCRW and TLW is required to explain the data, it is important to compare these models against simpler ones. Therefore we used likelihood functions for two simpler models: the BW and CRW (Table 1). The BW is a null model representing an individual moving randomly in space, while the CRW represents movement with directional persistence. These models are closely related to the null models used in Lévy walk studies (Bartumeus et al., 2005; Edwards et al., 2007) and use the same circular and exponential PDFs as the observation PDFs of the CCRW (eqn 1 and eqn 2). The truncated version of the exponential distribution is sometimes used as a null model in Lévy walk studies (e.g., Edwards et al., 2007) and we present analyses of the truncated version of these two models in Appendix A.
2.1.2 Statistical measures
To assess the support for each search model, we used the likelihood functions described above with a set of statistical measures. First, we computed the maximum likelihood estimates (MLEs) of the model parameters and calculated their confidence intervals through likelihood surface analyses. Second, we compared the fit of the models with Akaike Information Criterion () and Akaike weights. Finally, we tested the absolute fit of the models through analyses of pseudoresiduals. We performed these analyses with R 3.1.1 (R Core Team, 2014). The R code and Rcpp source code for the R package we have developed is available on GitHub (https://github.com/MarieAugerMethe/CCRWvsLW).
We used maximum likelihood to estimate the parameters of the models described above (Table 2). We used known analytical solutions for the MLE of and (Edwards et al., 2012). For the remaining parameters, we used numerical optimizing functions and, in the case of our CCRW, we used the ExpectationMaximization (EM) algorithm described by Zucchini & MacDonald (2009). We used the EM algorithm for our CCRW because it could be readily coded with Rcpp (Eddelbuettel & François, 2011). The resulting Rcpp algorithm was orders of magnitude faster than using R’s numerical optimizers to directly maximize the likelihood. Given that we fit our CCRW to 77 700 simulations, computational efficiency was an important consideration (see next section). A disadvantage of using the EM algorithm over the direct numerical maximization is the need to estimate (Zucchini & MacDonald, 2009), a parameter with little biological relevance. While both methods generally produce similar results, the EM algorithm is harder to code than the numerical maximization of the likelihood (MacDonald, 2014). Thus, while our fast Rcpp EM algorithm was required for our simulation study, the direct numerical maximization of the CCRW likelihood would have been easier to implement and would be an adequate solution to fit a CCRW to empirical data. Newly developed HMMs may be more easily implemented using direct numerical maximization of the likelihood (e.g., Langrock et al., 2012).
To estimate the confidence intervals of the parameters, we used the quadratic approximation described by Bolker (2008). This method uses the Hessian of the negative log likelihood at its minimum value. As the analytical solution of and is to use the minimum and maximum observed step lengths (Edwards et al., 2012) and the estimated value from the EM algorithm for depends only on the observations of the first step (Zucchini & MacDonald, 2009), it is difficult to estimate confidence intervals for these three parameters. We only provide point estimates for them.
The main goal of our likelihood functions is to identify which model fits the data best. To do so, we compared the relative fit of the models using and Akaike weights (Burnham & Anderson, 2002). The model with the lowest is considered to be the best model. To measure the weight of evidence the best model has over the other models, we calculated Akaike weights, , from the values of the models (Burnham & Anderson, 2002). Akaike weight values vary between 0 and 1, with a weight close to 1 suggesting that the data strongly support the model over the other models investigated.
As the best model according to and Akaike weights can still be a poor representation of the data, it is important to verify its absolute fit (AugerMéthé et al., 2011). In the context of Lévy walk analyses, the suggested test of absolute fit is a Gtest (Edwards et al., 2007; Edwards, 2011), a test that assumes that observations are independent of one another. This assumption is violated in the case of the CCRW because this model incorporates temporal autocorrelation. Hence, we modified the test of absolute fit by applying the Gtest to pseudoresiduals rather than to observations. We used ordinary uniform pseudoresiduals, which are residuals that account for the interdependence of observations and are uniformly distributed when the model adequately describes the data (Zucchini & MacDonald, 2009). We performed a Gtest that compares the observed frequency of these pseudoresiduals to a discretized uniform distribution. To reduce the potential bias associated with bins that have small expected values, we used William’s correction and ensured that each bin had 10 expected pseudoresiduals (Sokal & Rohlf, 1981). We applied the Gtest to the pseudoresiduals of step length and turning angle independently and subsequently combined their pvalues using Fisher’s method (Sokal & Rohlf, 1981). One can further investigate the absolute fit of the models by looking for the presence of autocorrelation in the pseudoresiduals. Appendix B describes pseudoresiduals and the test of absolute fit in more detail.
2.2 Simulation study
We used simulations to assess whether our method can differentiate between the TLW and CCRW (code also available on Github: https://github.com/MarieAugerMethe/CCRWvsLW). Because parameter values affect the resemblance of these models (AugerMéthé et al., 2014), we simulated the CCRW and TLW on a range of parameter values. For each set of parameters, we ran 50 simulations. Each simulation created a movement path of 500 biologically relevant steps (i.e., representing animal movement decisions, for which a constant time interval is not assumed, see Appendix C for simulations investigating alternative conditions). For each simulation, we used our proposed method to estimate the parameter values and calculate the Akaike weights of all models. This allowed us to verify that the method could accurately estimate parameters and appropriately differentiate between models. To assess whether the true model was rejected at the appropriate level, we also calculated the pvalue of the absolute fit test associated with the simulated model.
To simulate the CCRW, we initialized the movement path by selecting the starting behavior, either or , using a Bernoulli distribution with probability of being in the intensive search behavior defined by . If the behavior was the intensive search, we randomly selected a turning angle from a circular uniform distribution and a step length from an exponential distribution with . If the behavior was the extensive search, we randomly selected a turning angle from a von Mises distribution with and a step length from an exponential distribution with . After selecting the turning angle and step length for the first step, we selected the next behavioral state with a Bernoulli distribution that used the transition probability appropriate for the current behavior (i.e., if in intensive search and if in extensive search). As for the first step, we then selected a step length, a turning angle and the behavioral state for the next step from the appropriate distributions. This process was continued until the last step of the movement path. Our CCRW has seven parameters (Tables 1 and 2). We fixed the values of , , and , to 0.5, 0.01, and 1, respectively. We varied the value of to , that of to , that of to , and that of to . By choosing , the step lengths from the extensive search behavior were either the same length or longer on average than those from the intensive search. We chose the values of to be because the intensive search of the CCRW is efficient only if the animal remains multiple steps in a food patch. In contrast, we allowed to be because an efficient extensive search for a food patch can be produced in one step. All 720 combinations of these parameters were simulated.
For each step of the TLW simulations, we randomly selected a turning angle from a circular uniform distribution, and a step length from a truncated Pareto distribution. The TLW has three parameters (Tables 1 and 2). We set and varied the value of to (1.1, 1.2, …, 2.9) and to (100, 1000, 10000). All 57 combinations of these parameters were simulated.
2.3 Application to empirical data
To demonstrate its usefulness, we applied our method to the movement path of two polar bears from the Western Hudson Bay, Manitoba, Canada (data available on the University of Alberta Education & Research Archive: http://hdl.handle.net/10402/era.40993). These two adult females were captured in September 2010 using the standard immobilization techniques (Stirling et al., 1989, approved by the University of Alberta BioSciences Animal Policy and Welfare Committee  Protocol #6001004) and were collared with Gen IV collars from Telonics (Telonics Inc., Mesa, AZ, U.S.A). The collars were programmed to collect GPS locations at varying frequencies throughout the year. We used data from April 2011, the longest period with high frequency locations (location taken every 30 minutes) and a period where bears search for food (Pilfold et al., 2012; Thiemann et al., 2006). These two bears were on the sea ice during this period.
We applied our method to the data from each individual separately after estimating biologically relevant steps from the raw GPS data. Multiple techniques can be used to transform locations collected at regular time intervals into a timeseries of biologically relevant steps (e.g., Turchin, 1998; Codling & Plank, 2011; Humphries et al., 2013). In part for its ease of use, we used the local turn technique, which creates one step out of all consecutive sampled steps with a turning angle smaller than a threshold angle (see Codling & Plank, 2011). We have shown elsewhere that using these types of techniques can results in misidentifying CCRWs for the Lévy walk (Codling & Plank, 2011; Plank et al., 2013). However, such misidentification occurs mainly when high threshold angles are used (Codling & Plank, 2011; Plank et al., 2013). We chose a threshold angle of , meaning that any sampled step within the forward sector is interpreted as part of a biologically relevant step. Thus resulting steps are created from movement in the same general direction and the threshold is small enough that it is unlikely to result in misidentification. We applied our method to empirical data to demonstrate how to interpret results and to show the performance of our method with real animal movement paths, which, unlike simulated movement, are complicated by factors such as missing data. See Appendix C for a simulation study exploring the effects of the local turn method on movement paths and a description of how missing data are handled.
3 Results
3.1 Simulation results
The Akaike weights could differentiate the Lévy walk from a CCRW. When the CCRW was simulated, of the Akaike weight values of the CCRW exceeded 0.99 and the Akaike weight values of TLW never exceeded 0.01 (Fig. 1A). Although the CCRW simulations were never misidentified as a TLW, of the summed Akaike weight values of the null models, , exceeded 0.5. This only occurred when the step length distribution of the extensive search was relatively close to that of the intensive search, . In addition, this was generally limited to cases when the tendency to continue in the same direction was relatively low, . When the TLW was simulated, 96.7% of the Akaike weight values of the TLW exceeded 0.99 (Fig. 1B). While 3.3% of the Akaike weight value of the CCRW exceeded 0.01, only 0.1% exceeded 0.5. Note that, due to underflow, we were unable to estimate the value of the CCRW for 0.3% of the simulations. The Akaike weights results presented above and MLE results below ignore all problematic simulations.
In addition to differentiating between the two models, our method was capable of recovering the parameter values of the CCRW and TLW. As some parameter estimates can help identify whether the data are consistent with the Lévy walk or with a CCRW with an efficient arearestricted search behaviour, it is important for our method to adequately estimate their values. The CCRW requires specific values for and to be an efficient search model. The values of and can help further characterize the CCRW used by the animal. The TLW requires specific values for to be an efficient Lévy walk. For most parameters of the simulated CCRW and TLW, the median of the estimated values was close to their true value (Figs. 2 and 3). There were three exceptions. First, the estimated values of the initial probability of being in the intensive search of the CCRW, , approached either 0 or 1, not 0.5 (Fig. 2F). Second, some estimates of the minimum step length, , were positively biased, and those of maximum step length, , were negatively biased (Figs. 2G and 3BC). Third, similar to the Akaike weights, the estimates of most parameters of the CCRW were less accurate when the movement patterns of the two behaviors were similar. Specifically, the estimates were less reliable when the simulations values of were relatively close to . The estimated values of most parameters were much closer to the true value when simulations with were excluded. While the point estimates were generally reliable, the confidence intervals, as estimated with the quadratic approximation, had a tendency to be too narrow and excluded the correct simulation value more that of the time (Appendix C Table C.1) and often overlapped with the parameter space boundary (CCRW: , TLW: ), we thus used the more computationally intensive profile likelihood for the empirical results.
Finally, our tests of absolute fit had rejection rates adequate for the selected level of 0.05 (pvalue ). The proportion of simulated CCRWs that were rejected from being CCRW was 0.062. Similarly, the proportion of simulated TLWs that were rejected from being TLW was 0.067.
Appendix C shows that our method is affected by the local turn method, a technique used to transform raw GPS data into biologically relevant steps. Because the local turn method amalgamates all consecutive sampled steps with less than a threshold angle, and thus removes small turning angles, parameter estimates were heavily affected. In particular, estimates were negatively biased. Using the local turn method also affected the test of absolute fit. When movement paths are transformed with such method, as in the case of our polar bear data, the test of absolute fit should only use the pseudoresiduals associated with the step lengths. The results in Appendix C demonstrate that the local turn method strongly affects the parameter estimates and the test of absolute fit. However, when used with a small threshold angle, this technique did not decrease appreciably the capacity of our method to distinguish between the TLW and a CCRW.
3.2 Empirical results
The best model for the two empirical movement paths was our CCRW (Table 3). For Bear 2, the Akaike weights indicated that the CCRW was a much better model than the other alternatives. However, the Akaike weight of the CCRW for Bear 1 was only 0.55, with some evidence that the CRW may have been a more parsimonious description of the movement data (Table 3 and Fig. 4). While the best model was the CCRW, both movement paths were significantly different from it (Table 3). The movement path of Bear 1 was also significantly different from the CRW (). A visual representation of the fit of the models is presented in Fig. 4.
To identify whether the movement paths were consistent with the best model, we investigated the parameter estimates of the CCRW. For Bear 2, all parameters were consistent: and (Table 2). In contrast, not all parameters for Bear 1 were consistent with the CCRW. While as expected, .
4 Discussion
Through the analysis of TLW and CCRW simulations, we have demonstrated that our method can differentiate between the movement patterns of a Lévy walk and CCRW. The Akaike weights identified the correct search model, except for a few instances. The Akaike weights also distinguished the TLW and CCRW from our two null models. The rare exceptions occurred when both the intensive and extensive search behaviors of the CCRW simulations had similar step length and turning angle distributions. This was expected. Other methods developed to distinguish the intensive from the extensive search are also less efficient when the movement of these behaviors are similar (Knell & Codling, 2012). When the two behaviors are similar, models describing them as one behavior can be sufficient. The ability of our method to differentiate between a CCRW and null models would likely increase with sample size.
The simulation analyses also indicated that most parameter estimates of the TLW and CCRW were reliable. The estimates of the important parameters of both models (e.i., , , , , and ) were generally reliable and accurate. These are the only parameters that should be used to help identify whether the empirical data support the Lévy walk or the CCRW. No biological interpretation should be based on the probability of starting in the intensive search behavior, . As described by Zucchini & MacDonald (2009), the estimates from the EM algorithm for this parameter approached either 0 or 1 as will be one of the two unit vectors. Caution should be taken when interpreting the minimum, , and maximum, , step lengths. Even though using the minimum and maximum observed step lengths are the MLEs, and is the suggested method to estimate these values for TLW (Edwards et al., 2012), some of their estimates were biased. One likely explanation, is that 500 steps was too small a sample to accurately estimate these parameters. The estimates of most parameters of the CCRW suffered when the two search behaviors were not substantially different.
Because precise methods, such as the likelihood profile, become highly unpractical and computationally demanding when models have more than two or three parameters to be estimated, Bolker (2008) recommends the use of the quadratic approximation for estimating confidence intervals. Because the CCRW has seven parameters to be estimated, we investigated whether such approximation could be used. The simulation study showed that these approximated confidence intervals were often too narrow and excluded the simulation value. The quadratic approximation can be inaccurate when the parameter estimated is at the boundary of its parameter space (Zucchini & MacDonald, 2009). This approximation is symmetric around the MLE, thus might exceed the boundary of parameter space. This occurred for many simulations. For the polar bear data we estimated the confidence intervals using the likelihood profile.
The simulation results showed that our test of absolute fit was adequate, albeit with observed rejection rates that were marginally greater than the expected rate of 0.05. Thus our test had a slightly higher level of type I error than specified by the level. This problem could be associated with the known negative bias in pvalues of Gtests when sample size and expected values are small (Sokal & Rohlf, 1981). We have also explored the use of a number of other tests, such as tests of normality on normal pseudoresiduals (see Zucchini & MacDonald, 2009, for description of normal pseudoresiduals). None have outperformed the one presented here.
Some sampling procedures, in particular subsampling and the definition of steps by the significant turns, can cause Akaike weights to select Lévy walk models when CCRWs are simulated (Plank & Codling, 2009; Codling & Plank, 2011; Plank et al., 2013). Although our method is likely to be affected by such procedures, it has features that are known to decrease misidentification errors. In particular, it was shown that including an approximation of the CCRW and tests for the absolute fit mitigates the risks of such errors (Plank et al., 2013). Indeed, through a simulation study, we showed that, while parameter estimates and test of absolute fit were affected by the local turn method with a threshold angle of , our method’s capacity to distinguish between the movement patterns of CCRWs and TLWs remained almost unaffected. We have not fully explored the effects of data sampling and handling on the accuracy of our method. Future work should investigate how sampling procedures impact the capacity of our method to differentiate between the two models.
Overall, our simulation study showed that we can differentiate the Lévy walk from a strong alternative, such as a CCRW. The Akaike weights could differentiate the Lévy walk from a CCRW that used a combination of exponential distributions, something that is difficult to accomplish with other methods (Benhamou, 2007; Plank et al., 2013; AugerMéthé et al., 2014). Other alternative models, including other formulations of the CCRW, could result in movement patterns similar those of a Lévy walk. In many cases, it might be important to compare the Lévy walk to a wider range of alternative models. Because Akaike weights are relative measures of fit, it is important to verify that the best model describes the data adequately (AugerMéthé et al., 2011). Our simulation study demonstrates that our test of absolute fit can identify whether the model describes the data adequately. Finally, our simulation study shows that most parameter estimates are reliable and thus can be used to further investigate whether the data is consistent with the best model. Thus, the simulation study suggests that we could infer support for a model when: 1) it is compared to adequate alternatives, 2) has much higher Akaike weight values than the other models analysed, 3) sufficiently describes the data according to a test of absolute fit, and 4) has parameter estimates consistent with the hypothesis it represents.
We demonstrated how to interpret the results of our method by applying it to empirical data. Our results suggested that the two bears differed in their movement patterns. For Bear 2, the Akaike weights and parameter estimates suggested that the movement path was better represented by the CCRW. For Bear 1, the Akaike weights suggested that although the CCRW was the best model, the CRW, a onebehaviour null model, might be sufficient to explain the data. These two bears differed in their reproductive status: Bear 1 was accompanied by a yearling at capture while Bear 2 was accompanied by a cuboftheyear. Females with cubsoftheyear move smaller distances, avoid adult males to reduce the risk of infanticide, and use lower quality habitat than other bears (Stirling et al., 1993; Amstrup et al., 2000; Pilfold et al., 2014). Thus, it is possible that females with cubsoftheyear used different search strategies than other females and this difference could have resulted in the difference observed between the two bears.
An additional explanation for the difference between these two bears is that the quality of their movement data differ (Fig. 4). The results for Bear 1 demonstrated that our method can handle large amount of missing data. However, as with most analytical methods, missing data can impact biological interpretation. Specifically, reduced sample size likely hinders our method’s ability to differentiate between models and between the two behaviours of the CCRW. In addition, missing locations divides the path into smaller steps, which has the potential to impact model fit. Thus, we advise caution when interpreting results for movement paths with many missing locations.
The movement path of each bear was significantly different from the best model. This indicates that while better than the other alternatives, the best model is not sufficient to explain polar bear movement. One could easily extend the set of models explored by investigating multiple versions of the HMM. While our choices were made to reduce the number of parameters to be estimated or to ensure that certain characteristics of the CCRW were respected, there are many ways in which the CCRW could be modeled and certain changes could increase its absolute fit. Our CCRW used specific distributions for the frequency of step lengths and turning angles. The choice of such distributions can affect the movement behavior of random walks (Codling et al., 2010) and other distributions have been used in some multiphasic movement models (e.g., wrapped Cauchy and Weibull distributions, see: Morales et al., 2004; Langrock et al., 2012). In addition, by using a simple HMM, we are assuming that the number of steps an animal makes in each behavioral phase follows a geometric distribution. However, the autocorrelation in pseudoresiduals indicates that this assumption might be violated and that a firstorder Markov process might be an inaccurate representation of the switching probabilities for polar bears (see Appendix B). One could relax this assumption by using a hidden semiMarkov model (Langrock et al., 2012). While we explored only one version of a CCRW, our framework allow empiricists to explore a variety of models by simply altering the characteristics of the HMM.
While exploring a larger variety of CCRWs is likely to increase the absolute fit of our model, it is unlikely that we have sampled the movement paths at the exact scale at which the animals are making their decisions. Sampling scale affects behavioural inference made from movement data (e.g., Codling & Hill, 2005; Andersen et al., 2008; Plank & Codling, 2009). Thus, a lack of strong evidence for the Lévy walk and CCRW at the scale at which we have sampled our movement paths does not preclude the possibility for such evidence at different scales. Investigating the evidence for these movement models across multiple scales may be useful (e.g., Fryxell et al., 2008; Gautestad, 2013; Seuront & Stanley, 2014). Finally, it is possible that we are missing important characteristics of polar bear movement. For example, some polar bears move against sea ice drift and ignoring drift can impact interpretation of movement paths (Mauritzen et al., 2003; Gaspar et al., 2006). Thus, an important extension for polar bears might be the inclusion of drift in the analysis (e.g., Gaspar et al., 2006; Girard et al., 2006).
4.1 Conclusion
We have developed likelihood functions for models representing the Lévy walk and CCRW that make it possible to directly compare the evidence for these two prominent hypotheses. Unlike recently developed methods, our method uses information from both step lengths and turning angles, and incorporates the temporal autocorrelation inherent in the CCRW. Our simulation study showed that our method could differentiate between the two models. By applying our method to the movement path of two polar bears, we showed that our method can give easily interpretable results and handle complex movement paths. The specific model that we used for the CCRW is just one of the many CCRWs that could be created using a HMM. For example, alternate step length and turning angle distributions, such as Weibull and wrapped Cauchy distributions, could be used to create other multibehavior models with different characteristics (e.g., Langrock et al., 2012; Morales et al., 2004). We hope that application of this method to empirical data will further our understanding of the mechanisms used by animals to find resources.
5 Acknowledgments
We thank Devin Lyons, Craig DeMars, Roland Langrock, and anonymous reviewers for improving the manuscript, Ulrike Schlägel, and Alexei Potapov for technical support, and Dennis Andriashek, Nicholas Lunn, and Evan Richardson for fieldwork assistance. We received funding or scholarships from Alberta InnovatesTechnology Futures, Aquarium du Québec, ArcticNet, Canada Research Chairs program, Canadian Association of Zoos and Aquariums, Canadian Circumpolar Institute, Canadian Wildlife Federation, Environment Canada, Killam trust, Natural Sciences and Engineering Research Council of Canada, Northern Scientific Training Program, Polar Continental Shelf Project, Polar Bears International, Quark Expeditions, Steve and Elaine Antoniuk, University of Alberta, and World Wildlife Fund (Canada).
References
 Amstrup et al. (2000) Amstrup, S.C., Durner, G.M., Stirling, I., Lunn, N.J. & Messier, F. (2000) Movements and distribution of polar bears in the Beaufort Sea. Canadian Journal of Zoology, 78, 948–966.
 Andersen et al. (2008) Andersen, M., Derocher, A.E., Wiig, Ø. & Aars, J. (2008) Movements of two Svalbard polar bears recorded using geographical positioning system satellite transmitters. Polar Biology, 31, 905–911.
 AugerMéthé et al. (2014) AugerMéthé, M., Plank, M.J. & Codling, E.A. (2014) Distinguishing between Lévy walks and strong alternative models: Comment. Ecology, 95, 1104–1109.
 AugerMéthé et al. (2011) AugerMéthé, M., St. Clair, C.C., Lewis, M.A. & Derocher, A.E. (2011) Sampling rate and misidentification of Lévy and nonLévy movement paths: comment. Ecology, 92, 1699–1701.
 Bartumeus et al. (2005) Bartumeus, F., Luz, M.G.E.D., Viswanathan, G.M. & Catalan, J. (2005) Animal search strategies: a quantitative randomwalk analysis. Ecology, 86, 3078–3087.
 Benhamou (1992) Benhamou, S. (1992) Efficiency of areaconcentrated searching behaviour in a continuous patchy environment. Journal of Theoretical Biology, 159, 67–81.
 Benhamou (2007) Benhamou, S. (2007) How many animals really do the Lévy walk? Ecology, 88, 1962–1969.
 Bénichou et al. (2011) Bénichou, O., Loverdo, C., Moreau, M. & Voituriez, R. (2011) Intermittent search strategies. Reviews of Modern Physics, 83, 81–129.
 Bolker (2008) Bolker, B.M. (2008) Ecological models and data in R. Princeton University Press, Princeton.
 Burnham & Anderson (2002) Burnham, K.P. & Anderson, D.R. (2002) Model selection and multimodel inference: a practical informationtheoretic approach. 2nd edn. Springer, New York.
 Codling & Hill (2005) Codling, E. & Hill, N. (2005) Sampling rate effects on measurements of correlated and biased random walks. Journal of Theoretical Biology, 233, 573–588.
 Codling et al. (2010) Codling, E.A., Bearon, R.N. & Thorn, G.J. (2010) Diffusion about the mean drift location in a biased random walk. Ecology, 91, 3106–3113.
 Codling & Plank (2011) Codling, E.A. & Plank, M.J. (2011) Turn designation, sampling rate and the misidentification of power laws in movement path data using maximum likelihood estimates. Theoretical Ecology, 4, 397–406.
 Codling et al. (2008) Codling, E.A., Plank, M.J. & Benhamou, S. (2008) Random walk models in biology. Journal of the Royal Society Interface, 5, 813–834.
 Dragon et al. (2012) Dragon, A.C., BarHen, A., Monestiez, P. & Guinet, C. (2012) Comparative analysis of methods for inferring successful foraging areas from Argos and GPS tracking data. Marine Ecology Progress Series, 452, 253–267.
 Eddelbuettel & François (2011) Eddelbuettel, D. & François, R. (2011) Rcpp: Seamless R and C++ Integration. Journal of Statistical Software, 40, 1–18.

Edwards (2011)
Edwards, A.M. (2011) Overturning conclusions of Lévy flight movement
patterns by fishing boats and foraging animals.
Ecology, 92, 1247–1257.
 Edwards et al. (2012) Edwards, A.M., Freeman, M.P., Breed, G.A. & Jonsen, I.D. (2012) Incorrect likelihood methods were used to infer scaling laws of marine predator search behaviour. PLoS ONE, 7, e45174.
 Edwards et al. (2007) Edwards, A.M., Phillips, R.A., Watkins, N.W., Freeman, M.P., Murphy, E.J., Afanasyev, V., Buldyrev, S.V., da Luz, M.G.E., Raposo, E.P. & Stanley, H.E. (2007) Revisiting Lévy flight search patterns of wandering albatrosses, bumblebees and deer. Nature, 449, 1044–1048.
 Fauchald & Tveraa (2003) Fauchald, P. & Tveraa, T. (2003) Using firstpassage time in the analysis of arearestricted search and habitat selection. Ecology, 84, 282–288.
 Forbes et al. (2011) Forbes, C., Evans, M., Hastings, N. & Peacock, B. (2011) Statistical distributions. 4th edn. WileyInterscience, Hoboken.
 Fryxell et al. (2008) Fryxell, J.M., Hazell, M., Börger, L., Dalziel, B.D., Haydon, D.T., Morales, J.M., McIntosh, T. & Rosatte, R.C. (2008) Multiple movement modes by large herbivores at multiple spatiotemporal scales. Proceedings of the National Academy of Sciences, 105, 19114–19119.
 Gaspar et al. (2006) Gaspar, P., Georges, J.Y., Fossette, S., Lenoble, A., Ferraroli, S. & Maho, Y.L. (2006) Marine animal behaviour: neglecting ocean currents can lead us up the wrong track. Proceedings of the Royal Society B, 273, 2697–2702.
 Gautestad (2012) Gautestad, A.O. (2012) Brownian motion or Lévy walk? Stepping towards an extended statistical mechanics for animal locomotion. Journal of the Royal Society Interface, 9, 2332–2340.
 Gautestad (2013) Gautestad, A.O. (2013) Animal space use: distinguishing a twolevel superposition of scalespecific walks from scalefree Lévy walk. Oikos, 122, 612–620.
 Girard et al. (2006) Girard, C., Sudre, J., Benhamou, S., Roos, D. & Luschi, P. (2006) Homing in green turtles Chelonia mydas: oceanic currents act as a constraint rather than as an information source. Marine Ecology Progress Series, 322, 281–289.
 Humphries et al. (2012) Humphries, N.E., Weimerskirch, H., Queiroz, N., Southall, E.J. & Sims, D.W. (2012) Foraging success of biological Lévy flights recorded in situ. Proceedings of the National Academy of Sciences, 109, 7169–7174.
 Humphries et al. (2013) Humphries, N.E., Weimerskirch, H. & Sims, D.W. (2013) A new approach for objective identification of turns and steps in organism movement data relevant to random walk modelling. Methods in Ecology and Evolution, 4, 930–938.
 James et al. (2011) James, A., Plank, M.J. & Edwards, A.M. (2011) Assessing Lévy walks as models of animal foraging. Journal of the Royal Society Interface, 8, 1233–1247.
 Jansen et al. (2012) Jansen, V.A.A., Mashanova, A. & Petrovskii, S. (2012) Comment on âLévy walks evolve through interaction between movement and environmental complexityâ. Science, 335, 918–918.
 Jonsen et al. (2007) Jonsen, I.D., Myers, R.A. & James, M.C. (2007) Identifying leatherback turtle foraging behaviour from satellite telemetry using a switching statespace model. Marine Ecology Progress Series, 337, 255–264.
 Kareiva & Odell (1987) Kareiva, P. & Odell, G. (1987) Swarms of predators exhibit “preytaxis” if individual predators use arearestricted search. American Naturalist, 130, 233–270.
 Knell & Codling (2012) Knell, A.S. & Codling, E.A. (2012) Classifying arearestricted search (ARS) using a partial sum approach. Theoretical Ecology, 5, 325–339.

Langrock et al. (2012)
Langrock, R., King, R., Matthiopoulos, J., Thomas, L., Fortin, D. & Morales,
J.M. (2012) Flexible and practical modeling of animal telemetry data: hidden
markov models and extensions.
Ecology, 93, 2336–2342.
 MacDonald (2014) MacDonald, I.L. (2014) Numerical maximisation of likelihood: a neglected alternative to EM? International Statistical Review, 82, 296–308.
 Mauritzen et al. (2003) Mauritzen, M., Derocher, A.E., Pavlova, O. & Wiig, Ø. (2003) Female polar bears, Ursus maritimus, on the Barents Sea drift ice: walking the treadmill. Animal Behaviour, 66, 107–113.
 Morales et al. (2004) Morales, J.M., Haydon, D.T., Frair, J., Holsinger, K.E. & Fryxell, J.M. (2004) Extracting more out of relocation data: building movement models as mixtures of random walks. Ecology, 85, 2436–2445.
 Pilfold et al. (2014) Pilfold, N.W., Derocher, A.E. & Richardson, E.S. (2014) Influence of intraspecific competition on the distribution of a wideranging, nonterritorial carnivore. Global Ecology and Biogeography, 23, 425–435.
 Pilfold et al. (2012) Pilfold, N.W., Derocher, A.E., Stirling, I., Richardson, E. & Andriashek, D. (2012) Age and sex composition of seals killed by polar bears in the eastern beaufort sea. PloS one, 7, e41429.
 Plank et al. (2013) Plank, M.J., AugerMéthé, M. & Codling, E.A. (2013) Lévy or not? Analysing positional data from animal movement paths. Dispersal, individual movement and spatial ecology: a mathematical perspective (eds. M.A. Lewis, P.K. Maini & S.V. Petrovskii), Lecture Notes in Mathematics, vol. 2071. SpringerVerlag, Berlin.
 Plank & Codling (2009) Plank, M.J. & Codling, E.A. (2009) Sampling rate and misidentification of Lévy and nonLévy movement paths. Ecology, 90, 3546–3553.
 Plank & James (2008) Plank, M.J. & James, A. (2008) Optimal foraging: Lévy pattern or process? Journal of the Royal Society Interface, 5, 1077–1086.
 Pyke (2015) Pyke, G.H. (2015) Understanding movements of organisms: it’s time to abandon the Lévy foraging hypothesis. Methods in Ecology and Evolution, 6, 1–16.
 R Core Team (2014) R Core Team (2014) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
 Reynolds (2012) Reynolds, A. (2012) Distinguishing between Lévy walks and strong alternative models. Ecology, 93, 1228–1233.
 Seuront & Stanley (2014) Seuront, L. & Stanley, H.E. (2014) Anomalous diffusion and multifractality enhance mating encounters in the ocean. Proceedings of the National Academy of Sciences, 111, 2206–2211.
 Sims et al. (2012) Sims, D.W., Humphries, N.E., Bradford, R.W. & Bruce, B.D. (2012) Lévy flight and Brownian search patterns of a freeranging predator reflect different prey field characteristics. Journal of Animal Ecology, 81, 432–442.
 Sokal & Rohlf (1981) Sokal, R.R. & Rohlf, F.J. (1981) Biometry. 2nd edn. W.H. Freeman, New York.
 Stirling et al. (1993) Stirling, I., Andriashek, D. & Calvert, W. (1993) Habitat preferences of polar bears in the western Canadian Arctic in late winter and spring. Polar Record, 29, 13–24.
 Stirling et al. (1989) Stirling, I., Spencer, C. & Andriashek, D. (1989) Immobilization of polar bears (Ursus maritimus) with Telazol^{®} in the Canadian Arctic. Journal of Wildlife Diseases, 25, 159.
 Thiemann et al. (2006) Thiemann, G.W., Iverson, S.J. & Stirling, I. (2006) Seasonal, sexual and anatomical variability in the adipose tissue of polar bears (Ursus maritimus). Journal of Zoology, 269, 65–76.

Turchin (1998)
Turchin, P. (1998) Quantitative analysis of movement: measuring and
modeling population redistribution in animals and plants.
Sinauer Associates, Inc., Sunderland.
 Viswanathan et al. (1999) Viswanathan, G.M., Buldyrev, S.V., Havlin, S., Luz, M.G.E.D., Raposo, E.P. & Stanley, H.E. (1999) Optimizing the success of random searches. Nature, 401, 911–914.
 Viswanathan et al. (2008) Viswanathan, G.M., Raposo, E.P. & Luz, M.G.E.D. (2008) Lévy flights and superdiffusion in the context of biological encounters and random searches. Physics of Life Reviews, 5, 133–150.
 Zucchini & MacDonald (2009) Zucchini, W. & MacDonald, I.L. (2009) Hidden Markov models for time series: an introduction using R, Monographs on Statistics and Applied Probability, vol. 110. Chapman and Hall/CRC, Boca Raton.
Model  Likelihood function  

CCRW  7  
TLW  3  
BW  3  
CRW  4  

[b]

Description  Bear 1  Bear 2  


Minimum step length of all four models  21  2  

Maximum step length of the TLW  12614  11789  
Probability of starting in the CCRW’s intensive search  0  0  
Probability of remaining in the CCRW’s intensive search 



Probability of remaining in the CCRW’s extensive search 



Size of the directional correlation of the CRW 



Size of the directional correlation of the CCRW’s extensive search 




Rate parameter of the exponential distribution of the BW and CRW 




Rate parameter of the CCRW’s intensive search 




Rate parameter of the CCRW’s extensive search 



Scale parameter of the truncated Pareto distribution of the TLW 


Individual  n  Akaike weight  pvalue  

CCRW  TLW  BW  CRW  CCRW  TLW  BW  CRW  Best model  
Bear 1  235  0  302.4  170.3  0.4  0.55  0.45  
Bear 2  887  0  1479.7  648.6  137.2  1.00 