Unveiling the roles of temporal periodicity, the spatial environment and behavioural modes in terrestrial animal movement

Abstract

Background

Animal movement arises from complex interactions between animals and their heterogeneous environment. To better understand the movement process, it can be divided into behavioural, temporal and spatial components. Although methods exist to address those various components, it remains challenging to integrate them in a single movement analysis.

Methods

We present an analytic workflow that integrates the behavioural, temporal and spatial components of the movement process and their interactions, which also allows for the assessment of the relative importance of those components. We construct a daily cyclic covariate to represent temporally cyclic movement patterns, such as diel variation in activity, and combine the three components in a multi-modal Hidden Markov Model framework using existing methods and R functions. We compare the trends and statistical fits of models that include or exclude any of the behavioural, spatial and temporal components, and perform variance partitioning on the model predictions that included all components to assess their relative importance to the movement process, both in isolation and in interaction.

Results

We apply our workflow to a case study on the movements of plains zebra, blue wildebeest and eland antelope in a South African reserve. Behavioural modes impacted movement the most, followed by diel rhythms and then the spatial environment (viz. tree cover and terrain slope). Interactions between the components often explained more of the movement variation than the marginal effect of the spatial environment did on its own. Omitting components from the analysis led either to the inability to detect relationships between input and response variables, resulting in overgeneralisations when drawing conclusions about the movement process, or to detections of questionable relationships that appeared to be spurious.

Conclusions

Our analytic workflow can be used to integrate the behavioural, temporal and spatial components of the movement process and quantify their relative contributions, thereby preventing incomplete or overly generic ecological interpretations. We demonstrate that understanding the drivers of animal movement, and ultimately the ecological phenomena that emerge from it, critically depends on considering the various components of the movement process, and especially the interactions between them.

Background

Movement is ubiquitous to animal life and central to emergent ecological phenomena such as territoriality and home range formation [1, 2], predator-prey dynamics [3, 4], distribution patterns [5] and population dynamics [6, 7]. To understand such phenomena it is thus key to study what drives the movement of animals through their environment [7, 8]. However, this understanding is hampered by the complex relationships between animal movement and a multitude of behavioural and spatial-temporal environmental drivers, since a specific state of a driver can give rise to different movement patterns depending on the state of other drivers. For example, animals might move slowly through thick vegetation when foraging, or fast if they are translocating from one preferred high-visibility habitat to another. Temporal variation in the response of animals to their environment in diel, lunar or seasonal patterns [e.g. 9], complicates the analysis of animal movement drivers even further. In order to structure the study of animal movement and advance the interpretation of movement patterns, a movement ecology framework has been put forward that divides the movement process into multiple components: an internal motive regarding why to move and the temporally and spatially varying environment influencing when and where to move respectively [10].

First, ‘why move’ relates to the internal motives of movement (e.g. reproduction, resource acquisition and threat avoidance), which are often manifested in different behaviours (e.g. mating, foraging and fleeing) that produce distinct movement modes (i.e., movement trajectory geometries that can be quantified by metrics like speed and turning angle distributions) [10,11,12]. Distinguishing different movement modes in animal trajectories has yielded insights in different foraging strategies within individuals [13], predator-prey interactions [14] and effects of anthropogenic disturbance on behaviour [15]. Internal motives and their emergent movement modes should thus be acknowledged to further the interpretation of animal movement, specifically regarding its iteraction with the environment [11,12,13, 16].

Second, ‘when to move’ deals with temporal variation in animal movement, which often exhibits cyclic patterns. For example, seasonal migratory movements [e.g. 17, 18] have long been recognized, as well as the influence of lunar rhythms [e.g. 19, 20]. Moreover, most animals show strong diel activity patterns [21,22,23,24,25], which can arise directly due to variation in light and thermal conditions [26] or indirectly due to predation risk and human activity [22, 27,28,29]. Temporal periodicity is thus a key characteristic of animal movement [9]. Several methods of varying complexity have been developed to address temporal variation in movement [24, 25, 30,31,32,33]. Basis functions can be used to accommodate temporal heterogeneity in continuous-time movement models [32], for example by modelling processes related to energy discharge and recharge with a temporally explicit recharge function [33], or by augmenting spatial movement models with a generalised additive model that describes temporal movement heterogeneity (such as seasonal and daily patterns) [25]. In discrete-time descriptions of movement, resource selection can be weighted with time-specific distributions to reveal time-varying effects of resource covariates on selection [31], or transition probabilities between movement modes can be modelled with temporal heterogeneity [30].

