Approximate Bayesian inference for a "steps and turns" continuous-time random walk observed at regular time intervals

TL;DR

This paper develops an approximate Bayesian inference method for a continuous-time animal movement model observed at regular intervals, highlighting the importance of observation frequency relative to movement decision times.

Contribution

It introduces a state-space model combining continuous-time movement decisions with discrete-time observations and applies ABC methods for inference, addressing temporal scale mismatches.

Findings

Accurate estimates when observations are less than 5 times the average change interval.

The model effectively links observation and decision scales in animal movement.

High-frequency data may not always be necessary for reliable movement estimates.

Abstract

The study of animal movement is challenging because it is a process modulated by many factors acting at different spatial and temporal scales. Several models have been proposed which differ primarily in the temporal conceptualization, namely continuous and discrete time formulations. Naturally, animal movement occurs in continuous time but we tend to observe it at fixed time intervals. To account for the temporal mismatch between observations and movement decisions, we used a state-space model where movement decisions (steps and turns) are made in continuous time. The movement process is then observed at regular time intervals. As the likelihood function of this state-space model turned out to be complex to calculate yet simulating data is straightforward, we conduct inference using a few variations of Approximate Bayesian Computation (ABC). We explore the applicability of these methods as a function of the discrepancy between the temporal scale of the observations and that of the movement process in a simulation study. We demonstrate the application of this model to a real trajectory of a sheep that was reconstructed in high resolution using information from magnetometer and GPS devices. Our results suggest that accurate estimates can be obtained when the observations are less than 5 times the average time between changes in movement direction. The state-space model used here allowed us to connect the scales of the observations and movement decisions in an intuitive and easy to interpret way. Our findings underscore the idea that the time scale at which animal movement decisions are made needs to be considered when designing data collection protocols, and that sometimes high-frequency data may not be necessary to have good estimates of certain movement processes.

Figures

Tables

Table 1. Table 1: Values for the turning angles and duration of steps for the example of Figure 1

Duration of steps ($t_i^{\prime }s$)	0.2	0.2	0.7	0.4	0.4	0.8
Turning angle (${\omega }_i^{\prime }s$)	0.32	5.65	5.81	0.02	0.11	5.81

Table 2. Table 2: Summary statistics selected to fit the model using the three ABC algorithms

Summary Statistic	Formula
(1) Point estimate for $\lambda$	$1/\sum _{j=0}^{N_{obs}}{t}_{obs,j}$, with $t_{obs,j}=\sqrt{{\left(o_{j+1,1}-o_{j,1}\right)}^2+{\left(o_{j+1,2}-o_{j,2}\right)}^2}$
(2) Point estimate for $\kappa$	$A^{-1}\left(\frac{\sum _{j=1}^{N_{obs}}\cos \left({\omega }_{obs,j}\right)}{N_{obs}}\right)$ Where ${\omega }_{obs,j}=\arctan \left(\frac{o_{j+1,1}-o_{j,1}}{o_{j+1,2}-o_{j,2}}\right)$ and $A\left(x\right)=\frac{I_1\left(x\right)}{I_0\left(x\right)}$ (Kurt Hornik and Bettina Grün, 2014)
(3) Standard deviation of the turning angle	$\sqrt{\frac{\sum _{j=0}^{N_{obs}}{\omega }_{obs,j}-\overline{{\omega }_{obs}}}{N_{obs}-1}}$
(4)Standard deviation of the step length	$\sqrt{\frac{\sum _{j=0}^{N_{obs}}{t}_{obs,j}-\overline{t_{obs}}}{N_{obs}-1}}$

Equations

$$\begin{array}[]{cc}\mu_{i,1}=&\mu_{i-1,1}+cos(\phi_{i-1})t_{i-1}\\ \mu_{i,2}=&\mu_{i-1,2}+sin(\phi_{i-1})t_{i-1}\\ \phi_{i}=&\sum_{k=1}^{i}\omega_{i}\end{array}$$

$$\begin{array}[]{cc}o_{j,1}=&\mu_{N_{j},1}+cos(\phi_{N_{j}})\left(j\Delta t-\displaystyle\sum_{k<N_{j}-1}t_{k}\right)=h_{1}(M_{0:N_{j}})\\ o_{j,2}=&\mu_{N_{j},2}+sin(\phi_{N_{j}})\left(j\Delta t-\displaystyle\sum_{k<N_{j}-1}t_{k}\right)=h_{2}(M_{0:N_{j}})\end{array}$$

$$\begin{array}{rl}L(\kappa ,\lambda ;M,O)& =P\left(O_0=o_0,O_1=o_1,\cdots ,O_{N_{obs}}=o_{N_{obs}}\right)\\ & =P\left(h(M_{0:N_1})=o_1,h(M_{0:N_2})=o_2,...,h(M_{0:N_{N_{obs}}})=o_{N_{obs}}\right)\end{array}$$

$$\sqrt{\frac{{\sum }_i(\widetilde{{\theta }_i}-{\theta }_i{)}^2}{N_{rep}}}\theta =(\lambda ,\kappa );i=1,...,N_{rep}$$

$$MD=\left[\sum _i\frac{|\widetilde{{\theta }_i}-{\theta }_i|}{{\theta }_i}\right]/N_{rep}\theta =(\lambda ,\kappa );i=1,...,N_{rep}$$

