Our analytic framework quantifies the contribution of an environmental variable to animal movement by utilizing the multivariate richness of movement data. Instead of building a model to predict a simplified animal movement descriptor from a set of environmental predictors, i.e., the route of causal inference, we turn this around and build a model to predict an environmental variable from a large number of animal movement variables. By using animal movement variables, the model of this framework predicts a perceived environmental variable by the animal [6,20]. Although predicting an environmental variable from movement data is the goal of the model, it is an intermediate step of the framework in order to quantify environmental contribution to animal movement. In this framework it is key to use as many informative movement variables as possible, which could be meaningful human-constructed ecological (e.g., variables related to multiple classified animal activities), mathematical and/or physical variables, or abstract variables from an automated (deep learning) feature extraction algorithm. When effort is made to extract as many informative variables as possible from the animal movement data, chances are maximized that most of the variation of the environmental variable under scrutiny that is present in the data is captured. Furthermore, instead of creating the model as the end product during the analysis, the environmental variable should be predicted on a separate test dataset as the final step of the analysis. This follows from a data-driven and machine learning philosophy, in which complex multivariate models can be built that are not overfitted and therefore generalize better to new datasets. When distinguishing the train and test dataset, the test set used in the prediction of the environmental variable needs to be from a different temporal range than the train set that is used in the model building phase, due to autocorrelation in animal movement data that can otherwise cause the model to overfit [21]. The range of values in the test set of environmental variables (whether or not these are under scrutiny) should be comparable to the range of values in the train set, to prevent incorrect extrapolation. After generating model predictions on the test set, the coefficient of determination (R2) quantifies the fit of this predicted environmental variable from animal movement data to the measured environmental variable on a known scale and can thus be considered a metric on how much of the variation in the environmental variable contributed to animal movement in a multivariate fashion (see Additional file 1) [22]. This is in contrast to approaches where the correlation between a set of environmental variables and one or several simplified animal movement variables are studied [2,3,9], because here the environmental dependency of (instead of the contribution to) movement is often quantified, i.e., the variation in animal movement that is dependent on environmental variables. The measure of fit of the null model (i.e., no environmental contribution) should be chosen depending on the algorithm that is used, which is R2 = 0 for algorithms that are able to always predict the mean of the response variable (e.g., Support Vector Regression and Random Forest Regression), even when the input variables are white noise. The measure of fit of this null model will then form the baseline value for which there is a 0% environmental contribution and an R2 of 1 can always be interpreted as 100% environmental contribution. Of course R2 should only be used as the measure of fit when modelling a continuous environmental response variable. With a discrete environmental variable, a classification approach should be undertaken, which is outside the scope of this study. However, to compare the contribution of different environmental variables with each other fairly, the same measure of fit should be used.
In order to demonstrate the usefulness of the proposed analytic framework, we applied this framework to a case study about the contribution of resource availability (here grass biomass), time since milking, and wind speed to the movement of eight dairy cows in a pasture (Fig. 1). When animals are facing resource depletion, movement characteristics (through the landscape and of body parts), and emergent patterns like group (herd) characteristics, and time allocated to specific activities (e.g., foraging) often change, because animals need to invest more time and/or energy in searching for and acquiring resources [23,24]. Cows in a pasture are a good model for such a case study, because this provides a relatively homogenous foraging arena. Time since milking is another variable that could substantially influence the movement of dairy cows, because it has been shown that the lactation stage of cows (a variable that is intuitively linked to time since milking regarding its effect on cow behaviour) influences the relative distribution of their activity patterns [25,26]. Wind speed provided a good test case for our framework, because it was moderately correlated (r = 0.37) with grass biomass. We expected this correlation to be spurious and the effect of wind speed on cow movement to be negligible, because conditions were mild during the experiment (0-9 m s-1).
The exact methodological approach that we describe for this case study is one possible implementation of our proposed analytical framework (Fig. 1). However, there are numerous possible implementations of this framework for other studies, which may be influenced by the problem statement, experimental setup, animal movement sensors, environmental data types, data quantity, etc. However, the property that all implementations should have in common is that environmental contribution to animal movement is quantified by predicting environmental variables from movement descriptors in a data-driven (viz., machine learning) approach, which uses the coefficient of determination as a measure to quantify this contribution.
Data collection
For this case study eight adult female Holstein-Friesian dairy cows were kept in controlled pastures that were small enough so that foraging lead to resource depletion over the course of several days. The experiment ran from 25 April until 11 May 2017. During the experiment, the cows’ movements were recorded continuously with e-Track neck collars (Noldus InnovationWorks, Wageningen, Netherlands), containing an EGNOS-augmented GPS receiver and a tri-axial accelerometer sensor. The cows were continuously kept on pasture at Carus animal facility in Wageningen, Netherlands (51°59’8’’ N, 5°39’11’’ E), and could move freely around as a single group during the experiment. Over the course of this period, we relocated the cows between three 0.32 ha pasture plots (sequentially five, six and six consecutive days in each plot). At every pasture switch the cows were housed inside the Carus facility for one night where they were offered fodder, so that they were not hungry at the start of a new pasture plot session. Furthermore, the cows were taken inside for milking and feeding every morning between 7:30 and 8:30 CEST and solely for milking every afternoon between 16:30 and 17:00 CEST. The time the cows spent on pasture was short enough to assume that the pasture did not increase in grass quality because of re-growth after grazing and only decreased in grass availability [27]. The short duration of the pasture sessions (approximately one day longer than when a commercial farmer would have moved the cows, as judged by the farm manager) ensured that the cows were not hungry, but only had to put more effort into foraging when time progressed. Furthermore, the collaring process did not put the cows under noticeable stress, more so because they were accustomed to continuously wearing a neck collar.
The sensors in the cows’ neck collars recorded GPS and accelerometer data during the experiment. The data were stored with a millisecond-accurate timestamp on a local SD memory card, which was replaced every one to five days together with the battery. GPS data were stored on the SD card with a one second interval. The accelerometer data were sampled with a variable frequency of 25-500 Hz, which were later down-sampled and linearly interpolated to a constant 32 Hz signal. Both the GPS and the accelerometer did not record data during some hardcoded multi-hour periods of inactivity, which were variable in duration and time of day, to save battery power. However, the time between GPS fixes was exactly 1 second in more than 99% of the cases. The precision of the GPS fixes was high, with 98% of the fixes having a Horizontal Dilution of Precision (HDOP) of less than two (a dimensionless unit; two is considered “excellent” precision). All GPS fixes with an HDOP of more than five, which were 0.5% of all fixes, were considered to be untrustworthy and filtered out of the final dataset. We also tested the accelerometer data for precision by placing the sensor on a stable, non-moving surface while it recorded for several minutes. The fluctuations in the recorded signal of all three accelerometer axes were small, 0.06 m s-2 between the lowest and highest value, and were considered negligible and thus ignored.
Activity (or behaviour) observations were conducted on work days from 25 April to 9 May 2017. A single person visually classified the activities using focal-animal sampling with a pre-defined ethogram (Table 1). All activity types in the ethogram (grazing, walking, standing, standing while ruminating, lying, lying while ruminating) were mutually exclusive. Each individual cow was observed continuously for ten minutes in the morning (10:00-13:00 CEST) and ten minutes in the afternoon (13:00-17:00 CEST), in random order, resulting in a total observation time of 1760 minutes. During the observations, the start and end times of each displayed activity type from the ethogram were recorded. We conducted these observations to acquire annotations for an activity classification model. Representative acceleration plots of the three axes for the different activity types are provided (see Additional file 2).
Table 1: Ethogram. Descriptions of the recorded, mutually exclusive activity types.
Activity
|
Description
|
Grazing
|
Foraging behaviour by chewing grass from the pasture whilst standing still or slowly moving with the head down
|
Walking
|
Taking at least two steps without grazing, either with the head up or down
|
Standing without ruminating
|
Standing on all four legs with head erect, without swinging its head from side to side and without ruminating
|
Lying down without ruminating
|
All four legs tucked underneath the torso or lying down on one side of its body without ruminating
|
Ruminating while standing
|
Masticating regurgitated feed, swallowing masticated feed or regurgitating feed while standing with head erect
|
Ruminating while lying down
|
Masticating regurgitated feed, swallowing masticated feed or regurgitating feed while lying down
|
We measured resource availability as dry matter grass biomass in kilograms per hectare, excluding stubble biomass. We determined time-varying biomass levels using a combination of field-measured biomass levels at specific time points, satellite-based biomass estimates derived from the Normalized Difference Vegetation Index (NDVI), and modelling of grass dynamics (see Additional file 3). Wind speed (m s-1, mean speed 10 m above ground) were recorded at 10 minute resolution during the experiments with a weather station on a grass pasture at the Veenkampen, Wageningen, Netherlands. This weather station is located one kilometre west of the pasture plots used for the experiments.
Data processing
We used the pre-processed 32 Hz, tri-axial accelerometer signal as input for the accelerometer feature extraction. First, we converted all the records in the three-dimensional accelerometer dataset to 21 dimensions using multiple geometric transformations, i.e., resultant vectors, angles, solid angles, volumes and areas (Table 2). These dimensions constitute all geometric transformations of angles and distances in one, two and three dimensions. Considering that tri-axial accelerometer readings describe the movement forces in three dimensions, geometric transformations make sense from a physics perspective. More transformations could be considered, but these may lack to provide additional information to the feature set. Second, we divided the resulting dataset into non-overlapping time windows. We tried all window sizes in the range of 1 until 30 seconds and optimized this window size as a hyperparameter regarding the activity classification performance, where 3 s turned out to be the optimal window size (see Additional file 5). For every time window we computed multiple statistics per accelerometer dimension per cow, e.g., mean, standard deviation, quantiles and Fast Discrete Fourier Transform (FFT) parameters (Fig. 2). These statistics were chosen to provide summary statistics about both the time-invariant and sequential aspects of the data, given that accelerometer data also includes patterns in the frequency domain regarding animal activity (e.g., head movement of cows during grazing has a strong cyclic behaviour). We computed the FFT with the base R stats package [28], of which we used the maximum FFT value as the dominant amplitude, the corresponding period of the dominant amplitude as the dominant period, and finally the sum of all squared FFT values as the spectral energy. Our list of computed statistics is not all-encompassing and more statistics can be thought of to describe patterns in the data. However, as these statistics were mainly used in the activity classification part of the analysis, we deemed the computed statistics sufficient when it resulted in a high performance during activity classification. Overall, computing all statistics for each dimension resulted in 210 accelerometer features per time window per cow.
[Please see the supplementary files section to view Table 2.]
We used the filtered 1 Hz GPS data as input for the GPS feature extraction. First, we transformed all the latitude, longitude coordinates to Cartesian coordinates by projecting them to zone 31N of the UTM system (EPSG 32631). Second, we extracted a number of individual GPS features from the projected GPS coordinates per time window per cow, related to speed, turning angle, tangential velocity, mean squared displacement, and first passage time (Table 3). The time windows were exactly the same as the time windows used in the extraction of the accelerometer features. Third, we extracted a number of group GPS features from the projected GPS data per time window per cow, related to group shape, group area, and distances and directions to other cows (Table 4). We determined which individual and group GPS features to compute by drawing fake GPS trajectories and animal clusters, after which we discussed which geometrical properties (e.g., tangential velocity: the linear speed of an animal moving along a circular path) could be extracted from these patterns. Furthermore, we computed ecological properties of animal trajectories that were known to us (e.g., Mean Squared Displacement: a measure of the deviation of the position of an animal with respect to a reference position over time) and searched the literature and animal movement related R packages for other ecological properties (e.g., First Passage Time: the time required for an animal to cross a circle with a given radius). We do not suggest that the provided list of computed features is all-encompassing, but we do suggest that spending time and effort in the engineering of features (or optimizing the architecture of a neural network in a deep learning approach) is an important part of our suggested framework. The more informative variation that is extracted from the raw data, the better the model could potentially perform and thus the better the quantified environmental contribution to animal movement matches reality. Overall, computing both the individual and group GPS features resulted in 38 GPS features per time window per cow.
Table 3: Individual GPS features extracted per time window and cow.
Dimension
|
Statistic
|
Description
|
Distance
|
Net gross ratio
|
Distance between first and last position divided by sum of distances of all segments
|
Speed
|
Mean
|
|
Standard deviation
|
|
Median
|
|
Minimum
|
|
Maximum
|
|
First quartile
|
|
Third quartile
|
|
Autocorrelation function index
|
Autocorrelation value at a lag of 1 second
|
Brownian motion scaling parameter
|
See Equation 1
|
Turning angle
|
ρ
|
Length of the mean resultant vector
|
Autocorrelation function index of the absolute turning angles
|
Autocorrelation value at a lag of 1 second
|
Absolute tangential velocity
|
Mean
|
|
Standard deviation
|
|
Median
|
|
Minimum
|
|
Maximum
|
|
First quartile
|
|
Third quartile
|
|
Autocorrelation function index
|
Autocorrelation value at a lag of 1 second
|
Mean Squared Displacement
|
Diffusion coefficient
|
The value of a in the fitted model on MSD values for τ from 1 to 6
|
Diffusion power coefficient
|
The value of b in the fitted model on MSD values for τ from 1 to 6
|
First Passage Time
|
Mean, 5m radius
|
|
Variance of log, 5m radius
|
|
Autocorrelation function index, 5m radius
|
Autocorrelation value at a lag of 1 second
|
Radius with maximum variance of log (integers from 1 to 10m)
|
|
Linear regression coefficient log radius vs. log mean FPT
|
|
Table 4: Group GPS features extracted per time window and cow.
Dimension
|
Statistic
|
Description
|
Net distances to other cows
|
Mean
|
|
Median
|
|
Minimum
|
|
# cows within 2m radius
|
|
# cows within 4m radius
|
|
# cows within 8m radius
|
|
# cows within 16m radius
|
|
All mean cow coordinates
|
Group elongation index, φ
|
Variance explained by the first principal component through all cow coordinates, minus 0.5 and times 2
|
Group area proxy
|
; where σ is the standard deviation of the first principal component values
|
Directions to other cows
|
ρ
|
Length of the mean resultant vector
|
Periphery index
|
Maximum difference between consecutive directions, minus and divided by
|
Data analysis
We used the accelerometer features and individual GPS features per time window per cow for which activity observations were undertaken as input data for the activity classification models (Fig. 2, Table 3), which we first converted to principal components. We linked the time-matched activity observations to these input data and used the activity type as output variable for the classification models. We trained a multi-class classification model for the activity types: grazing, walking, standing and lying down. As a second step after the main activity classification we also trained a binary classification model for ruminating, with an extra input variable that indicated standing versus lying down. We tried for both classification models a Support Vector Machine (SVM) with a Radial Basis Function (RBF) kernel and a one-against-one approach, implemented in the e1071 package for R 3.6.2 [28,29], and a Random Forest (RF) with 500 trees, implemented in the randomForest package [30]. To prevent overfitting due to autocorrelation in the data we randomly assigned each hour of the dataset into a train (80%) or test set (20%) and performed 5-fold cross-validation on the train set, which was also split per hour at each of the 5 cross-validation iterations [21]. To find the optimal hyperparameters for the models (number of principal components and time window size for both SVM and RF; cost, gamma and class weights for SVM; and mtry, sample size and node size for RF), we used an extensive grid search on a High Performance Cluster of Wageningen University, Netherlands (see Additional file 5). We started the grid search with a coarse resolution search that covered a large range of all hyperparameters, to make sure that the global optimum was covered and to get a feel for the performance landscape. We zoomed in with a finer resolution during a second grid search and finished with an even more zoomed in and finer resolution during a final grid search. We determined the optimal classification model and hyperparameters by selecting for the highest mean balanced accuracy during cross-validation (Equation 2). The classification models with the highest performance during cross-validation were then evaluated for performance on the test dataset. Finally, we used the models to predict the displayed activity type (grazing, walking, standing or lying down) and whether or not the cows were ruminating, for all the time windows and cows with available sensor data.
We computed the dataset for the environmental variable predictions per cow over one-hour time windows. The window size that is chosen has of course an influence on the results, as the effect of an environmental variable on animal movement data varies with temporal scales [31]. In short, the window size that is chosen represents the scale at which the animals’ behavioural decisions are made [31]. The choice of this temporal scale should therefore be chosen in line with the study’s aim and based upon ecological considerations, which are different for every study. We chose a window size of one hour for a semi-illustrative purpose and because it traded off the number of resulting data records (number of rows in the dataset after applying the one hour window) and the convergence of variables well; meaning that the resulting dataset consists of hundreds of records (thereby being enough for a data-driven machine learning approach) and each record was based on 1200 (one hour divided by 3 seconds) underlying records or more (thereby making sure that the inherent heterogeneity of animal movement is taken into account by averaging it out over a large enough period). The calculated variables consisted of multiple variable sets, based on the source of the data (GPS or accelerometer), organizational level (group or individual), transformation type, and variables conditional on foraging (Table 5). We did not consider variables conditional on other activity types than foraging, because the cows sometimes did not display one of the other activity types during a one hour time window. This resulted in a total of 548 variables per cow per one-hour time window. We standardized these variables (to zero mean and unit variance) per combination of day/night and cow ID to account for differences in nocturnal and diurnal activities of cows and individual differences in movement characteristics, group characteristics, and activities. These standardized variables were used as input for a principal component analysis, but were first one by one visually checked for symmetric unimodality by inspecting the histograms and normal Q-Q plots. Two of the 548 variables displayed signs of bimodality and eight variables appeared to be somewhat heavy-tailed. Due to the low number of variables that showed these deviations and due to the small severity of these deviations, we decided not to correct these ten variables and thus left all standardized variables untransformed. Moreover, symmetric unimodality is not an actual requirement of a principal component analysis, but it does result in a better centring and scaling of the principal components. After that we converted the standardized variables to principal components separately for the GPS and accelerometer variables and linked these principal components to the mean grass biomass, time since milking, and wind speed values per hour (see Additional file 4). To prevent overfitting of the model due to autocorrelation of the time series, we trained the model on the data of all cows from two of the three pasture plot sessions (n = 600, viz., number of rows in the train set) and tested the model on the data of all cows from the other pasture plot session (n = 259, viz., number of rows in the test set). We used the second pasture plot session as our test set, because its range of biomass values fell within the range of biomass values of the first and third pasture plot session.
Table 5: Calculated variable sets per cow over one-hour time windows.
Variable set
|
Statistic
|
Transformed data
|
Individual GPS
|
All statistics from Table 3
|
1 Hz GPS data
|
Proportion activity
|
Proportion
|
Predicted activity per three-seconds window (Table 1)
|
Individual GPS distribution parameters while grazing
|
Mean and standard deviation of log-transformed data
|
Median speed and median absolute tangential velocity per three-seconds window while grazing (Table 3)
|
Median group GPS
|
Median
|
Group GPS features per three-seconds window (Table 4)
|
SD group GPS
|
Standard deviation
|
Group GPS features per three-seconds window (Table 4)
|
Median individual GPS while grazing
|
Median
|
Individual GPS features per three-seconds window while grazing (Table 3)
|
SD individual GPS while grazing
|
Standard deviation
|
Individual GPS features per three-seconds window while grazing (Table 3)
|
Median group GPS while grazing
|
Median
|
Group GPS features per three-seconds window while grazing (Table 4)
|
SD group GPS while grazing
|
Standard deviation
|
Group GPS features per three-seconds window while grazing (Table 4)
|
Median accelerometer while grazing
|
Median
|
Accelerometer features per three-seconds window while grazing (Fig. 2)
|
SD accelerometer while grazing
|
Standard deviation
|
Accelerometer features per three-seconds window while grazing (Fig. 2)
|
To predict the environmental variables we built a Support Vector Regression (SVR) model with a RBF kernel and a Random Forest Regression (RFR) with 1000 trees on the train set with both GPS and accelerometer principal components, with only GPS components, and with only accelerometer components. These models are time-invariant, as they assume independence between the data records, and are particularly well-suited to model complex interactions between a large number of variables. To find the optimal hyperparameters for the models (number of principal components for both SVR and RFR; and cost, gamma and epsilon for SVR), we used a grid search (following the same procedure as during the grid search of the activity classification) on a High Performance Cluster of Wageningen University, Netherlands (see Additional file 4). We did not optimize any other RFR hyperparameter, because the performance improved barely compared to the default values during a trial analysis. We determined the optimal hyperparameters by selecting for the highest R2 on the test set (Equation 3). Ideally, (cross-)validation is performed before a test set evaluation to prevent overfitting in hyperparameter space, but the limited quantity of data records in our case study prevented us from setting aside more data from the train set. However, we prevented overfitting in hyperparameter space by not optimizing the hyperparameters of the RFR and by limiting the amount of hyperparameter values that were tested for the SVR.