Third, ‘where to move’ relates the movement of animals to the characteristics of their spatial surroundings. For example, topography influences movement through energetic costs associated with moving through sloping or rugged terrain [34,35,36,37]. Daily cycles can interact strongly with the spatial surrounding, for example in animals that select areas with high vegetation cover during daytime and more open areas during the night to decrease predation risk [38, 39]. These examples illustrate that relationships between spatial environmental characteristics and movement patterns are conditional on the context of internal motive and time, as this context determines the costs and benefits that the spatial environmental characteristics represent [e.g. 40].

Internal motives, temporal periodicity and spatial environmental heterogeneity are thus three intricately linked components of the animal movement process. However, animal movement research has often been hampered by the focus on univariate analyses (Fig. 1) [41], while ideally the relationships between movement and each of these three components are considered within each other’s context [10]. Although the multi-component movement ecology framework has proven helpful in elucidating drivers of animal movement [42] and a plethora of methods exists to study the various components, challenges remain regarding how to (1) include the components in a single coherent modelling workflow, (2) quantify the relative contributions of the components in shaping animal movement, and (3) assess the influence of interactions between the components on the quantitative and qualitative inferences made about the movement process.

Fig. 1

Results of a systematic literature search (Scopus, 06-06-2024; see Supporting Information), demonstrating the numbers of studies that focused on the different components of the animal movement process (modes, time and space) separately or combined

Here, we demonstrate an integrative analytic workflow, drawing from and combining existing methods, to evaluate the role of movement modes, time (temporal rhythms) and space (environmental heterogeneity) in shaping animal movement. We apply this workflow to a case study of large savanna ungulates. We quantify the relative contribution of these three components to the movement process in isolation as well as in interaction, and we assess how ecological interpretations of movement change when omitting components and interactions from the analysis. We use fine-scale (10-minute resolution) GPS tracking data of eland antelope (Taurotragus oryx), blue wildebeest (Connochaetes taurinus) and plains zebra (Equus quagga) in a South-African game reserve, and fit Hidden Markov Models (HMMs) with input from tree cover, terrain slope and a constructed daily cyclic covariate (which captures the diel effect on the movement behaviour of these crepuscular species). Since HMMs can describe the movement process as a multi-modal random walk with input from spatial-temporal environmental covariates [43], they provide an ideal method to quantify the effects of modes (as a proxy for internal motives) and temporal and spatial environmental heterogeneity on movement. We include effects of the time of day and the spatial environment in both the mode-specific movement properties speed and turning angle, as well as the switching probabilities between modes. Besides demonstrating our analytic workflow, we use our case study of savanna ungulates to demonstrate that inferences about the behavioural ecology of animals are made best when explicitly modelling the movement process with these three components in interaction.

Methods

Animals and study area

From June to August 2017, 34 eland antelope, 34 blue wildebeest and 35 plains zebra were captured in Welgevonden Game Reserve, Limpopo, South Africa (24°13’S, 27°54’E), and equipped with neck collars containing GPS and accelerometer sensors. For more details about the collaring process and sensors, see [44]. The animals were released in a fenced study area of 1200 hectares at the northern edge of the reserve (Fig. 2). As a part of the Waterberg mountain massif and Bushveld ecoregion [45], the study area is characterised by nutrient poor sandy soils in the lower northern part and rocky, flat hilltop plateaus dissected by ravines extending towards the south. Woody vegetation varies from mixed to open dry deciduous woodland, with the rocky plateaus being the most open.

Fig. 2

Overview of the study area with (a) terrain slope (radians) and (b) fractional tree cover. Numbers at the isoclines denote elevation from sea level (m)

Animal movement data

The GPS sensors collected positions with an interval ranging between 15 and 2 min, depending on the amount of body activity as determined by the accelerometer [44]. The position data were corrected, smoothed and modelled to regular 10 min resolution trajectories as described in [44]. Due to sensor failures, we selected from each species only the ten individuals with the most recorded data in the period from September to December 2017 (the first four months of the study period, during which the sensors collected most data). Our subset contained a total of 357,815 animal positions across 30 individuals of the three species, comprising 481 bouts (i.e. an uninterrupted and regular sequence of positions of a single individual) with a median bout length of 496 positions (3.4 days).