$$\begin{array}{rl}L(\kappa ,\lambda ,M,O)& =P\left(O_0=o_0,O_1=o_1,\cdots ,O_{N_{obs}}=o_{N_{obs}}\right)\\ & =P\left(h(M_{0:N_1})=o_1,h(M_{0:N_2})=o_2,...,h(M_{0:N_{N_{obs}}})=o_{N_{obs}}\right)\end{array}$$

$${\mu }_{i,1}={\mu }_{i-1,1}+cos({\varphi }_{i-1})t_{i-1}$$

$$f_V=\left(f_{\varphi }(-acos(v))+f_{\varphi }(acos(v))\right)\frac{1}{\sqrt{1-v^2}}{I}_{-1\le v\le 1}(v)$$

$$f_Z(z_1)=\int f_V(\frac{z_1}{z_2})f_t(z_2)I_{\{-z_2\le z_1\le z_2\}}(z_1)I_{\{z_2>0\}}(z_2)\cdot d{z}_2$$

$$\begin{array}{rl}{o}_{j,1}& =h({\mu }_{0:N_j,1})\\ & ={\mu }_{N_j,1}+\cos ({\varphi }_{N_j})\left(j\Delta t-\sum _{k<N_j-1}{t}_k\right)\\ & ={\mu }_{N_j,1}+V_{N_j}(c_j-W_{N_j-1})\end{array}$$

$$o_{j,1}={\mu }_{N_j,1}+(V_{N_j})({\tilde{W}}_{N_j-1})={\mu }_{N_j,1}+S_{N_j-1}$$

$$f_{S_{N_j-1}}(s_1)=\int f_{V_{N_j}}(\frac{s_1}{s_2})f_{{\tilde{W}}_{N_j-1}}(s_2)I_{\{-s_2\le s_1\le w_2\}}(s_1)I_{\{s_2\le -c_j\}}(s_2)\cdot d{s}_2$$

Full text

To account for the temporal mismatch between observations and movement decisions, we used a state-space model where movement decisions (steps and turns) are made in continuous time. The movement process is then observed at regular time intervals.

As the likelihood function of this state-space model turned out to be complex to calculate yet simulating data is straightforward, we conduct inference using a few variations of Approximate Bayesian Computation (ABC). We explore the applicability of these methods as a function of the discrepancy between the temporal scale of the observations and that of the movement process in a simulation study. We demonstrate the application of this model to a real trajectory of a sheep that was reconstructed in high resolution using information from magnetometer and GPS devices.

Our results suggest that accurate estimates can be obtained when the observations are less than 5 times the average time between changes in movement direction.

The state-space model used here allowed us to connect the scales of the observations and movement decisions in an intuitive and easy to interpret way. Our findings underscore the idea that the time scale at which animal movement decisions are made needs to be considered when designing data collection protocols, and that sometimes high-frequency data may not be necessary to have good estimates of certain movement processes.

Introduction

The way in which animals move is of fundamental importance in ecology and evolution. It plays important roles in the fitness of individuals and in gene exchange (Nathan et al., 2008), in the structuring of populations and communities (Turchin, 1998, Matthiopoulos et al., 2015, Morales et al., 2010), and in the spread of diseases (Fèvre et al., 2006). The study of movement is challenging because it is a process modulated by many factors acting at different spatial and temporal scales (Gurarie and Ovaskainen, 2011, Mevin B. Hooten et al., 2017). In general, the process of animal movement occurs in continuous time but we observe individual locations at time intervals dictated by logistic constrains such as battery life. It is necessary to be aware of this fact in order to avoid drawing conclusions about the movement process that depend on the time scale in which the observations were taken.

In this context, state-space models provide a convenient tool for movement data analysis (Patterson et al., 2008). The main idea is to estimate the latent movement process given the observational process. Thus, they consist of two stochastic models: a latent model and an observation one. The first describes the state of the animal (it could be the location, behaviour, etc) and the second one describes the observation of the state, possibly with some measurement error. Several state-space models have been proposed to model animal movement differing primarily in the temporal conceptualization of the movement process, namely discrete and continuous time formulations (McClintock et al., 2014). On one hand, discrete-time models describe movement as a series of steps and turns (or movement directions) that are performed at regular occasions (Morales et al., 2004, Jonsen et al., 2005, McClintock et al., 2012). Typically in these models the temporal scales of both the latent and the observation process are the same. Thus, the observation times coincide with the times in which the animal are assumed to make movement decisions. The advantage of this approach is that it allows the dynamics involved in the movement process to be conceptualized in a simple and intuitive way, which facilitates implementation and interpretation. On the other hand, continuous-time models have been proposed (Blackwell, 1999, Johnson et al., 2008, Harris and Blackwell, 2013) in which the movement process is defined for any time and expressed through stochastic differential equations that account for the dependence between successive locations. The observation process is then truly independent from the movement process and does not need to be recorded at regular time intervals. The continuous-time approach has the advantage of being more realistic and that the inference is not affected by the choice of scale of observation. The main drawback is probably in the interpretation of instantaneous movement parameters (e.g., those related to Ornstein-Uhlenbeck processes and other diffusion models). In such manner both approaches have advantages and disadvantages, the discrete-time models are more intuitive and easy to interpret but can be considered less realistic than continuous ones (McClintock et al., 2014).

