Using approximate Bayesian inference for a “steps and turns” continuous-time random walk observed at regular time intervals

Movement model with random time between movement decisions

Many animal movement studies assume discrete time correlated random walks (CRW) as building blocks to more complex models (Turchin, 1998; Morales et al., 2004). At every time step, an individual chooses a turning angle (difference between the previous and new movement direction) and a speed. The turning angle distribution is concentrated around zero, resulting in short-term correlations in the direction of movement (persistence in movement direction). Here, we allow changes in direction to occur at any point in time. 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. Our movement model is a form of the velocity jump process in which the motion consists of a sequence of ”runs” separated by reorientations, during which a new velocity is chosen (Othmer, Dunbar & Alt, 1988; Codling, Plank & Benhamou, 2008). In our model, 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.

Since, in practice, the path of an animal is usually observed at particular sampling occasions, we consider that the observation occurs at regular time intervals. Therefore the observation lies in the location of the individual every time Δt. As a simplification, we assume that there is no observation error.

Assuming constant movement speed, let the variable M_i describes the position of the state process at step i, presented in x–y coordinates, i.e.,  M_i = (m_i,1, m_i,2) where i represents an index of the time over the steps for i = 0, …, N_steps. Given, m_0,1 = 0 and m_0,2 = 0, we have for i = 1, …, N_steps that, (1)${\displaystyle \begin{array}{cc}\hfill m_{i,1}=\hfill & \hfill m_{i-1,1}+cos\left({\varphi }_{i-1}\right)t_{i-1}\hfill \\ \hfill m_{i,2}=\hfill & \hfill m_{i-1,2}+sin\left({\varphi }_{i-1}\right)t_{i-1}\hfill \\ \hfill {\varphi }_i=\hfill & \hfill \sum _{k=1}^i{\omega }_i\hfill \end{array}}$ where t_i is the duration of step i, and ω_i is the turning angle between steps i and i + 1, so that ϕ_i represents the direction of the step i. Each t_i is assumed to be independently drawn from an exponential distribution with parameter λ and each ω_i from a von Mises distribution with a fixed mean ν = 0 and parameter κ for the concentration  (Wu et al., 2000; Codling, Plank & Benhamou, 2008). While the model can be extended to allow κ and λ 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 and its links with the state movement process (Eq. (2)). Let O_j = (o_j,1, o_j,2) denote the position of observation j in x–y coordinates, with j = 0, …, N_obs. A second index j is used for the time over the observations. We defined T_i as the time in which the change of direction i took place: ${\displaystyle \begin{array}{c}{T}_0=0\hfill \\ T_i=\sum _{k=0}^{i-1}{t}_k\text{for}i=1,\dots ,N_{steps}.\hfill \end{array}}$

In addition 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Δt.

O₀ = M₀

And for j = 1, …, N_obs (2)${\displaystyle \begin{array}{c}\hfill o_{j,1}=m_{N_j,1}+cos\left({\varphi }_{N_j}\right)\left(j\Delta t-T_{N_j}\right)\hfill \\ \hfill o_{j,2}=m_{N_j,2}+sin\left({\varphi }_{N_j}\right)\left(j\Delta t-T_{N_j}\right).\hfill \end{array}}$

Note that N_j is the index that corresponds to the maximum time T_i less or equal to jΔt, i.e., $N_j=max\left\{m∕T_m\le j\Delta t\right\}$. Therefore, the location j is the last location of the state process given by N_j plus the difference between jΔt and the time at which the step N_j was produced in the direction ϕ_Nj. To better understand this relationship consider a minimal example of a few steps. Assuming the duration of steps and turning angles of Table 1 and Δt=0 .5, in Fig. 1 we present a short trajectory showing the locations of changes in direction and observations at regular times. In that case N₁ = 2, because T₁ = t₀ = 0.2 ≤ 1Δt, T₂ = t₀ + t₁ = 0.4 ≤ 1Δt but T₃ = t₀ + t₁ + t₂ = 1.1⁄ ≤ 1Δt. With the same reasoning N₂ = 2, N₃ = 4, N₄ = 5, N₅ = 5, etc.

Expression for the likelihood function