Environmental data

A digital elevation model of the study area with resolution 2 × 2 m was available, from which terrain slope (hereafter referred to as slope) was calculated (Fig. 2a). We used aerial images of the study area from August 2013 to create an orthomosaic with a 0.25 × 0.25 m resolution, from which pixels were classified into trees, grass, bare ground and other/built up area using the Semi-Automatic Classification Plugin for QGIS [46]. After having applied an isotropic Gaussian kernel smoother with a 10 m standard deviation to the classified pixels, we derived fractional tree cover (hereafter referred to as tree cover) and fractional grass cover maps of the study area at a 2 × 2 m resolution (Fig. 2b). Tree cover was strongly negatively correlated with grass cover (r = 0.95), so we used only tree cover in our analyses and considered it a measure for the inverse of grass cover as well. We extracted tree cover and slope values for each animal position.

Movement model structure

We constructed a Hidden Markov Model (HMM) for each species separately to examine how their movement was shaped by movement modes, diel patterns and spatial environmental heterogeneity. We described the movement process as a random walk consisting of steps between subsequent animal positions with lengths s and turning angles φ. Since intervals between steps were constant at ten minutes, we refer to s as speed throughout. We fitted gamma (Γ) and von Mises (vM) distributions to s and φ, respectively, which were unique to each discrete movement mode z. Thus, z was defined as a unique combination of three movement parameters: the mean μ and standard deviation σ of the speed (converted from the shape and scale parameter of the Γ distribution for easier interpretation) and the vM concentration parameter κ of the turning angle (the vM mean was fixed at 0 since there was no reason to assume a movement bias to either the right or left direction). μ, σ and κ were in turn modelled as functions of several covariates (see section Movement model covariates) using a log-link. Mode switching between the steps was governed by a switching matrix M, where the off-diagonal elements are logit-linear functions of covariates. We considered a three-mode, two-mode and a ‘single-mode’ HMM (i.e. 3 × 3, 2 × 2 and 1 × 1 dimensions for M respectively), where the single-mode model is analogous to a movement model that does not include the mode component. Thus, the HMM had the general form:

$$M = \left[ {\begin{array}{*{20}{c}}{{\gamma _{1 \to 1}} = 1 - \sum _{z = 2}^n{\gamma _{1 \to {\rm{z}}}}}& \cdots &{{\gamma _{1 \to {\rm{n}}}} = {f_a}\left( {{\theta _a},{\xi _a}} \right)}\\\vdots & \ddots & \vdots \\{{\gamma _{{\rm{n}} \to 1}} = {f_b}\left( {{\theta _b},{\xi _b}} \right)}& \cdots &{{\gamma _{{\rm{n}} \to {\rm{n}}}} = 1 - \sum _{z = 1}^{n - 1}{\gamma _{{\rm{n}} \to {\rm{z}}}}}\end{array}} \right]$$

$${s_{z = i}} \sim \Gamma \left( {\mu = {f_c}\left( {{\theta _{c,i}}{\xi _c}} \right),\sigma = {f_d}\left( {{\theta _{d,i}},{\xi _d}} \right)} \right)$$

$${\phi _{z = i}} \sim vM\left( {\mu = 0,\kappa = {f_e}\left( {{\theta _{e,i}},{\xi _e}} \right)} \right)$$

Where γ_1→2 is the probability of switching from mode z = i to z = j between consecutive steps, f_n are functions for mode switching probabilities and Γ and vM parameters, θ_n are sets of to-be-fitted coefficients of function f_n (and for mode z = i in case of θ_{n, i}), and ξ_n are sets of covariates (see next section).

Movement model covariates