In this work, we present a state-space model that formulates the movement process in continuous time and the observation process in discrete time (regular intervals). As a compromise between the ease of interpretation of models based on steps and turns and the realism of continuous-time models, we use a random walk where the movement decisions (steps and turns) are made in continuous time. The movement process is then observed at regular time intervals. In this model there are two different time scales: one for the latent process and one for the observation process. The advantage here is that this model allows us to differentiate between the times in which the animals make movement decisions and the times in which the observations are made. One challenge faced with the proposed model formulation is that the resulting likelihood function seems to be computationally intractable. However, simulating the movement and observation process is straightforward suggesting that likelihood-free methods such as Approximate Bayesian Computation (ABC) could be useful (Beaumont, 2010, Csilléry et al., 2010).

Here we describe, formalize, and expose the possible complications of a state-space movement model with two different temporal scales. We use stochastic simulations to evaluate the ability of three ABC techniques to recover the parameter values driving the movement process. Keeping in mind the ecological purpose behind implementing such a model, we assess the quality of these estimations with regard to the relationship between the two temporal scales. Finally, we apply the model to a high resolution trajectory of sheep to evaluate the performance of the ABC inference with real data.

Methods

Movement Model with Random Time between movement decisions

The general structure of the model is based on a correlated random walk. An individual moves in a certain direction for a certain period of time, and then it makes a turn and starts moving in a new direction for another period of time. Since in practice the path of an animal is usually observed at particular sampling occasions, we consider that the observation process occurs at regular time intervals. Therefore the observation process lies in the location of the individual every time $\Delta t$. As a simplification, we assume that there is no observation error. Our movement model is a form of the velocity jump process (Othmer et al., 1988) where the speed of movement during the active phase is constant, and the temporal scale of the waiting time of the reorientation phase is considered instantaneous. Assuming constant movement speed, let the variable $M_{i}$ describe the position of the latent process at step $i$, presented in x-y coordinates, i.e. $M_{i}=(\mu_{i,1},\mu_{i,2})$ where $i$ represents an index of the time over the steps for $i=0,\dots,N_{steps}$. Given, $\mu_{0,1}=0$ and $\mu_{0,2}=0$, we have for $i=1,...,N_{steps}$ that,

[TABLE]

where $t_{i}$ is the duration of step $i$, and $\omega_{i}$ is the turning angle between steps $i$ and $i+1$, so that $\phi_{i}$ represents the direction of the step $i$. Each $t_{i}$ is assumed to be independently drawn from an exponential distribution with parameter $\lambda$ and each $\omega_{i}$ from a von Mises distribution with a fixed mean $\nu=0$ and parameter $\kappa$ for the concentration. While the model can be extended to allow $\kappa$ and $\lambda$ to depend on the landscape, environment, or animal behaviour, for this work we only consider the initial case as an starting point.

Next, we define the observation process and its links with the latent movement process. Let $O_{j}=(o_{j,1},o_{j,2})$ denote the position of observation $j$ in x-y coordinates, with $j=0,\dots,N_{obs}$. A second index $j$ is used for the time over the observations. Equations (2) show the relationship between observational and latent process. For this, it is necessary to determine the number of changes in direction that occurred before a given observation, we define $N_{j}$ as the number of steps (or changes in direction) that the animal took from time $1$ to time $j\Delta t$.

$O_{0}=M_{0}$

[TABLE]

Then, $O_{j}=h(M_{0:N_{j}})$, a function of all positions $M_{i}$ from $i=0$ to $i=N_{j}$ with $M_{0:D}=(M_{0},M_{1},M_{2},\dots,M_{D})$. Note that $N_{j}$ is the maximum index such that the sum over all the duration times of the steps less or equal to it are at most $j\Delta t$, is possible to express it as $N_{j}=\max\left\{m/\sum_{s=0}^{m}t_{s}\leq j\Delta t\right\}$. Therefore, the location $j$ is the last location of the latent process given by $N_{j}$ plus the difference between $j\Delta t$ and the time at which the step $N_{j}$ was produced in the direction $\phi_{N_{j}}$. To better understand this relationship consider a minimal example of a few steps (Figure 1). Lets assume the duration of steps and turning angles of Table 1 and $\Delta t=0.5$. In that case $N_{1}=1$, because $t_{0}=0.2\leq 1\Delta t$, $t_{0}+t_{1}=0.4\leq 1\Delta t$ but $t_{0}+t_{1}+t_{2}=1.1\nleq 1\Delta t$. With the same reasoning $N_{2}=1$, $N_{3}=3$, $N_{4}=4$, $N_{5}=4$, etc.

Expression for the likelihood function

The likelihood function for the defined process is a function of the number of changes in direction that occurred, $\{N_{j}\}_{j=1}^{N_{obs}}$, and the parameters that determine the underlying movement process.

To construct this expression recall that $O_{j}$ is a function of all positions $M_{i}$ from $i=0$ to $i=N_{j}$, $O_{j}=h(M_{0:N_{j}})$ (Eq. 2). We first suppose that we know the number of changes of direction that the animal took between consecutive observations, i.e. we know $N_{j}\forall j$. Therefore, we can express the likelihood as a function of the underlying movement parameters given both the observational process, $\{o_{j}\}_{j=1}^{N_{obs}}$, and the number of changes in direction, $\{N_{j}\}_{j=1}^{N_{obs}}$, known as the complete-data likelihood.