In order to construct the likelihood we first present the complete-data likelihood, i.e., we assume that we know all the values for t_i and ω_i. In that case, as O|M is deterministic, we have ${\displaystyle L\left(\kappa ,\lambda ;M,O\right)=P\left(O,M|\kappa ,\lambda \right)=P\left(O|M,\kappa ,\lambda \right)P\left(M|\kappa ,\lambda \right)=P\left(M|\kappa ,\lambda \right)=P\left(t_1,...,t_{N_{steps}}|\lambda \right)P\left({\omega }_1,...,{\omega }_{N_{steps}}|\kappa \right)}$

Suppose that we do not know the values of t_i and ω_i, but do know the number of steps that the animal took between consecutive observations, N_j∀j. In that case, it would be necessary to obtain the distributions of M_i (Eq. (1)) (or at least of a proportion of them), which are not available in closed form.

Finally, to formulate the marginal likelihood, L(κ, λ), it is further necessary to integrate over all the possible values of N_j, by determining P(N_j = r) for r ∈ ℕ, 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.

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 & Beaumont, 2018) such as molecular genetics (Marjoram & Tavar, 2006), epidemiology (Tanaka et al., 2006; McKinley, Cook & Deardon, 2009; Lopes & Beaumont, 2010), evolutionary biology (Bertorelle, Benazzo & Mona, 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 θ denote the vector of parameter of interest and y denote the observed data. The posterior distribution p(θ∣y) is proportional to the product of the prior distribution p(θ) and the likelihood function p(y∣θ)

p(θ∣y) ∝ p(θ)p(y∣θ)

The basic idea of ABC methods is to obtain simulations from the joint distribution, p(y, θ) and retain the parameter values that generate simulated data close to the observed data, (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é et al., 1997; Pritchard et al., 1999). It can be described as follows:

1. Compute a vector of summary statistics with observed data, S(y).

2. Simulate parameters θ_∗ sampled from p(θ) and data y_∗ sampled from p(.∣θ_∗).

3. Compute a vector of summary statistics with simulated data, S(y_∗).

4. θ_∗ is accepted as a posterior sample, if ρ(S(y_∗), S(y)) < δ, for some distance measure ρ and threshold δ.

5. Repeat 2–4 K times.

The above rejection algorithm produces samples from p(θ∣ρ(S(y), S(y_∗))) < δ which is an approximation of p(θ∣y). In particular, when the summary statistics are sufficient or near-sufficient for ρ, the approximate posterior distribution converges to the true posterior distribution as δ goes to 0 (Marjoram et al., 2003). Instead of selecting a value for δ, it is a common practice to set a threshold ϵ as a tolerance level to define the proportion of accepted simulations. For a complete review of ABC methods and techniques see (Csilléry et al., 2010; Beaumont, 2010; Sisson, Fan & Beaumont, 2018).

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 & Balding, 2002), and the other is based on neural networks (Blum & François, 2010). To make the correction, both methods use the regression equation given by

θ_i = r(S(y_i)) + ξ_i

where r is the regression function and the ξ_i’s are centered random variables with equal variance. For the linear correction r is assumed to be linear function and for the neural network correction r is not necessary linear. A weight K[d(S(y_i), S(y₀))] (for K a statistical kernel) is assigned to each simulation, so those closer to the observed summary statistics are given greater weight. The r and ξ 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{r}\left(S\left(y_0\right)\right)+{\hat{\xi }}_i$

where $\hat{r}\left(.\right)$ 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 S1. Finally, the four selected summaries were: 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 κ, calculated as the inverse function of the ratio of the first and zero 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, Franois & 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 κ (dispersion parameter for the turning angles) and λ (parameter for the times between changes of direction) drawn from the priors p(κ) = U[0, 100] and p(λ) = U[0, 50]. The choice of uniform priors was because the aim of our simulation study is to analyze the performance of the ABC methods to fit the model under a diversity of scenarios. Using uniform priors it is possible to consider different values for the parameters with equal probability. The number of simulated steps was such that all the trajectories had at least 1500 observations. All trajectories were observed at regular times of Δt=0 .5. The R code is available in https://github.com/sofiar/ABC-steps-and-turns-.

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 (Δt) and the temporal scale for changes in directions (1∕λ) (Fig. 2). For instance, if R = 1 then λ = 1∕Δ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 (over-sampling case), and the opposite occurs if R > 1 (sub-sampling case) (Fig. 2). We considered different values of R (between 0.06 and 5) and for each simulated 50 trajectories with values of κ ∈ {10, 20, 30, …, 70}. Then, using the original million trajectories, the estimations for the three methods of ABC were computed considering these new trajectories as the true observations. We calculated the predictor error for κ and λ for every combination of R and κ.

