Metrics for describing dyadic movement: a review

Abstract

In movement ecology, the few works that have taken collective behaviour into account are data-driven and rely on simplistic theoretical assumptions, relying in metrics that may or may not be measuring what is intended. In the present paper, we focus on pairwise joint-movement behaviour, where individuals move together during at least a segment of their path. We investigate the adequacy of twelve metrics introduced in previous works for assessing joint movement by analysing their theoretical properties and confronting them with contrasting case scenarios. Two criteria are taken into account for review of those metrics: 1) practical use, and 2) dependence on parameters and underlying assumptions. When analysing the similarities between the metrics as defined, we show how some of them can be expressed using general mathematical forms. In addition, we evaluate the ability of each metric to assess specific aspects of joint-movement behaviour: proximity (closeness in space-time) and coordination (synchrony) in direction and speed. We found that some metrics are better suited to assess proximity and others are more sensitive to coordination. To help readers choose metrics, we elaborate a graphical representation of the metrics in the coordination and proximity space based on our results, and give a few examples of proximity and coordination focus in different movement studies.

Introduction

Collective behaviour has been the object of study of many disciplines, such as behavioural ecology, psychology, sports, medicine, physics and computer sciences [7, 13, 19, 56, 57]. In multiple contexts, individuals – in a very wide sense of the word – adapt their behaviour as a function of their interaction with others. In movement ecology, where movement is regarded as an expression of behaviour [43], collective behaviour should be considered as a key element given that collective dynamics and individual movement are intricately intertwined [7]. Accordingly, mechanistic movement models should account for these dynamics. The vast majority of movement models neglect this aspect, with a few exceptions (e.g., [29, 44, 47, 53]). The consequence has been that the forms that these dynamics take in the few existing works rely on very simple theoretical assumptions.

Collective behaviour can be produced at large group scales (flocks, colonies, schools) but also at small group scales (triads, dyads). Regardless of the actual group scale, global patterns of collective behaviour originate from local interactions among neighbouring members [11], so analysing dyad interaction as a first step is a pertinent choice. Concerning dyadic interaction, here we focus on what we call ‘joint movement’, where two individuals move together during the total duration or a partial segment of their paths. Dyadic movement behaviour has been mostly studied in a data-driven approach, using several metrics to quantify it. In movement ecology, few works have applied and compared some of these metrics [38, 41]. However, their theoretical properties, and thus the similarities and differences in their construction and in what they actually assess, have not been thoroughly analysed yet.

This manuscript reviews a series of metrics used to assess pairwise joint-movement and proposes some modifications when appropriate (Table 1). Two criteria are taken into account for the review of these metrics: practical use and dependence on parameters; they are evaluated through both a theoretical (conceptual) as well as a practical approach. Metrics found in the literature essentially measured two aspects of joint movement: proximity and coordination. Proximity refers to closeness in space-time, as in how spatially close simultaneous fixes (individual locations recorded) are in a dyad; a point pattern perspective. The notion of proximity is thus subjective, since a judgement on proximity involves a threshold in distance whether local or global, or the definition of a reference zone (where encounters may be observed). Coordination, on the other hand, refers to synchrony in movement, which can be assessed through measures of similarity or correlation in movement patterns such as speed or direction. There might be a thin line between proximity and coordination, and some metrics may be associated with both at some degree, as we show through the description of their theoretical properties and the practical analysis of case scenarios.

Table 1 Metrics for measuring dyad joint movement

The manuscript is thus organized as follows. We first describe the criteria used to evaluate the metrics as indices of dyadic joint movement. We then present the different metrics and their theoretical properties with special attention to their dependence towards parameters. Next, we define case scenarios to evaluate the practical properties of the metrics. In the last section, we discuss the overall suitability of the metrics for assessing joint movement in ecology and give some practical guidelines for their use.

Evaluation criteria

We categorized the desirable properties of metrics for assessing dyadic joint movement into three criteria: practical use, considered the most important one; dependence on parameters; and computational cost:

• C1 Practical use [50, 52, 58]: 1) A metric is useful if it is interpretable and reflects a marked property of collective behaviour. 2) It should also be sensitive to changes in patterns of joint movement (e.g. higher values for high joint movement and lower values for independence in movement). 3) Being able to attain the theoretical range of values would also be important, as not doing so makes it harder to interpret empirical values. C1 is therefore a three dimensional criterion comprising interpretation, sensitivity and attainable range. Attainable range is covered in the theoretical properties section; we highlight the difficulties or implausibility to attain minimum and maximum values for the metrics when this is true. How to interpret each metric is also explained in this section; evidently, a metric without an attainable range is difficult to interpret. Sensitivity is addressed in the case-scenario section.