[TABLE]

In order to derive the complete-data likelihood, $L(\kappa,\lambda;M,O)$, it is necessary to obtain the distributions of $M_{i}$ (Eq. 1) and of $O_{j}=h(M_{0:N_{j}})$ (Eq. 2), which are not available in closed form. For derivation of the marginal likelihood, $L(\kappa,\lambda)$, it is further necessary to integrate over all the possible values of $N_{j}$, by determining $P(N_{j}=r)$ for $r\in\mathbb{N}$, which can be, in principle, infinite. Obtaining the expression for and evaluation of the likelihood results to be a complex task. Likelihood-free methods that allow one to circumvent the need to evaluate the likelihood, such as ABC, have proven to be useful in these cases. es.

Inference Using Approximate Bayesian Computation

Approximate Bayesian Computation (ABC) is a family of simulation-based techniques to obtain posterior samples in models with an intractable likelihood function. In recent years, ABC has become popular in a diverse range of fields (Sisson, Fan, and Beaumont, 2018) such as molecular genetics (Marjoram and Tavaré, 2006), epidemiology (Tanaka, Francis, Luciani, and Sisson, 2006, McKinley, Cook, and Deardon, 2009, Lopes and Beaumont, 2010), evolutionary biology (Bertorelle et al., 2010, Csilléry et al., 2010, Baudet et al., 2015), and ecology (Beaumont, 2010, Sirén et al., 2018). This approach is also useful when the computational effort to calculate the likelihood is large compared to that of the simulation of the model of interest. The likelihood function described earlier turns out to be complex to calculate, yet it is easy to simulate trajectories from the statistical model, based on independent draws from exponential and von Mises distributions combined with the observations at regular time intervals.

Let $\theta$ denote the vector of parameter of interest and $y$ denote the observed data. The posterior distribution $p(\bm{\theta}\mid\bm{y})$ is proportional to the product of the prior distribution $p(\bm{\theta})$ and the likelihood function $p(\bm{y}\mid\bm{\theta})$

$p(\bm{\theta}\mid\bm{y})\propto p(\bm{\theta})p(\bm{y}\mid\bm{\theta})$

The basic idea of ABC methods is to obtain simulations from the joint distribution, $p(\bm{y},\bm{\theta})$ and retain the parameter values that generate synthetic data close to the observed data, $\bm{(}y)$. In this way, ABC methods aim to replace the likelihood function with a measure of similarity between simulated data and actual data.

The rejection algorithm is the simplest and first proposed method to perform ABC (Tavaré, Balding, Griffiths, and Donnelly, 1997, Pritchard, Seielstad, Perez-Lezaun, and Feldman, 1999). It can be described as follows:

Compute a vector of summary statistics with observed data, $S(\bm{y})$. 2. 2.

Simulate parameters $\bm{\theta_{*}}$ sampled from $p(\bm{\theta})$ and synthetic data $\bm{y_{*}}$ sampled from $p(.\mid\bm{\theta_{*}})$. 3. 3.

Compute a vector of summary statistics with simulated data, $S(\bm{y}_{*})$. 4. 4.

$\bm{\theta}_{*}$ is accepted as a posterior sample, if $\rho(S(\bm{y}_{*}),S(\bm{y}))<\delta$, for some distance measure $\rho$ and threshold $\delta$. 5. 5.

Repeat 2-4 $\bm{K}$ times.