We modelled the parameters μ, σ and κ and the mode switching probabilities γ_1→2 each as a function of three covariates. Two of those were the spatial environmental covariates tree cover and slope. The third covariate was a constructed daily cyclic covariate (DCC) that describes the diel patterns in μ, σ and κ (driven by e.g. ambient temperature, solar position, perceived predation risk and internal diel rhythm). To determine those diel patterns, we first divided all the animal positions in 10-minute time bins based on the time of day that they were taken (i.e., there were a total of 144 time bins of 10 min, totalling 24 h). We then fitted Γ distributions (parameters μ and σ) to the movement speeds and vM distributions (parameter κ) to the turning angles of the positions within a centred moving window of three time bins wide. This gave us empirical diel patterns in parameters μ, σ and κ (Fig. S1). We rescaled and averaged those empirical diel patterns to one aggregate diel pattern per species (Fig. 3), because the separate patterns in μ, σ and κ were strongly correlated (Fig. S2). Finally, we fitted to that aggregate diel pattern a mixture of a cosine and four wrapped normal distributions with 24 h periods, separately per species, with the wrapped normal distributions representing two peaks and two lows per day that were fixed at the peaks and lows of the aggregate diel patterns to capture the crepuscular movement activity of these species (Fig. 3). We defined this mixture distribution as our daily cyclic covariate (DCC). Thus, DCC had the following form:

$$DCC = \alpha + \beta {\rm{cos}}\left( {{{\rm{t}}_{{\rm{day}}}}{\rm{ + }}\gamma } \right){\rm{ + }}\sum _{{\rm{n = 1}}}^{\rm{4}}{\delta _{\rm{n}}}{{\rm{f}}_{{\rm{WN}}}}\left( {{\mu _{\rm{n}}},{\sigma _{\rm{n}}}} \right)$$

Where t_day is the time of day (range 0–2π), f_WN is the probability density function of the wrapped normal distribution with a 24 h period (describing 1 out of 4 peaks or lows per day in the aggregate diel pattern), and δ_n are scalars for those distributions. Note that μ_n and σ_n here denote the mean and standard deviation of the peaks and lows in the aggregate diel pattern, and are not related to the mean and standard deviation of movement steps. As a result, each animal position had a corresponding DCC value, which depended on the time of day of that position and the species (Fig. 3).

By its design, DCC captures both the purely temporal (space-invariant) effect of time of day on animal movement and the combined spatial-temporal effect of diel variation in habitat selection, but it excludes the purely spatial (time-invariant) environmental effects. Using DCC in our analytic workflow therefore enables quantifying which fraction of the animals’ movement is tied directly or indirectly to time of day. If we were to quantify the influence of all three components in one HMM modelling step, without using the constructed DCC, that would not allow us to distinguish between variation caused by purely temporal (space-invariant) effects, purely spatial (time-invariant) effects and their combined spatial-temporal effects.

We standardized tree cover, slope and DCC to zero-mean and unit variance, and included all three variables and their two-way interactions as covariates of μ, σ and κ in the HMMs. We included tree cover, slope and DCC without interactions as covariates of the mode switching probabilities γ_1→2.

Fig. 3

Aggregate diel patterns averaged from μ, σ and κ for eland, wildebeest and zebra (solid line). We fitted a daily cyclic covariate (DCC, dashed line) to those aggregate diel patterns, which summarizes the diel activity patterns per species as expressed in the movement properties μ, σ and κ. DCC is a mixture distribution of a cosine and four wrapped normal distributions with 24 h periods

Model fitting

We fitted a series of HMMs per species separately. The HMM series consisted of models containing any combination of the three components of the movement process (modes, time and space), ranging from a null model that contained none, to a full model that contained all three in interaction (Table 1). The models were fit by numerical maximisation of the log-likelihood function, using the fitHMM routine from the R package momentuHMM [47]. We fitted the models in order of increasing complexity, every step adding covariates and/or interactions between covariates to the distribution parameters and switching probabilities. The full model included all three movement components modes, time and space, namely: DCC, tree cover and slope as covariates to the mode switching probabilities γ_1→2, and these three variables and their two-way interactions as covariates to μ, σ and κ. This way, we could determine whether the covariates influenced movement within the separate movement modes, or through affecting the probability of switching to another mode, or both. The fitted parameter values from each fit were used as starting values for the next; starting values for other parameters were set to 0 (note that all covariates had been scaled to mean 0 and standard deviation 1). For each iteration, four sets of perturbations of those starting values β were sampled from $\sim {\rm{N}}(\beta ,\sigma = 0.25 \cdot \beta )$. The model fit with the highest likelihood out of all starting sets was the optimal fit of that iteration, and it was used to initiate the model fit of the next iteration. However, in practice all iteration sets of fitted models converged to almost identical coefficient values and model likelihoods. We used the Akaike information criterion [AIC; 48] to assess model parsimony throughout the workflow. The model structure and fitting sequence were identical for all three species.