• C2 Dependence on parameters: A metric that depends on few parameters and hypotheses is more robust and generic than one that strongly relies on many parameters and hypotheses, since the former can produce more easily comparable results and interpretations. In addition, an ideal metric can be defined in such a way that the user can easily see how a change in the values of the parameters or in the components related to movement assumptions conditions the metric derivations and interpretations. In the next section, we describe the assumptions underlying each metric and the parameters needed to be fixed by the user. This description will allow distinguishing user-tractable parameter-dependent metrics from those that are not.

Definition and theoretical properties of the metrics

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)

Fig. 1

Example of Prox for δ=3 (left panel) and Cs (right panel). Circles and squares represent locations of two different individuals. Left panel: The numbers inside as well as the arrows represent the time sequence of both tracks. Grey lines correspond to the distances between simultaneous fixes; their values are shown. At the bottom: a dummy variable indicating if distances are below δ for each pair of simultaneous fixes, then the derived Prox and D_O (average of observed distances). Right panel: Grey lines represent the distances of all permuted fixes; D_E is their average

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.

Coefficient of Sociality (Cs)

The Coefficient of Sociality (Cs) [26] compares the mean (Euclidean) distance between simultaneous pairs of fixes (D_O) against the mean distance between all permutations of all fixes (D_E).

$$ Cs = \frac{D_{E} - D_{O}}{D_{E} + D_{O}} = 1 - 2 \frac{D_{O}}{D_{E} + D_{O}}, $$

(3)

where

$$D_{O} = \left(\sum\limits_{t=1}^{T} d_{t}^{A,B} \right)/T, $$

and

$$D_{E} = \left(\sum\limits_{t_{1}=1}^{T}\sum\limits_{t_{2}=1}^{T}d_{t_{1},t_{2}}^{A,B}\right)/T^{2}. $$

Kenward et al. [26] stated that Cs belongs to [−1,1], and it has been used as a symmetrical index since. Nevertheless, that is not true. Cs equals 1 if and only if D_O=0 and D_E≠0, which occurs only when the two individuals always share the exact same locations. However, Cs equals −1, if and only if D_E=0 and D_O≠0, which is impossible (it could asymptotically approach to 1 for very large series when DO approaches infinity). Cs equals 0 when D_O=D_E.

If all simultaneous fixes are very proximal but not in the same locations, Cs would approach 1 (how close to 1 would depend on the value of D_E as illustrated in the right hand side of Eq. 3). Moreover, only if D_E<D_O, Cs can take a negative value. For Cs to take a largely negative value, the difference in the numerator should be very large compared to the sum in the denominator; in Additional file 2 we show how implausible that situation is and how sensitive it is to the length of the series. The latter makes Cs from dyads of different length difficult to compare, because their real range of definition would differ. This fact is neither evoked in the work that introduced the metric [26] nor in the ones that evaluated this and other metrics [38, 41], despite the fact that in those works no value lower than −0.1 was obtained.

Indeed, [26] assumed that the permutation of all fixes is a way to represent locations of independent individuals. While this is questionable, some modified versions, as the one proposed by [62], use correlated random walks as null models and simulated independent trajectories under these models to replace D_E by a more realistic reference value. Thus, a generalized version of Cs would be:

$$ Cs = \frac{ D_{chance} - D_{O}}{D_{chance} + D_{O}}, $$

(4)

where D_chance is defined through a user-chosen movement model for independent trajectories.

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.

Fig. 2

Two examples of the derivation of L_ixnT and HAI. L_ixnT was computed using expected frequencies. HAI was computed with $K_{\delta }(t) = \mathbbm {1}{\{ d_{t}^{A,B} < 5 \}}$. Circles and squares represent locations of two different individuals. The numbers inside as well as the arrows represent the time sequence of both tracks. Grey lines correspond to the distances between simultaneous fixes; their values are shown. The dashed lines circle an arbitrary reference area

Coefficient of Interaction (L _ixn and L _ixn T)

Minta [42] proposed a Coefficient of Interaction (L_ixn) that assesses how simultaneous are the use and avoidance of a reference area S_AB by two individuals:

$$ L_{ixn} = \ln{\left(\frac{n_{AB}/p_{AB} + n_{00}/p_{00}}{n_{A0}/p_{A0} + n_{0B}/p_{0B} }\right)}, $$

(6)