Sheep data example

In the case that there is only location data available, i.e., GPS data at each time Δt, we have concluded that as it is not possible to know how many changes of directions were between consecutive observations, the likelihood function resulted to be intractable, making ABC techniques useful methods to make inference about the model. Suppose that there exists a source of information that permits estimation of the values of t_i and ω_i ∀i. Considering these estimated values as the truth, we have samples t₁, …, t_n and ω₁, …, ω_n and therefore the estimation of κ and λ is direct.

We compared both estimations for a real trajectory. As we have presented thus far, we fit the model using only information from the GPS locations with ABC techniques. In addition, we used information from DailyDiary devices in order to reconstruct the high-resolution trajectory. This allowed us to infer the times at which the animal changes direction, obtain samples of t and ω and use MCMC to obtain draws from the posterior distributions of the parameters. 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, http://www.perthold.de; USA), that was programmed to record location data every five minutes, and a DailyDiary (DD, Wilson, Shepard & Liebsch, 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. In some cases there is no such data, and only information of the location at certain times is available. The goal here is to compare the ability of the ABC’s methods to estimate the parameters of the model when only information of the location at certain times is available, with the estimation of MCMC techniques when extra information is available.

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 & 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 the estimations using the data from the GPS (Liu et al., 2015). This correction was made using the R package “BayesianAnimalTracker” (https://cran.r-project.org/web/packages/BayesianAnimalTracker/index.html). In this manner, we obtained a trajectory sampled with a resolution of 1 s. 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.

To determine the points at which there was a change of movement direction, we applied the algorithm proposed by 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. With that information it is straightforward to infer the parameter’s values via MCMC techniques and obtain samples from the joint posterior distribution directly using the software Stan as the likelihood can be written in closed form (Carpenter et al., 2017).

We also 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. The data used is available at Figshare with DOI (10.6084/m9.figshare.9971642.v2).

Assessment of the inference capacity of the ABC methods

Figure 3 shows the values of the prediction errors and the dispersion index (DI) for each method and ϵ value obtained from the ABC cross validation analysis. In all cases the prediction errors decreased when the value of the threshold (ϵ) decreased. However, for the algorithms corrected via linear regression and neural networks larger threshold levels (ϵ) can produce lower prediction errors. Something similar happened with the DI values: lower threshold values imply lower values of this index (Fig. 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 ϵ = 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 κ and in 0.3 for λ in the best case (ϵ = 0.001). We also analyzed the relationship between the true parameter values and the median of the estimated posterior (Fig. 4). The estimate of λ improves when it takes lower values, especially for the algorithm corrected via linear regression. That behaviour is not so clear for the parameter κ. We further elaborate this point in the discussion. Based on these results, the algorithm corrected via linear regression seems to perform the best.

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

Finally, in order to check the coverage property we performed a coverage test for both parameters (Fig. 5). In most cases, the distributions obtained do not show a clear approximation to a U(0, 1). However there is an evident difference between the histograms obtained with the simple rejection ABC and those 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 in the case of the lowest ϵ 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 and the time at which changes occur in the movement process, we evaluated how well the two parameters fit in relation to the R ratio. The prediction errors for λ increased as the value of R increased (Fig. 6). For the case of the prediction errors for κ this relation can be seen when the true value of this parameter takes larger values (Fig. 7).

For the case of the prediction errors for κ 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 results shown in Figs. 6 and 7, 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 be approximately less than 5 times the time-scale at which the animals decide to change direction. For higher values of Δt it would be more difficult to make inferences using this technique.

Sheep data

The selected trajectory was from February 27th, 2019 from 19:01:21hs to 20:02:00hs, a total of 1.01 hours (Fig. 8). For the estimations made using the reconstruct high resolution trajectory, the posterior distribution of each parameter was estimated from a sample of 3 × 1000 MCMC draws. As they were in the lower limits of the prior distributions, to conduct inference via ABC with only the location data, we simulate a new set of trajectories with priors of p(κ) = U[0, 10] and p(λ) = U[0, 10].

Draws from the posterior distributions obtained through MCMC and through the ABC algorithms gave similar results (Fig. 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.