Metrics for describing dyadic movement: a review

In the following subsections the metrics are defined and their theoretical properties are described. A summary is proposed in Table 1. Considering two individuals named A and B, the position of A (resp. B) at time t is denoted by \(X_{t}^{A}\) (resp. \(X_{t}^{B}\)). The distance between A at time t₁ and B at time t₂ will be referred to as \(d_{t_{1},t_{2}}^{A,B}\). When the distance between two individuals is regarded at simultaneous time, this will be shortened to \(d_{t}^{A,B}\). Whenever possible, metrics introduced by different authors but that are actually very similar in their definition, are grouped under a unified name and a general definition.

Proximity index (Prox)

The proximity index (Prox in [5]) is defined as the proportion of simultaneous pairs of fixes within a distance below an ad hoc threshold (Fig. 1). Other metrics in the literature are actually analogous to Prox: the coefficient of association (Ca) [12] and the I_AB index [4]. Denoting by T the number of pairs of fixes in the dyad, we propose a unified version of those metrics using a kernel K (formula 1):

$$ Prox_{K,\delta} = \frac{1}{T} \sum\limits_{t=1}^{T} K_{\delta}\left(X_{t}^{A}, X_{t}^{B}\right), $$

(1)

where δ is a distance threshold parameter.

Choosing \(K_{\delta }(x,y) =\mathbbm {1}_{\{ \| x-y \| < \delta \}}\) (\(\mathbbm {1}_{\{\}}\) represents the indicator function) as a kernel leads to the Prox metric in [5], denoted by \(Prox_{\mathbbm {1},\delta }\) henceforward. Instead, choosing K_δ(x,y)= exp(−∥x−y∥²/(2δ²)) gives the I_AB index. Regarding Ca, for simultaneous fixes, its definition becomes exactly the same as \(Prox_{\mathbbm {1},\delta }\) (using Ca’s adaptation to wildlife telemetry data shown in [38]).

Most of the proximity-related metrics are based on symmetric kernels and depend only on the distance between A and B; therefore, the formula notation (1) can be simplified as:

$$ Prox_{K,\delta} = \frac{1}{T} \sum\limits_{t=1}^{T} K_{\delta}\left(X_{t}^{A}, X_{t}^{B}\right) = \frac{1}{T} K_{\delta}^{+}. $$

(2)

If the distance between two individuals is below the threshold δ during their whole tracks, \(Prox_{\mathbbm {1},\delta }\) will be 1 (and 0 in the opposite case). \(Prox_{\mathbbm {1},\delta }\) might be interpreted as the proportion of time the two individuals spent together. This interpretation is, of course, threshold dependent. The I_AB index provides a smoother measure of the average proximity between two individuals along the trajectory. Proximity is thus dependent on the choice of a δ parameter and of a kernel function. Graphical examples illustrating the differences in \(K_{\delta }(x,y) =\mathbbm {1}_{\{ \| x-y \| < \delta \}}\) and K_δ(x,y)= exp(−∥x−y∥²/(2δ²)) are in Additional file 1.

The Half-weight Association Index (HAI)

The Half-weight Association Index (HAI) proposed by [10] measures the proportions of fixes where individuals are close to each other (within a user-defined threshold). By that definition, HAI is exactly the same as \(Prox_{\mathbbm {1},\delta }\). However, HAI was popularized by [2] in another form that did not consider all fixes for the computation of the metric, but used counts with respect to a reference area (called overlapping zone in the original paper):

$$ HAI = \frac{K_{\delta}^{+}}{K_{\delta}^{+} + \frac{1}{2}(n_{A0}+n_{0B})} $$

(5)

where n_AB (resp n_A0; n_0B; n₀₀) is the number of simultaneous occurrences of A and B in the reference area S_AB (resp. simultaneous presence of A and absence of B; simultaneous absence of A and presence of B; simultaneous absence of A and absence of B), and where \(K_{\delta }^{+}\) is computed over the reference area.

It is worth noticing that the HAI adaptation proposed by [2] does not correctly account for spatial joint movement, as would do a \(Prox_{\mathbbm {1},\delta }\) version constraint to the reference area; i.e. the denominator should be equal to n_AB + n_A0 + n_0B, which is the total number of simultaneous fixes where at least one individual is in the reference area.

The dependence to the definition of an overlapping zone or reference area is discussed in the following subsection dedicated to L_ixnT, which also relies on the definition of a static reference area.

If the individuals remain together (i.e. in the reference area and closer than δ) all the time, HAI is close to 1, and 0 in the opposite case. An example of the computation of HAI under [2]’s definition is given in Fig. 2.