where p_AB is the probability, under some reference null model, of finding A and B simultaneously in S_AB (the same interpretation as for n when a subscript is 0; see HAI subsection). Attraction between individuals would cause greater simultaneous use of S_AB than its solitary use, which would give positive values of L_ixn. Conversely, avoidance would translate into negative values of L_ixn, since use of S_AB would be mostly solitary. A logistic transformation of the metric (L_ixnT) produces values between 0 (avoidance) and 1 (attraction), making the interpretation easier:

$$ L_{ixn}T = logistic(L_{ixn})= \frac{1}{1+e^{-L_{ixn}}}. $$

(7)

Minta [42] proposed two different approaches for computing the associated probabilities conditionally to the fact that the reference area is known (see examples in Fig. 2 and Table in Additional file 3). In both cases, the probabilities are estimated under the assumptions of independence in movement among the individuals and of uniform utilization of the space. Indeed this latter assumption can be relaxed and p_AB can be derived from any kind of utilization distribution (see for instance [20] for the estimation of utilization distribution).

HAI and L_ixnT (thus L_ixn as well) rely heavily on a static reference area – either known or estimated – and on the probabilities of presence within this reference area. The static reference area could be defined, for instance, as the intersection of the respective home ranges of A and B. However, there are many approaches for estimating home ranges, each one relying on particular assumptions about the spatial behaviour of the studied populations [9]. Thus, S_AB is not a simple tuning parameter. The way it is defined may completely modify the output. If the reference area is equal to the whole area of movement of the two individuals, then both the numerator and the denominator in the logarithm are equal to infinity and L_ixnT cannot be derived. That problem could arise for extremely mobile individuals, such as tuna, turtles and seabirds [8], or fishing vessels [6], and avoiding it would require the computation of multiple dynamic reference areas. Therefore, L_ixnT may be better used for specific cases where the definition of the reference area relies on a deep knowledge of the spatial behaviour of the populations.

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.

Fig. 3

Example of the derivation of the joint potential path area (when ϕ=10). Circles and squares represent locations of two different individuals; the numbers inside represent the time sequence. The grey scales of the ellipses correspond to the time intervals used for their computation: from light grey for the [1,2] interval to dark grey for the [3,4] interval. The black regions with white dashed borders correspond to the potential meeting areas

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)

Fig. 4

Example of the derivation of CSE and CSEM when the compared features correspond to the positions of the individuals and δ=3. Circles and squares represent positions of two different individuals. The grey scales and arrows represent the time sequence of both tracks. Dotted lines represent the distances between simultaneous fixes; their values are shown. Values for all steps for CSEM computation are also shown

with the convention that max{∅}=0.

Correlations (r _V)

Pearson and Spearman correlations between variables such as longitude, latitude, distance, velocity, acceleration and turning angles from pairs of tracks, have been used as measures of synchrony in several studies (e.g. [16],). Correlations are easy to interpret. Pearson correlation coefficients (Table 1) assess linear correlations, while Spearman correlation coefficients based on ranks statistics capture any functional correlation. The correlation in a given V variable between dyads is denoted by r_V.

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.

Fig. 5

Example of a dyad for which correlations in longitude, latitude and an average of both (r_Lon,r_Lat and r_Lonlat, respectively), DI_d,DI_θ and DI are derived. Circles and squares represent locations of two different individuals; the numbers inside represent the time sequence. Displacement lengths and absolute angle values are also shown

Conclusions on the theoretical properties of the metrics

Practical use (C1): While each metric concerns a concrete aspect of joint-movement behaviour, some of them, such as Cs and DI, are harder to interpret. DI mixes up the coordination in displacement and direction. When DI is close to 1, it is certainly explained by high values in both components. When it is close to −1, it is an indication of overall high displacement coordination but in opposite directions. With values around zero, however, it is impossible to know if it is because of displacement or direction or both. For Cs, because obtaining values close to −1 is extremely rare, values around zero and, more particularly, slightly negative values are difficult to interpret. In addition, the maximum attainable value depends on the length of the series, which is likely to vary from dyad to dyad (Additional file 2).

Dependence on parameters (C2): Almost every metric depends on the ad hoc definition of a parameter or component, as summarized in Table 1. This is consistent with the fact that, since there is no consensus on the definition of behaviour [34], and much less on that of collective behaviour, its study depends heavily on the definition that the researcher gives to it. It should be noted that behind each choice of a parameter value, there is also an underlying assumption (e.g. that a distance below a δ value means proximity); the difference is that parameters can be tuned, and a variety of values can be easily tested. HAI and L_ixnT make a critical assumption of a static reference area, and its definition, which may be tricky for highly mobile individuals, is a key issue for the computation of both metrics. On the other hand, r_V and DI_θ are the only metrics that do not depend on parameter tuning or assumptions for its derivation; except for the assumptions of correlations being linear, or of linear movement between two successive positions when deriving directions, respectively.