Variance partitioning of movement process components

Based on the fitted HMM with all predictors included (Model 3 in Table 1), we partitioned the explained variance of the modelled parameters μ, σ and κ for each species into contributions from the three components of the movement process: movement mode, time (the DCC covariate) and space (the covariates tree cover and slope). To quantify what proportion of the total variance in the predictions of μ, σ and κ was accounted for by each of the three components, we used linear regressions with the predicted values for μ, σ and κ as response variables (on the log-scale). The predicted values for μ, σ and κ were conditional on the most probable movement mode for each step along the movement path, which was determined using Viterbi decoding [49]. As predictors, we alternately omitted or included the movement mode component (i.e. the Viterbi path with the most likely hidden state sequence, included as predictor via a dummy variable), the time component (the DCC covariate), the space component (the tree cover and slope covariates) or their interactions as predictors and summarised the regression fit using R² values. By doing so, we partitioned the variance of the fitted model parameters μ, σ and κ into the independent model components (i.e. mode, time and space), as well as their combined effects. All analyses were performed in R 3.6.1 [50].

Results

Hidden Markov models

We fitted Hidden Markov Models (HMMs) that included combinations of the three components of the movement process: modes, time and space (Table 1). When including the modes component (Models 1a, 2a, 2b and 3), we explored models with two and three movement modes. The full three-mode models achieved a higher model fit than the full two-mode models, albeit with a relatively small difference (ΔAIC = 32,093; 29,170 and 27,381 for eland, wildebeest and zebra respectively). However, AIC as a model selection criterium generally favours the inclusion of more modes, so the number of included modes should rather be informed by qualitative inspection of interpretability and biological meaningfulness of the modes’ properties [51]. Both the two- and three-mode models detected one faster, more directed mode and one slower, more tortuous mode (Fig. S3). We labelled these as ‘transit’ and ‘encamped’ modes, considering the labels to be reasonable abstractions while retaining spatially descriptive meaning and avoiding to project too specific behavioural annotations on the modes. The three-mode HMMs detected an additional slow and directional mode (Fig. S3). This third mode responded relatively weak to the temporal and spatial environmental covariates, while the response displayed by the other two modes was qualitatively similar for both the two- and three-mode HMMs (Figs. S4-S6). Therefore, we do not further consider the three-mode HMMs, and hereafter refer to the two-mode HMMs when discussing models that include movement modes.

The AIC values of all HMMs are presented in Table 1. The inclusion of movement modes led to much higher model fits, with the models including modes but not time and space (Model 1a) outperforming the models including time and space but not modes (Model 2c; ΔAIC = 63,185; 39,136 and 41,029 for eland, wildebeest and zebra respectively). Furthermore, the models including time but not space were superior to those including space but not time (Model 1c versus 1b; ΔAIC = 9,797; 14,153 and 5,507 for eland, wildebeest and zebra respectively; Model 2b versus 2a; ΔAIC = 9,335; 7,301 and 3,087 for eland, wildebeest and zebra respectively), indicating that DCC is a dominant covariate accounting for substantial variation, which is not accounted for by tree cover and slope. For each species the full model, which includes modes, time and space (Model 3), was the most parsimonious model according to AIC value.

Table 1 Overview of the fitted HMMs, each including a combination of the predictors modes, time and space

Variance partitioning of movement process components

Partitioning the variance in the modelled parameters μ, σ and κ for each species between the modes, time and space components showed that modes explained most of the variation (see Fig. S7). The contributions of the model components to the within-mode variation in model parameters for the three species are shown in Fig. 4. For most combinations of species, parameter and modes, the temporal component explained most of the variance. However, the interactive effects of spatial-temporal predictors contributed notably to the variance in κ, especially for zebra (Fig. 4). By themselves, the spatial predictors did not explain much variation in movement parameters except for κ in zebra, and σ in eland (encamped mode) and zebra (transit mode; Fig. 4).

Fig. 4

Variance partitioning, separately per movement parameter, species and movement mode (enc.= encamped, trans = transit). For variance partitioning that includes modes as a partition, see Fig. S7