Joint Potential Path Area (jPPA)

Long et al. [39] computed the relative size of the potential encounter area at each time step of two individuals’ tracks. Assuming a speed limit ϕ, the potential locations visited between two consecutive fixes define an ellipse (Additional file 4). Then, the potential encounter area corresponds to the intersection between the ellipses of the two individuals (at simultaneous time steps; see Fig. 3). The overall potential meeting area is given by the spatial union of all those potential encounter areas. This area is then normalized by the surface of the spatial union of all the computed ellipses to produce the joint Potential Path Area (jPPA) metric ranging from 0 to 1 (see formula in Table 1). jPPA values close to 0 indicate no potential spatio-temporal overlap, while values close to 1 indicate a strong spatio-temporal match.

Several issues can be discussed here. First, no movement model is assumed and therefore the method confers the same probabilities of presence to every subspace within the ellipse regions. This is clearly unrealistic as individuals are more likely to occupy the central part of the ellipse because they cannot always move at ϕ, i.e. maximal speed. Second, the computation of the ellipses relies strongly on the ϕ parameter. If ϕ is unrealistically small, it would be impossible to obtain the observed displacements and the ellipses could not be computed. By contrast, if ϕ is too large, the ellipses would occupy such a large area that the intersected areas would also be very large (hence a large jPPA value). Alternatively, [36] proposed a dynamic computation of ϕ as a function of the activity performed by the individual at each fix. Within this approach, additional information or knowledge (i.e. other data sources or models) would be required for the computation of ϕ.

Cross sampled entropy (CSE and CSEM)

Cross sampled entropy (CSE) [51] comes from the time series analysis literature and is used for comparing pairs of motion variables [3; 18, e.g.]. It evaluates the similarity between the dynamical changes registered in two series of any given movement measure. Here we present a simplification of the CSE for simultaneous fixes and position series. A segment of track A would be said to be m-similar to a segment of track B if the distance between paired fixes from A and B remain below a certain threshold during m consecutive time steps. If we define N_m as the number of m-similar segments within the series, then CSE can be defined as (the negative natural logarithm of) the ratio of N_m+1 over N_m and might be understood as (the negative natural logarithm of) the probability for an m-similar segment to also be (m+1)-similar. Formally, CSE is defined as:

$$ \begin{aligned} CSE_{\delta}(m) &= -\ln\left\{\frac{ \sum\nolimits_{t=1}^{T-m} \mathbbm{1}{\left\{ \left(\max_{k \in [0,m]} \left| X_{t+k}^{A} - X_{t+k}^{B} \right| \right) < \delta \right\} } }{\sum\nolimits_{t=1}^{T-m} \mathbbm{1}{\left\{ \left(\max_{k \in [0,m-1]} \left| X_{t+k}^{A} - X_{t+k}^{B} \right| \right) < \delta \right\} }}\right\}\\ &=-\ln\frac{N_{m+1}}{N_{m}}, \end{aligned} $$

(8)

A large value of CSE corresponds to greater asynchrony between the two series, while a small value corresponds to greater synchrony.

CSE relies on an ad hoc choice of both m and δ. In practice, it is expected that the movement series of A and B will not be constantly synchronous and that, for a large value of m, N_m could be equal to 0, in which case CSE would tend to ∞. Therefore, the largest value of m such that N_m>0, i.e. the length of the longest similar segment, could be an alternative indicator of similarity between the series (do not confuse with the longest common subsequence LCSS; see [60]). We propose to use this measure (standardized by T−1 to get a value between 0 and 1) as an alternative index of joint movement (formula 9), which we denote by CSEM. An example of a dyad and the computation of its CSEs and CSEM is shown in Fig. 4.

$$ CSEM = \frac{\max \left\{m; N_{m} >0\right\} }{T-1}, $$

(9)

with the convention that max{∅}=0.

Dynamic Interaction (DI, D I _d and D I _θ)

Long and Nelson [37] argued that it is necessary to separate movement patterns into direction and displacement (i.e. distance between consecutive fixes or step length), instead of computing a correlation of locations [55] which may carry a mixed effect of both components. To measure interaction in displacement, at each time step, the displacements of simultaneous fixes are compared (formula 10).

$$ g_{t}^{\beta} = 1- \left(\frac{\left| d_{t,t+1}^{A} - d_{t,t+1}^{B} \right|}{d_{t,t+1}^{A} + d_{t,t+1}^{B}}\right)^{\beta} $$