The above rejection algorithm produces samples from $p(\bm{\theta}\mid\rho(S(\bm{y}),S(\bm{y}_{*}))<\delta$ which is an approximation of $p(\bm{\theta}\mid\bm{y})$. In particular, when the summary statistics are sufficient or near-sufficient for $\rho$, the approximate posterior distribution converges to the true posterior distribution as $\delta$ goes to [math] (Marjoram, Molitor, Plagnol, and Tavaré, 2003). Instead of selecting a value for $\delta$, it is a common practice to set a threshold $\epsilon$ as a tolerance level to define the proportion of accepted simulations. For a complete review of ABC methods and techniques see (Csilléry, Blum, Gaggiotti, and François, 2010, Beaumont, 2010, Scott. A Sisson and Beaumont, 2018).

For this work we consider two regression-based correction methods. These implement an additional step to correct the imperfect match between the accepted and observed summary statistics. One of these use local linear regression (Beaumont, Zhang, and Balding, 2002), and the other is based on neural networks (Blum and François, 2010). To make the correction, both methods use the regression equation given by

$\theta_{i}=m(S(y_{i}))+\xi_{i}$

where $m$ is the regression function and the $\xi_{i}$’s are centered random variables with equal variance. For the linear correction $m$ is assumed to be linear function and for the neural network correction $m$ is not necessary linear. A weight $K[d(S(y_{i}),S(y_{0}))]$ (for $K$ a statistical kernel) is assigned to each simulation, so those closer to the observed summary statistics are given greater weight. The $m$ and $\xi$ values can be estimated by fitting a linear regression in the first case and a feed-forward neural network regression in the second case. Then, a weighted sample from the posterior distribution is obtained by considering $\theta_{i}^{corr}$ as follows

$\theta_{i}^{corr}=\hat{m}(S(y_{0}))+\hat{\xi}_{i}$

where $\hat{m}(.)$ is the estimated conditional mean and the $\hat{\xi}_{i}$s are the empirical residuals of the regression.

After a preliminary analysis, in which $20$ summary statistics were assessed, we choose four that characterize the trajectories according to parameter values. Looking for summaries that capture diverse features of the movement, we plotted the proposed summaries against known parameters and decided to keep those summary statistics that changed monotonically with parameters values. The plots of all the summaries assessed are provided in the Appendix (Supporting Information). Finally, the four selected summaries were: a point estimator for $\lambda$, calculated as the inverse of the observed average step length (where an observed step is the distance between positions of consecutive observed times); a point estimator for $\kappa$, calculated as the inverse function of the ratio of the first and zeroth order Bessel functions of the first kind evaluated at the mean of the cosine of the observed turning angles (where observed turning angles were determined as the difference between consecutive directions in the observations); the standard deviation of the observed turning angles and lastly, the standard deviation of the observed step lengths (Table 2).

We used the R package “abc” (Csilléry, François, and Blum (2012), http://cran.r-project.org/web/packages/abc/index.html) to perform the analysis. This package uses a standardized Euclidean distance to compare the observed and simulated summary statistics. We present results for the two regression-based correction methods and for the basic rejection ABC method.

Simulations

We did two simulation experiments. First we assessed the performance of the three ABC methods for our model. Then, we evaluated how well these methods approximate posterior probabilities depending on the relation between the temporal scales of simulated trajectories and their observations. For both experiments we used a set of one million simulated trajectories, with parameters $\kappa$ (dispersion parameter for the turning angles) and $\lambda$ (parameter for the times between changes of direction) drawn from the priors $p(\kappa)\sim U[0,100]$ and $p(\lambda)\sim U[0,50]$. The number of simulated steps was such that all the trajectories had at least $1500$ observations. All trajectories were observed at regular times of $\Delta t=0.5$.

Assessment of the inference capacity of the ABC methods

We assessed the performance of the three ABC versions: simple rejection, rejection corrected via linear regression and rejection corrected via neural network. For different values of threshold ($\epsilon$) and for each algorithm version we conducted an ABC cross validations analysis. That is, we selected one trajectory from the million reference set and used it as the real one. We did this selection in a random manner but with the condition that the parameters chosen were not close to the upper limit of the prior distribution. We consider $\lambda\leq 25$ and $\kappa\leq 70$. Then, parameters were estimated using different threshold values ($\epsilon$) with the three algorithms and using all simulations except the chosen one. This process was replicated $N_{rep}=100$ times. For each method and $\epsilon$ value, we recorded the posterior samples obtained for both $\lambda$ and $\kappa$. We then calculated the prediction error as

[TABLE]

where $\theta_{i}$ is the true parameter value of the $i$th simulated data set and $\tilde{\theta}_{i}$ is the posterior median of the parameter. We also compute a dispersion measure of the errors in relation to the magnitude of the parameters for each method and tolerance value. We call $MD$ to this index and we calculate it as

[TABLE]

Furthermore, in order to assess whether the spread of the posterior distributions were not overly large or small, we computed the empirical coverage of the $\alpha=95\%$ credible interval for the two parameters and for different thresholds ($\epsilon$). The empirical coverage is the proportion of simulations for which the true value of the parameter falls within the $\alpha\%$ highest posterior density interval (HPD). If the nominal confidence levels were accurate, this proportion should have been near $0.95$. If this is true for all $\alpha$, it is said that the analysis satisfies the coverage property. A way to test this property is by performing the Coverage Test and it is also a useful way to choose the threshold value $\epsilon$. This test was first introduced by Prangle et al. (2014) and the basic idea is to perform ABC analyses under many data sets simulated from known parameter values and for each of them compute $p$, the proportion of the estimated posterior distribution smaller than the true parameter. Ideally these values should be distributed as a $U(0,1)$. For a complete description of this test see (Prangle et al., 2014). In order to analyze all possible $\alpha$ values, we performed this test using the package “abctools” (Nunes and Prangle (2015); https://cran.r-project.org/web/packages/abctools/index.html).

Relative scale of observations and accuracy of the posterior density

We continued the analysis evaluating how well these methods approximate posterior probabilities as a function of the ratio, $R$, between the temporal scale of observation ($\Delta t$) and the temporal scale for changes in directions ($1/\lambda$) (Figure 2). For instance, if $R=1$ then $\lambda=1/\Delta t$, which means that the time between consecutive observations is equal to the mean time between changes in direction. Conversely, if $R<1$ then the time scale between consecutive observations is smaller than the time scale at which animals decide to change directions, and the opposite occurs if $R>1$ (Figure 2). We considered different values of $R$ (between $0.06$ and $5$) and for each simulated $50$ trajectories with values of $\kappa\in\{10,20,30,\ldots,70\}$. Then, using the original million set, the estimations for the three methods of ABC were computed considering these new trajectories as the true observations . We calculated the predictor error for $\kappa$ and $\lambda$ for every combination of $R$ and $\kappa$.

Real data example

We evaluate the performance of this model using data from a real trajectory reconstructed in high-resolution using information from DailyDiary and GPS devices. With a high resolution trajectory it is possible to infer in such manner the moments in which the animal changes the direction and use that information to estimate the mean of the step length and the dispersion of the turning angles. In this manner it is possible to obtain certain “true” parameter values to then see if they can be recovered using the ABC techniques.

The data used were collected from one sheep in Bariloche, Argentina, during February and March of 2019. The sheep was equipped with a collar containing a GPS (CatLog-B, Perthold Engineering, www.perthold.de; USA), that was programmed to record location data every five minutes, and a DailyDiary (DD, (Wilson et al., 2008)), that was programmed to record $40$ acceleration data per second (frequency of $40$Hz) and $13$ magnetometer data per second (frequency of $13$Hz). The DD are electronic devices that measure acceleration and magnetism in three dimensions, which can be described relative to the body of the animal. Such data allow the Dead-Reckoned (DR) path of an animal to be reconstructed at high resolution. The goal here is to use the detailed observed trajectories to decompose them in steps and turns and then test if the ABC approach can estimate the parameters of the distributions for steps and turns.

From the original data, we randomly selected one segment of $6$ hours. Using the DD information, we first estimated the path traveled (pseudotrack) by the sheep using the dead-reckoning technique. (Wilson and Wilson, 1988, Wilson et al., 2007). In this step we made use of the R package “TrackReconstruction” (https://cran.r-project.org/web/packages/ TrackReconstruction/index.html). After that, we corrected the bias of those estimation using the data from the GPS (Liu, Y. et al., 2015). This correction was made using the R package “BayesianAnimalTracker” (https://cran.r-project.org/web/packages/ BayesianAnimalTracker/index.html). In this way, we obtained a trajectory sampled with a resolution of 1 second. To satisfy the hypotheses of the model, we selected part of that trajectory that appeared to come from the same behaviour, i.e we selected a piece of the trajectory that visually appeared to have the same distribution of turn angles and step lengths.

In order to estimate the parameters of the trajectory it was necessary to determine the points in which there was a change of movement direction. We applied the algorithm proposed by Potts. et al. (Potts et al., 2018) that detects the turning points of the trajectory using data of the animal headings and subsequently calculated steps and turning angles. So, in that way we not only obtained the values of the $N_{j}$, but also we obtained samples for the step’s length and the turning angles. With that information it was really simple to infer the parameter’s values via MCMC techniques and obtain samples from the joint posterior distribution using the Stan software (Carpenter et al., 2017).

Then, we calculated the summary statistics of the trajectory observed at $dt=50$ secs ($1$ observation every $50$ of the reconstructed trajectory) and applied the three ABC algorithms. We finally compared both estimations.

Results

Assessment of the inference capacity of the ABC methods

Figure 3 shows the values of the prediction errors and the dispersion index ($MD$) for each method and $\epsilon$ value obtained from the ABC cross validation analysis. In all cases the prediction errors decreased when the value of the threshold ($\epsilon$) decreased . However, for the algorithms corrected via linear regression and neural networks larger threshold levels ($\epsilon$) can produce lower prediction error. Something similar happened with the $MD$ values: lower threshold values imply lower values of this index (Figure 3). These values give us an idea of the width of the posterior distributions. It is evident that for the case of the rejection algorithm the posterior distributions are quite wide, especially for $\epsilon=0.1$. However, for the corrected algorithms we can assume that the difference between the estimated and true parameters are up to approximately $1.3$ units for $\kappa$ and in $0.3$ for $\lambda$ in the best case ($\epsilon=0.001$). Figure 4 shows that the estimate of $\lambda$ improves when it takes lower values, especially for the algorithm corrected via linear regression. We discuss this point in the Discussion. Based on these results, the algorithm corrected via linear regression seems to be the one that presents the best performance.

We estimated the empirical coverage of the $95\%$ HPD intervals for both parameters ($\kappa$ and $\lambda$), for the three ABC algorithms and for different threshold ($\epsilon$) values. Almost always these indices were greater than $95$, except in the case of the highest threshold value ($\epsilon=0.1$) for the simple rejection algorithm and the one corrected via neural network, for which the empirical coverages were a little below $0.95$. The plots for this analysis are provided in the Appendix.

Finally, in order to check the coverage property we performed a coverage test for both parameters (Figure 5). In most cases, the distributions obtained do not show a clear approximation of a $U(0,1)$. However there is an evident difference between the histograms obtained with the simple rejection ABC and the histograms obtained with the other two algorithms. The shapes of the rejection ABC are those that are farthest from being uniform: for both parameters the distributions of the p values are left skewed indicating that the algorithm tends to overestimate the parameters. For the other two algorithms the left skewed is much moderate, and even in the case of the lowest $\epsilon$ values for the linear algorithm the histograms are more uniform, indicating that its coverage could be being reached. However, not rejecting that coverage holds does not unequivocally demonstrate that the ABC posterior approximation is accurate. If the empirical data is uninformative, the ABC will return posterior distributions very similar to the priors, which would produce uniform coverage plots.

Relative scale of observations and accuracy of the posterior density

In order to evaluate the importance of the relationship between the time scale of the observation process and the time in which changes occur in the latent process, we evaluated how well the two parameters fit in relation to the R ratio. The prediction errors for $\lambda$ increased as the value of $R$ increased (Figure 6). For the case of the prediction errors for $\kappa$ this relation can be seen when the true value of this parameter takes larger values (Figure 7).

For the case of the prediction errors for $\kappa$ this relation can be seen when the true value of this parameter takes larger values. Again, the corrected algorithms have the smallest errors for both parameters.

According to the above results, it is evident that there is a relationship between the ratio $R$ and the capacity of these methods to estimate the parameters. For rates approximately less than $5$ the errors are small and it is possible to obtain good estimates. This necessitates that the time scale of the observation process be approximately less than $5$ times the time-scale at which the animals decide to change direction. For higher values of $\Delta t$ it would be more difficult to make inferences using this technique.

Sheep data

The selected trajectory was of February $27$, $2019$ from $19:01:21$hs to $20:02:00$hs, a total of $1.01$ hours (Figure 8). The posterior distribution of each parameter was estimated from a sample of $3\times 1000$ independent MCMC observations. As they were in the lower limits of the prior distributions we simulate a new set of trajectories with priors of $\kappa\sim U[0,10]$ and $\lambda\sim U[0,10]$ and compute the ABC inferences with those simulations.

The posteriors obtained through MCMC and through the ABC algorithms gave similar results (9). Again, the rejection ABC algorithm produced the estimation which is less exact, i.e the posterior is the furthest from the one obtained by MCMC. Although this trajectory is just a simple example, it shows that it is possible to apply this model to actual animal trajectories.

Discussion

Animal movement modelling and analysis is either considered in continuous or discrete time. Continuous time models are more realistic but often harder to interpret than the discrete versions (McClintock et al., 2014). A compromise between these approaches is to model movement as steps and turns but to have step duration (or the times at which turns are made) occurring in continuous time. In this way one could get the best of both worlds, so to speak. Here we considered that the underlying movement process, evolving in continuous time, is observed at regular time intervals as would be standard for a terrestrial animal fitted with a GPS collar. The likelihood function resulted to be complex to calculate, but it is feasible to quickly generate simulations from the process and observation models. Thus, we proposed to use ABC methods. Even though these techniques showed certain limitations, it was possible to obtain accurate parameter estimates when the temporal scale of observations was not too coarse compared to the scale of changes in direction.

Our simulation study showed that simple rejection ABC does not perform well for the proposed state-space model but the two corrected version of this algorithm really improve the estimations (Table 3 and Figure 4). Overall, the best performance was obtained with the linear correction. However, the applicability of these methods depends strongly on the rate between the observation process’s scale and the mean time between changes in movement direction. We found that when this ratio is smaller than $5$ it is possible to make inferences about the parameters (Figures 6 and 7). That is, it would be necessary that the observations are less than $5$ times the average of the times between changes of directions in order to be able to generate good estimations.

Beyond our findings about the capacity to make inference with these techniques in a simulation study, it is important to note that in an applied case more informative priors could be considered. Here, our aim was to evaluate the performance of the ABC techniques considering several parameter combinations generating trajectories and then sampling from those trajectories. In order to optimize computing time, we simulated a million trajectories sampling their parameters from uniform distributions and then we randomly choose one of them as observed data while the rest of the simulations was used to perform the ABC computations. That justify the use of uniform priors for our parameters. As we did in our real data example, in applied cases it wold be relatively straightforward to come up with more informative priors, especially for the expected time for changes in movement direction.

The movement model presented here is quite simple as we assume constant movement speed and turning angles with zero mean. Nevertheless, the model is an improvement over discrete-time versions where the temporal scale of movement has to match the scale of observations. Further developments of these methods should consider additional features that are common in movement studies such as including more than one movement behavior and the effect of habitat features on both movement parameters and changes between behaviors (Morales et al., 2004, Mevin B. Hooten et al., 2017). Even though this means estimating several parameters, such models will imply further structure in the trajectories that could be used as part of the summary statistics used to characterize the data. Hence, it might reduce the combination of parameters values capable of reproducing features present in the data, allowing for ABC inference.

As new technologies allow us to obtain very detailed movement data, we can have better estimates of the temporal scales at which animals make movement decisions. As we did in our real data example, high-frequency data from accelerometers and magnetometers combined with GPS data can be used to obtain trajectories with sub-second temporal resolution to then detect fine-scale movement decisions such as changes in direction. These detailed trajectories could be used to elicit informative priors to use when only coarser data is available.

In general, the processes behind the realized movement of an individual and the processes that affect how we record the trajectory are usually operating at different time scales, making it challenging to analyze and to understand the former given the latter. The state-space model used here allowed us to connect these two scales in an intuitive and easy to interpret way. Our findings underscore the idea that the time scale at which animal movement decisions are made needs to be considered when designing data collection protocols, and that high-frequency data may not be necessary to have good estimates of certain movement processes.

Acknowledgements

The authors are grateful to Agustina Di Virgilio, Teresa Moran Lopez and Florencia Tiribelli for their useful comments that substantially improved the manuscript.

Authors’ Contributions

JMM conceived the ideas; All the authors designed methodology; SRS simulated the data and performed the analysis; SRS and JMM led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

Supporting Information

Calculation of the complete data likelihood

Lets consider the variable $M_{i}=(\mu_{i,1},\mu_{i,2})$ for describing the position (in $x-y$ coordinates) of the latent process by step $i$ and the variable $O_{j}=(o_{j,1},o_{j,2})$ for the position of the observation $j$. Lets remember that we defined $N_{j}$ as the amount of steps (or changes of direction) that the animal took from time $1$ to time $j(\Delta t)$.

We have that $\mu_{0,1}=0$ and $\mu_{0,2}=0$

For $i=1,...,N_{steps}$

$\mu_{i,1}=\mu_{i-1,1}+cos(\phi_{i-1})t_{i-1}$

$\mu_{i,2}=\mu_{i-1,2}+sin(\phi_{i-1})t_{i-1}$

And then it is possible to parameterize the observation process as

$o_{0,1}=0$ and $o_{0,2}=0$ and for $j=1,...,N_{obs}$

$\displaystyle o_{j,1}=\mu_{N_{j},1}+cos(\phi_{N_{j}})\left(j\Delta t-\sum_{k<N_{j}-1}t_{k}\right)$

$\displaystyle o_{j,2}=\mu_{N_{j},2}+sin(\phi_{N_{j}})\left(j\Delta t-\sum_{k<N_{j}-1}t_{k}\right)$

So $o_{j}$ is a function of all positions $M_{i}$ from $i=0$ to $i=N_{j}$. Then $O_{j}=h(M_{0:N_{j}})$, where $M_{0:D}=(M_{0},M_{1},M_{2},...,M_{D})$. Lets suppose that we know the number of changes of direction that the animal took between consecutive observations, so we know the $N_{j}\forall j$. Then the likelihood of the SSM with different temporal scales for an individual trajectory is given as

[TABLE]

As $O_{j}=h(M_{0:N_{j}})$, in order to get a formulation of $L$ it is necessary to obtain the distributions of $M_{i}$ $(1)$ and for $O_{j}=h(M_{0:N_{j}})$ $(2)$.

Step 1: Formulation of $(1)$

We are looking for a formulation for $M_{i}=(\mu_{i,1},\mu_{i,2})$. We are going to consider just the variable corresponding to the x-coordinate ($\mu_{i,1}$), the second is analogous.

We have that

[TABLE]

with $\phi_{i}\sim vonMises(\phi_{i-1},\kappa)$ and $t_{i}\sim Exp(\lambda)$. To obtain the distribution of $\mu_{i}|\mu_{i-1}$ it is necessary to obtain the distribution form of $Z=cos(\phi)t$. Using the Change of variable Theorem it is possible to calculate this distribution. To do that, lets first consider $V=g(\phi)=cos(\phi)$. We want to obtain an expression for $f_{V}$. Splitting the domain of $g$ and applying the Transformation Method Theorem, is obtain:

[TABLE]

Now we can calculate $f_{Z}$ as $f_{Z}=Vt$. Again making use of the Transformation Method Theorem and using the fact that the times and angles are independent, it is possible to obtain the following expression

[TABLE]

Having $f_{Z}$ obtaining $p(\mu_{i}|\mu_{i-1})$ is immediate.

Step 2: Formulation of $(2)$

Now, we are looking for a formulation for $O_{j}=h(M_{0:N_{j}})$. Les write

$O_{j}=(o_{j,1},o_{j,2})=\left(h_{1}(\mu_{0:N_{j},1}),h_{2}(\mu_{0:N_{j},2})\right)$

Again, we are going to consider just the variable corresponding to the x-coordinate ($o_{j,1}$), the second is analogous.

We have that

[TABLE]

We already know the distribution of $V_{N_{j}}$. The distribution of $W_{N_{j}-1}$ is just a sum of $N_{j}-2$ $Exp(\lambda)$, a $\Gamma(N_{j}-2,\lambda)$. If we consider $\tilde{W}_{N_{j}-1}=c_{j}-W_{N_{j}-1}$ (which differs with $W_{N_{j}-1}$ just in a constant), we have that $f_{\tilde{W}_{N_{j}-1}}(v)=f_{W_{N_{j}-1}}(c_{j}-v)$

So, we can rewrite $o_{j,1}$ as

[TABLE]

To obtain the distribution of $S_{N_{j}-1}$, again is necessary to use the Transformation Method Theorem and the independence between the times and the angles:

[TABLE]

Summary Statistics

We provide the plots of the summary statistic analyzed versus the parameters (Figures 10, 11, and 12). We choose four that attempt to describe the trajectories in an integral way and characterize them according to parameter values. The selected were 10(a), 10(b),10(c) and 10(d)

Empirical Coverage

We present the results for the empirical coverage of the $95\%$ high posterior density(HPD) intervals for the two parameters. This value is the proportion of simulations for which the true parameter value falls within the $95\%$ HPD interval. If the posterior distributions were correctly estimated, this proportion should have been near $0.95$. We compute this index for both parameters ($\kappa$ and $\lambda$) and for the three ABC algorithms: Simple Rejection, Corrected via Linear Regression and Corrected via Neural Network. We did that for threshold ($\epsilon$) values of: $0.001$, $0.005$, $0.01$ and $0.1$.