Different model descriptions of the movement process

Incrementally expanding the predictor set according to Table 1 (from Model 0 to 3) yielded various interpretations of the movement process (Figs. S8-S14). Model 0, which included neither space, time nor mode as predictors, showed that wildebeest moved on average slowest and most tortuous, with eland and zebra moving equally directed and zebra fastest (Fig. S8). When we included movement modes (Model 1a), zebra showed the same speed as eland in the encamped mode and wildebeest moved faster than eland in the transit mode, instead of being slower overall (Fig. S9). Model 1b, which included only time as a predictor, showed large fluctuations in speed and tortuosity throughout the day, with both movement parameters peaking around sunrise and sunset and showing lows before sunrise and after sunset (Fig. S10). Model 1c, including only spatial predictors, showed that speed decreased with tree cover for eland, but increased for the other species (Fig. S11). Furthermore, speed decreased with slope irrespective of tree cover. Tortuosity decreased with tree cover but increased very slighty with slope. When we included both mode and time (Model 2a), we found the same temporal patterns in speed and tortuosity as the model including only time, but also indications that daily peaks in speed were much higher for animals in the transit mode than those in the encamped mode (Fig. S12). Model 2b, including both modes and space as predictors, showcased that often only in one of the two modes did animals respond substantially to spatial environmental heterogeneity (Fig. S13). For example, only in the transit modes did animals strongly decrease (eland) or increase (wildebeest and zebra) their speed. Zebra did not respond to slope, whereas the model including only spatial predictors suggested it slowed down on sloping terrain. Including space and time as predictors (Model 2c) reveiled that tree cover, slope and time of day interacted in their effect on speed and tortuosity, with their effects being most pronounced in the evening (Fig. S14). Model 3, including space, time and modes all as predictors, yielded the most complete description of movement, with aspects that were not apparent from all previous models (Fig. 5a, b). The movement properties displayed a strong diel pattern, with the effects of tree cover and slope on speed being largest in the evening when animals generally moved faster. Wildebeest and zebra in their transit mode moved faster with increasing tree cover whereas eland moved slower, but only so on steep terrain. In the encamped mode, tree cover had almost no effect on speed, but a strong effect on tortuosity. The animals also generally moved more tortuously with higher tree cover, even more so on sloping terrain. The model indicated that the spatial or temporal environment did not impact movement in all modes, and that the interaction effect between slope and tree cover on movement was also mode-specific.

Fig. 5a

Model predictions of movement speed μ as a result of tree cover and slope in the morning (04:00, top row), at noon (12:00, middle row) and in the evening (18:00, bottom row) for eland, wildebeest and zebra. The times were chosen to capture a dip, baseline and peak in DCC (see Fig. 3). Blue and red indicate the transit and encamped mode respectively. Colour mapping within blue and red indicates a progression of slope from 0 (light) till 0.15 (the 90th percentile of slope values in the complete dataset; dark). Small ticks at the inside of tree cover axes mark the 10th, 20th … 90th percentile values of tree cover in the respective data sets. Note that the y-axis scales differ between species, but not between the times of day. Confidence intervals have been omitted for visual clarity, but these were generally small

Fig. 5b

Idem to Fig. 5a, instead with model predictions of turning angle concentration on the y-axis. Concentration parameter κ (range 0-∞) has been converted to ρ (range 0–1) in this figure for visualisation purposes

Discussion

We developed an integrative analytic workflow to quantify and understand the roles of movement modes (‘why move’), temporal rhythms (‘when to move’) and spatial environmental heterogeneity (‘where to move’) in shaping animal movement. In this workflow, multi-modal (why) Hidden Markov Models (HMMs) are fitted on animal movement data, with spatial environmental covariates (where) and our constructed daily cyclic covariate (DCC; when, which aggregates the empirical diel patterns in the animals’ movement parameters) influencing both the movement parameters and the mode-switching probabilities. Finally, by using a variance partitioning approach on the modelled parameters, we quantified the relative importance of the why-, when- and where-components and their interactions in shaping animal movement.