(10)

where β is a scaling parameter meant to give more or less weight to similarity in displacement when accounting for dynamic interaction. As β increases, \(g_{t}^{\beta }\) is less sensitive to larger differences in displacement. Its default value is 1. When \(d_{t,t+1}^{A} = d_{t,t+1}^{B}, g_{t}^{\beta }=1\); and when the difference in displacement between A and B at time t is large, \(g_{t}^{\beta }\) approaches zero. For \(g_{t}^{\beta }\) to be 0, one (and only one) of the individuals in the dyad should not move; for a sum of \(g_{t}^{\beta }\) to be equal to zero, at every time t one of the two individuals should not move.

Interaction in direction is measured by

$$ f_{t} = \cos\left(\theta_{t,t+1}^{A} - \theta_{t,t+1}^{B}\right) $$

(11)

where θ_t,t+1 is the direction of an individual between time t and t+1. f_t is equal to 1 when movement segments have the same orientation, 0 when they are perpendicular and −1 when they go in opposite directions.

Long and Nelson [37] proposed 3 indices of dynamic interaction: 1) DI_d, dynamic interaction in displacement (average of all \(g_{t}^{\beta }\)); 2) DI_θ, dynamic interaction in direction (average of all f_t); and 3) DI, overall dynamic interaction, defined as the average of \(g_{t}^{\beta } \times f_{t}\) (Table 1). DI_d ranges from 0 to 1, DI_θ from -1 to 1, and DI from -1 (opposing movement) to 1 (cohesive movement). Figure 5 shows an example of the three indices.

In this section we used schematic, simple and contrasting case scenarios to evaluate the ability of the metrics to assess joint movement, in terms of proximity and coordination.

To build the case scenarios, we considered three levels of dyad proximity (high, medium and low); coordination was decomposed into two aspects, direction (same, independent and opposite) and speed (same or different). Eighteen case scenarios were thus built, with one example of dyad per scenario (Fig. 6; metrics in Additional file 5). The dyads for each case scenario were deliberately composed of a small number of fixes (∼10 simultaneous fixes, as in [37],) to facilitate interpretation of the metric values and the graphical representation of the arbitrarily constructed tracks (online access to tracks in github repository; see Availability of data and materials section). To assess the sensitivity of the metrics to changes in patterns of proximity and coordination, the case scenarios were grouped according to the categories in Table 2.

ID	Proximity	Coordination
		Direction	Speed
1	High	Same	Same
2	High	Same	Different
3	High	Independent	Same
4	High	Independent	Different
5	High	Opposite	Same
6	High	Opposite	Different
7	Medium	Same	Same
8	Medium	Same	Different
9	Medium	Independent	Same
10	Medium	Independent	Different
11	Medium	Opposite	Same
12	Medium	Opposite	Different
13	Low	Same	Same
14	Low	Same	Different
15	Low	Independent	Same
16	Low	Independent	Different
17	Low	Opposite	Same
18	Low	Opposite	Different

Due to the simplicity for its interpretation, Prox was defined as \(Prox_{\mathbbm {1},\delta }\). Three distance thresholds \(Prox_{\mathbbm {1},\delta }\) of 1, 2 and 3 distance units were used for Prox, HAI and CSEM, thus denoted for instance Prox₁, Prox₂ and Prox₃. For Cs, the original definition (Eq. 3) was used. jPPA, ϕ was arbitrarily fixed to 10. Regarding dynamic interaction, β was fixed to 1. The v variables for Pearson correlations (Table 1) were longitude (r_Lon), latitude (r_Lat) and speed (r_Speed). An average of correlations in longitude and latitude, denoted by r_Lonlat, was also computed. Boxplots of each metric were derived for each proximity and coordination category (Figs. 7, 8 and 9).

The values taken by Prox, jPPA, CSEM and, to a lesser degree, Cs, showed sensitivity to the level of proximity (Fig. 7). Conversely, no association was revealed between the proximity scenarios and the metrics based on correlation, dynamic interaction and reference area occupation.

Changes in direction were reflected in values taken by correlation metrics on location (r_Lonlat,r_Lon and r_Lat) and two dynamic interaction metrics, DI and DI_θ (Fig. 8). Cs took lower values in scenarios of opposite direction, but independent and same direction scenarios reflected no distinction for this metric. High correlation in speed was found for scenarios of opposite and same direction, while a large variability was found when direction was independent. r_speed showed differences when direction was independent between dyads, but no distinction was caught by the metric between same and opposite direction scenarios. The other metrics did not show distinguishable patterns related to changes in direction coordination.