Exploration of metrics through case scenarios

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.

Fig. 6

One example of dyad for each case scenario representing contrasting patterns of proximity and coordination (in direction and speed, C_Direction and C_Speed, respectively). Numbers correspond to scenario ID in Table 2. Solid lines represent the two trajectories, the solid points correspond to the start of the trajectories. The black dashed circumferences represent arbitrary reference areas; two circumferences correspond to an absence of a common reference area

Table 2 Case scenarios

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).

Fig. 7

Boxplots of each metric by category of proximity. Green, orange and purple correspond to case scenarios of high, medium and low proximity. For each category, the solid horizontal bar corresponds to the median, the lower and upper limit of the box correspond to the first and the third quartiles, while the solid vertical line joins the minimum to the maximum values. The green and purple boxplots are shifted to the left and right, respectively, to distinguish them better in case of overlap. X-axis: The metrics ranging from 0 to 1 are on the left (up to DI_d) while those ranging from -1 to 1 are on the right

Fig. 8

Boxplots of each metric by category of direction coordination. Green, orange and purple correspond to case scenarios of same, independent and opposite direction. For each category, the solid horizontal bar corresponds to the median, the lower and upper limit of the box correspond to the first and the third quartiles, while the solid vertical line joins the minimum to the maximum values. The green and purple boxplots are shifted to the left and right, respectively, to distinguish them better in case of overlap. X-axis: The metrics ranging from 0 to 1 are on the left (up to DI_d) while those ranging from -1 to 1 are on the right

Fig. 9

Boxplots of each metric by category of speed coordination. Green and orange correspond to case scenarios of same and different speed. For each category, the solid horizontal bar corresponds to the median, the lower and upper limit of the box correspond to the first and the third quartiles, while the solid vertical line joins the minimum to the maximum values. The green boxplots are shifted to the left to distinguish them better in case of overlap. X-axis: The metrics ranging from 0 to 1 are on the left (up to DI_d) while those ranging from -1 to 1 are on the right

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.

Synthesis of metric analysis

Table 3 summarizes the theoretical and case-scenario analyses. Most metrics reflected marked properties of dyadic joint movement, evidenced both theoretically and through the case scenario assessment. Exceptions were Cs, HAI and L_ixnT. Cs was sensitive to the null model for the distance expected by chance (D_chance; formula 4), it did not attain its whole range of definition, turned out to be asymmetric and dependent on the length of the series (Additional file 2), and was less sensitive than the other metrics to changes in patterns of joint movement. Perhaps a change in the null model for D_chance could improve Cs’s power to assess joint movement, though the new null model should be justified. HAI and L_ixnT, dependent on the reference area definition, were even less sensitive to changes in joint movement patterns. This supports our earlier statement that L_ixnT and HAI should only be used when a reference area exists and is known. Alternatively, Prox works as a simpler metric and is highly sensitive to changes in proximity. The only drawback of Prox is the need to choose a distance threshold parameter, eventually based on prior knowledge of the spatial dynamics of the population. Otherwise, a set of values can be tested, as shown here. jPPA presents the advantage of not requiring the knowledge of a reference area, but still relies on assumptions related to equal probability of presence in an ellipse, which strongly depends on a ϕ parameter whose tuning is not obvious.

Table 3 Evaluation of the two criteria for each metric

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).

Fig. 10

Representation of metrics in terms of their distance relative to proximity and coordination

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.

Further perspectives on collective behaviour

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.

Conclusions

The increasing availability of telemetry data for movement studies allow exploring patterns of collective movement. Here we reviewed metrics for assessing dyadic joint movement. We showed that some of the metrics were more suited for assessing proximity, others for coordination in direction or speed, and some others were not very sensitive to any of those aspects of joint movement. The results shown in this review offer guidelines to readers for choosing the metrics depending on which aspect of joint movement they would like to either describe or incorporate into movement models.

This study also contributes to highlighting the movement assumptions behind each metric as well as the parameters that need tuning. Users need to be able to decide whether these assumptions are realistic for their case studies, and to understand the consequences of their choice of parametrization. An accurate interpretation of movement patterns (here dyadic movement) relies on understanding the tools – in this case, metrics – used for obtaining those patterns.

Though the present work only concerns dyadic movement, further studies should concern the identification of larger groups moving together, where the size of the group would change in time, and metrics that would account for more than two individuals.