We demonstrated the usefulness of our workflow by applying it on a case study of African ungulates. The variance partitioning showed that movement modes impacted movement the most, followed by diel rhythms and the spatial environment (i.e., tree cover and terrain slope) being the least important. Furthermore, the interactions between these components were important contributors as well, often explaining more of the movement variation than the marginal effect of the spatial environment did. This contrasts recent movement ecology research, where most studies seem to focus on the spatial environment without considering movement modes and temporal rhythms (Fig. 1). That is not to say that the spatial environment is the least important driver of animal movement across all study systems, or that research should focus solely on the most important movement driver of a system. However, when studying a specific component of interest of the movement process (e.g. a spatial environmental variable), that component should still be considered properly in context of the full (why, when, where) movement process [10], to avoid drawing misleading conclusions.

To showcase that misleading conclusions can follow from models in which one or more components of the movement process are ignored (‘omitted-variable bias’), we applied the same analytic pipeline to all combinations of the movement components in our case study. We identified two types of misleading conclusions that followed from models that only considered a subset of the movement process components: (1) failure to detect relationships and (2) detecting questionable relationships. First, the failure to detect relationships leads to overgeneralisations when drawing conclusions about the movement process (resembling a type II error). For example, wildebeest appeared to move slowest of the three species when not considering any movement process component in the model, but when movement modes were considered, the transit mode of wildebeest turned out to be faster than that of eland (Figs. S8 and S9 respectively). While this does not equate an erroneous conclusion (i.e., on average wildebeest did indeed move slowest), it does inhibit uncovering ecologically relevant complexities in the movement process [52]. Second, various questionable (likely spurious) relationships between input and response variables were revealed in models that included only a subset of all movement components that were present in the most parsimonious model (as determined by model AIC value; resembling a type I error). For example, zebra appeared to slow down substantially with increasing slope in the model with only spatial components, but with the addition of movement modes this effect disappeared completely for both modes (Figs. S11 and S13 respectively). Apparently, mode occupancy was a confounding variable that wrongly implied a causal relationship between slope and speed (i.e. zebras were more likely to be in their encamped mode in steep areas than they were in level areas; Fig. S15). Finally, besides these two types of misleading conclusions, the omission of a movement process component that is uncorrelated to other components does not lead to overgeneralisations or questionable relationships, but can nonetheless still lead to simply ‘incomplete’ conclusions. For example, a model with both movement modes and temporal predictors (Fig. S12) displayed only additive effects when compared to the models that solely contained movement modes (Fig. S9) or temporal predictors (Fig. S10).

An important consideration when using our analytic workflow, is that we developed DCC to capture both the direct, space-invariant effects of time of day on the animals’ movement properties, as well as the indirect effects of time of day through diel variation in habitat selection (given that varying environmental conditions between habitats also influence movement properties). For example, crepuscular animals often move faster during early morning because they forage before the midday heat sets in, but on top of that they also move faster at that time of day since they select for open foraging habitat which provides less movement resistance. Therefore, DCC includes both a purely temporal and a spatial-temporal component, which allows to separate those from purely spatial, time-invariant effects on the movement process (as these are captured by the spatial environmental, time-invariant covariates tree cover and slope). The influence of the temporal and spatial-temporal components together on movement was large compared to that of purely spatial environmental heterogeneity. Future research could focus on further disentangling these temporal and spatial-temporal components. For example, in the construction of the daily cyclic covariate, instead of solely averaging the movement parameters μ, σ and κ within each time bin, these movement parameters could be modelled based on time bin-specific spatial environmental covariates, after which the model predictions for e.g. the mean of the spatial environmental covariates could be used as the DCC value of that time bin. In this way, the purely temporal variation in movement could be separated from the temporally varying spatial influence on movement. Moreover, we see opportunities to build upon our workflow with datasets that are larger in spatial-temporal scale. Namely, seasonality (e.g. through climatic variation) should come into play on top of diel rhythms as a temporal covariate when the movement data spans more months, and extra spatial environmental covariates such as vegetation quality or water availability should be added when animals roam over larger or more heterogenous areas.

Conclusions

Our integrative analytic workflow can be used to analyse the effect and relative importance of all components (why, when and where to move) of the movement process. Compared to analyses that focus on a subset of these components, our workflow prevents the drawing of misleading conclusions through overgeneralisations and spurious correlations. Understanding the drivers of animal movement, and ultimately of ecological phenomena that emerge from it, critically depends on considering the various components of the movement process in concert, and especially the interactions between them.