Concerning coordination in speed, the most sensitive metric was DI_d, which measures similarity in the distances covered by individuals at simultaneous fixes (Fig. 9). r_Speed took a wide range of values when speed was not coordinated, while it was equal to 1 when perfectly coordinated. DI_d is more sensitive to changes in the values of speed (similar to step length because of the regular step units) than r_speed which characterizes variations in the same sense (correlation), rather than correspondence in values. HAI and L_ixnT showed slight differences in their ranges of values with changes in speed-coordination scenarios. When analysing combined categories of proximity and speed-coordination, and proximity and direction-coordination, less distinctive patterns were found, probably due to the higher number of categories, each containing fewer observations (Figure in Additional file 6).

Overall, Prox, jPPA, CSEM, r_Lonlat,r_Speed,DI_d,DI_θ and DI were highly sensitive to changes in patterns of either proximity or coordination. For proximity scenarios, the variance of some metrics for each category was also sensitive to the δ chosen; i.e. for larger δ, the variance of Prox and CSEM decreased in high proximity, while it increased for low proximity cases. This pattern does not hold for HAI, probably due to the strong dependence of this metric on the arbitrary choice of the reference area. Cs showed a slight sensitivity to changes in direction and proximity scenarios, although the values taken for each type of case scenario did not show a clear separation.

CSEM evaluates the similarity between the dynamical changes in movement patterns within a δ bandwidth, and, because of that, was expected to be more sensitive to changes in proximity than in coordination. It should be further assessed if using other variables for deriving CSEM (i.e. using [51] generic definition) could make it more sensitive to coordination than proximity. As with Prox, it is in the hands of the user to tune the threshold parameter. Because we were using locations as the analysed series (so the dynamical changes assessed were in fact changes in distance), we used exactly the same threshold values as for Prox. By contrast, correlations in location (r_Lon,r_Lat,r_Lonlat) did show sensitivity to changes in coordination, as expected. The same occurred with DI_θ and DI. Correlation in speed was sensitive to changes in both coordination components, showing high variance when there was no coordination (independent direction or speed). DI_d, on the other hand, was only sensitive to changes in speed. Because the time-step was regular, identical speed was equivalent to identical covered distance (at simultaneous fixes), which explained how in those scenarios DI_d was equal to 1. While DI behaved more similarly to DI_θ, its definition makes it impossible to separate the effects of coordination in displacement and in azimuth, which makes the interpretation of the metric more difficult than interpreting DI_d and DI_θ independently.

We also analysed the computational cost associated to these metrics. We simulated 50000 dyads with trajectories following a Brownian motion, each one composed of 100 fixes. Using a parallelization procedure, we found low CPU times for all metrics (<1 s) except jPPA (∼68 s). CPU time for jPPA and CSEM increased when we increased the number of fixes to 1000, to ∼161 and ∼94 s, respectively. It should be noted that for jPPA, the areas of intersection and union of the ellipses were approximated by grid cells, so for smaller cell sizes (i.e. more accurate jPPA estimation), the computational cost would increase. Researchers with long series of trajectories and a large amount of dyads should take this into consideration (results for the computational cost and more details on its calculation are in Additional file 7).

Although this review is directed at trajectory data (i.e. time series of locations that allow for movement path reconstruction) and the metrics presented here were defined for simultaneous fixes at regular time steps, technically speaking, some of these metrics could be computed only based on the identification of individuals simultaneously observed in a certain area (e.g. L_ixnT). These cases, which may be extremely sensitive to the spatial accuracy and the time intervals between observed fixes, are out of the scope of this review. For the case scenarios built to illustrate the metrics, we assumed that the granularity was correct, i.e. that the temporal and spatial resolution of the data were coherent in respect to the dyadic behavioural patterns under scope. Likewise, for practical uses of the metrics, researchers should 1) make sure that the spatiotemporal data that they are analysing allow reconstructing the movement paths of a dyad and 2) that the sampled (discretised) version of these paths are characterized by locations estimated with high precision, and that the time steps are small enough so that movement between two points could be assumed to be linear, so that the derivation of distances, speed and turning angles could be reliable. Further discussions on the importance of scale and granularity in the analysis of movement patterns can be found in [14, 30, 31].

We expected to obtain a binary classification of the metrics into proximity and coordination, based on the theoretical and case scenario evaluations. This was not so straightforward and we ended up instead with a 3-dimensional space representation (Fig. 10). Prox and CSEM are the most proximity-like indices. jPPA would be the third one due to its sensitivity to changes in proximity in the case scenario evaluation. Cs would be somewhere between Prox and direction coordination because it showed certain sensitivity to both HAI and L_ixnT are almost at the origin but slightly related to speed coordination. Theoretically, both metrics should account for proximity, since when two individuals are together in the same area, they are expected to be at a relative proximity; in practice, this was not reflected in sensitivity to proximity from HAI and L_ixnT. Still, HAI is represented in the graphic slightly above L_ixnT since its formulation specifically accounts for proximity in solitary use of the reference area. They are both graphically represented in association with the speed coordination axis because of the case scenario results which reflected that being in the same area only simultaneously requires some degree of synchrony. DI_d was the most sensitive metric to speed coordination, followed by r_Speed. DI_θ and r_Lonlat are the most strongly linked to direction coordination, seconded by DI, which is also related to speed coordination. A principal component analysis (PCA) using the values obtained for the case scenarios gave very similar results to those in Fig. 10 (Additional file 8), but this schematic representation is more complete because: 1) the theoretical and case-scenario assessment were both taken into account; 2) the PCA was performed without L_ixnT and HAI that had missing values for case scenarios with no common reference area (data imputation as in [25] was not appropriate for this case).

Figure 10 and Table 3 could be used as guidelines to choose the right metrics depending on the user’s case study. For instance, in an African lion joint-movement study [4], proximity was the focus of the study; in that case, the I_AB (Prox) metric was used. For similar studies several proximity-related metrics could be chosen; the choice would depend on the assumptions that the researcher is willing to make. In other cases, researchers may want to assess collective behaviour in tagged animals (e.g. birds or marine mammals) that do not remain proximal during their foraging/migration trips. Then, the collective behaviour component that could be evaluated would be coordination. Whether it is in direction or speed would depend on the researcher’s hypotheses. Coordination, or synchrony, has already been observed in some animal species such as northern elephant seals [17, e.g.] and bottlenecked sea turtles (e.g. [46]), among others. The use of the metrics presented here would allow a quantification of the pairwise behavioural patterns observed, a first step towards a quantitative analysis of the factors explaining those behaviours (e.g. physiological traits, personality or environmental conditions). The metrics presented here are applicable to any organism with tracking data (not necessarily georeferenced).

If the aim is to evaluate all three joint-movement dimensions, we advice to consider for each dimension at least one metric that is highly sensitive to it, rather than a metric that is weakly related to two or three. The complementarity of the metrics (i.e. multivariate approach) has not been studied here, and should be the focus of a future study.

The assessment of a ‘lagged-follower’ behaviour, where one individual would follow the other, was out of the scope of this work and should be addressed in the future. The study of this type of interactions is rather challenging, since the lag in the following behaviour is probably not static, and could vary between tracks and also within tracks. A few works use entropy-based measures similar to CSE (transfer entropy [54] or a causation entropy [40]), to measure how much the movement dynamics of an individual (called the source individual, or the leader) influences the transition probabilities in the movement dynamics of another individual [45, 61]. Some other works have focused on this type of interaction regarding it as a delay between trajectories and transforming the problem into one of similarity between trajectories, where one is delayed from the other [22, 27]. Metrics based on the Fréchet distance [1, 21] or the Edit distance [33] are common choices for measuring those similarities in computer science studies. In terms of computational cost, assessing following behaviour should be much more expensive than assessing joint movement.

This study focused on dyadic joint movement. The next step would be to identify metrics to characterize collective behaviour with more than two individuals. A pragmatical approach to investigate this more complex issue could be to identify, within large groups of individuals, the ones that move together for each given segment of trajectories (as dyads, triads or larger groups), and to study those dynamics. A similar procedure could then be used to spot following behaviour and leadership. Movement could be then regarded as spatio-temporal sequences of joint, following, hybrid and independence movement with one or more partners. Dhanjal-Adams et al. [15] present a Hidden Markov modelling approach to identify joint-movement states using metrics of direction and amplitude of flight synchronization in long-distance migratory birds (and assuming proximity between individuals). A similar approach could be used to identify more stages of collective behaviour, using several metrics as observed variables in the movement process.

Finally, a robust assessment of the different patterns of collective behaviour (e.g. proximal joint movement, coordination movement, follower movement) at multiple scales would provide realistic inputs for including group dynamic into movement models, which until now have relied on strong assumptions on collective behaviour in the few cases where it was taken into account [23, 29, 44, 47, 53], mostly due to the lack of understanding of collective motion.