Changing measurements or changing movements? Sampling scale and movement model identifiability across generations of biologging technology

TL;DR

This study evaluates how differences in sampling scale across generations of biologging technology affect the ability to accurately identify seabird movement patterns, emphasizing the need for synthetic data assessments.

Contribution

It demonstrates the impact of varying observational regimes on movement model identifiability and provides guidelines for valid comparisons across different tracking technologies.

Findings

Model selection and parameter estimation can be biased when comparing data at different temporal scales.

Degrading data to the lowest resolution ensures more valid comparisons.

Synthetic data assessments are recommended before analyzing disparate datasets.

Abstract

1. Animal movement patterns contribute to our understanding of variation in breeding success and survival of individuals, and the implications for population dynamics. 2. Over time, sensor technology for measuring movement patterns has improved. Although older technologies may be rendered obsolete, the existing data are still valuable, especially if new and old data can be compared to test whether a behaviour has changed over time. 3. We used simulated data to assess the ability to quantify and correctly identify patterns of seabird flight lengths under observational regimes used in successive generations of tracking technology. 4. Care must be taken when comparing data collected at differing time-scales, even when using inference procedures that incorporate the observational process, as model selection and parameter estimation may be biased. In practice, comparisons may only be valid when degrading all data to match the lowest resolution in a set. 5. Changes in tracking technology that lead to aggregation of measurements at different temporal scales make comparisons challenging. We therefore urge ecologists to use synthetic data to assess whether accurate parameter estimation is possible for models comparing disparate data sets before conducting analyses such as responses to environmental changes or the assessment of management actions.

Tables

Table 1. Table 1: A glossary of terms describing movement paths used in this study.

Term	Definition
Trip	A trip is assumed to be one foraging excursion, beginning when the animal leaves the nest site and ending when it returns. A trip is comprised of flights interspersed with (water) landings.
Flight	Flights are the subcomponents of a trip, the units of space or time between prey capture attempts, in which the bird is actively flying.
Step	In tracking studies of terrestrial animals, this is more commonly used to describe distance, rather than time, and again, represents the sub-unit of a trip. Here, we use interchangeably with flight.
Segment	A discrete time unit over which the wet/dry status of the bird is measured. These segments may be aggregated into longer intervals.
Interval	The period over which aggregation of one or more wet-dry segments occurs. In the interval, the number of wet and dry segments are recorded. Flights are comprised of integer numbers of consecutive completely dry intervals. For data at high (time) resolution data the segment and interval time scales may be the same.
Record	The counts of flight lengths (in intervals) within or across trips.

Table 2. Table 2: Overview of immersion logger datasets used in this study

Study	Species	Year	Agg. interval (s)	$N_{deployments}$	$N_{flights}$
BBA2002	Thalassarche melanophris	2002	600	1	1503
walb2004	Diomedea exulans	2004	10	39	3604
walb1998	Diomedea exulans	1998	15	17	878
walb1993	Diomedea exulans	1993	720	11	298
walb1992	Diomedea exulans	1992	3600	21	340

Table 3. Table 3: Model selection for observational data. Here we show the difference in BIC from the best performing model ( Δ Δ \Delta BIC), such that the best model has a value of 0. We also show the calculated model probabilities, based on Equation ( 2 ). Dataset identifiers correspond to Table  2 .

BBA2002	$\Delta$BIC	MP
exp	81.0	0
pareto	3626.9	0
gamma	0	1
qexp	108.0	0

Table 4. Table 4: Estimated parameters (with approximate 95% confidence intervals) for the gamma distribution for each of the 3 observational data sets. Additionally, we show the calculated theoretical mean and variance based on the fitted parameters. Dataset identifiers correspond to Table  2 .

Gamma	$\alpha$ (shape)	$\beta$ (rate)	mean	variance
BBA2002	1.38 (1.29 – 1.47 )	1.15 (1.06 – 1.24)	1.20	1.05
walb2004	0.314 (0.293 – 0.334)	0.392 (0.363 – 0.422)	0.799	2.04
walb1998	0.0730 (0.0385 – 0.107)	0.170 (0.133 – 0.206)	0.430	2.53

Equations

$$\ell (\mathbf{\theta }|\mathbf{r})=\sum _{j=1}^J{d}_j\log [\mathrm{P}(j|\mathbf{\theta })]$$

$$P(M_i)\approx$$

$$=\frac{e^{-\frac{1}{2}[BIC(M_i)-BI{C}_{min}]}}{{\sum }_{r=1}^R{e}^{-\frac{1}{2}[BIC(M_r)-BI{C}_{min}]}}$$

$$f(x)=\{\begin{array}{ll}0,& \text{for }x<x_0\\ \lambda e^{-\lambda (x-x_0)}& \text{for }x\ge x_0\end{array}$$

$$\hat{\lambda }=\frac{1}{\bar{x}-x_0}$$

$$f(x)=\{\begin{array}{ll}0,& \text{for }x<x_0\\ \frac{\alpha x_0^{\alpha }}{x^{\alpha +1}}& \text{for }x\ge x_0\end{array}$$

$$\hat{\alpha }=\frac{n}{{\sum }_i(\ln x_i-\ln x_0)}$$

$$f(x)=\{\begin{array}{ll}0,& \text{for }x<x_0\\ \frac{{\beta }^{\alpha }}{\Gamma (\alpha )}(x-x_0{)}^{\alpha -1}{e}^{-\beta (x-x_0)}& \text{for }(x)\ge x_0.\end{array}$$

$$f(x)=\{\begin{array}{ll}0& \text{for }x<x_0\\ (2-q)\lambda e_q(\lambda (x-x0))& \text{for }(x)\ge x_0.\end{array}$$

$$e_q(x)=[1+(1-q)x{]}^{\frac{1}{1-q}}.$$

$$\ell (\mathbf{\theta }|\mathbf{r})=\sum _{j=1}^J{d}_j\log [P(j|\mathbf{\theta })]$$

$$\mathrm{P}(R=j|\theta )=(1-j)$$

$$+\frac{1}{s}\left[{\int }_{s(j-1)}^{sj}xf(x;\theta )dx+{\int }_{sj}^{s(j+1)}xf(x;\theta )dx\right]$$

$$\mathrm{P}(R=j_{\min }|\theta )=(1+j_{\min }){\int }_{s{j}_{\min }}^{s(j_{\min }+1)}f(x;\theta )dx+\frac{1}{s}{\int }_{s{j}_{\min }}^{s(j_{\min }+1)}xf(x;\theta )dx$$

$$\ell (\mathbf{\theta }|\mathbf{r})=\sum _{j=1}^J{d}_j$$

$$+\frac{1}{s}\left[{\int }_{s(j-1)}^{sj}xf(x;\theta )dx+{\int }_{sj}^{s(j+1)}xf(x;\theta )dx\right])$$

$$-n\log \left(1-P(R=j_{\min }|\theta )\right).$$

$$\mathrm{P}(R=j|\theta )={\int }_{sj}^{s(j+1)}(x-j)f(x;\theta )dx+{\int }_{s(j+1)}^{s(j+2)}(2-x+j)f(x;\theta )dx.$$

$$\ell (\mathbf{\theta }|\mathbf{r})=\sum _{j=1}^J{d}_j$$

$$-n\log \left(1-\mathrm{P}(R=0|\theta )\right).$$

Full text

Changing measurements or changing movements? Sampling scale and movement model identifiability across generations of biologging technology

Leah R. Johnson

Department of Statistics, Virginia Tech, Blacksburg, VA 24061, USA

Philipp H. Boersch-Supan

Department of Geography, University of Florida, Gainesville, FL 32601, USA

Emerging Pathogens Institute, University of Florida, Gainesville, FL 32610, USA

Department of Integrative Biology, University of South Florida, Tampa, FL, USA

Richard A. Phillips

British Antarctic Survey, Natural Environment Research Council, Cambridge, CB3 0ET, UK

Sadie J. Ryan

Department of Geography, University of Florida, Gainesville, FL 32601, USA

Emerging Pathogens Institute, University of Florida, Gainesville, FL 32610, USA

Summary

Animal movement patterns contribute to our understanding of variation in breeding success and survival of individuals, and the implications for population dynamics. 2. 2.

Over time, sensor technology for measuring movement patterns has improved. Although older technologies may be rendered obsolete, the existing data are still valuable, especially if new and old data can be compared to test whether a behavior has changed over time. 3. 3.

We used simulated data to assess the ability to quantify and correctly identify patterns of seabird flight lengths under observational regimes used in successive generations of wet/dry logging technology. 4. 4.

Care must be taken when comparing data collected at differing time-scales, even when using inference procedures that incorporate the observational process, as model selection and parameter estimation may be biased. In practice, comparisons may only be valid when degrading all data to match the lowest resolution in a set. 5. 5.

Changes in tracking technology, such as the wet/dry loggers explored here, that lead to aggregation of measurements at different temporal scales make comparisons challenging. We therefore urge ecologists to use synthetic data to assess whether accurate parameter estimation is possible for models comparing disparate data sets before planning experiments and conducting analyses such as responses to environmental changes or the assessment of management actions.

1 Introduction

Movement is an integral part of the foraging behavior of many animals (Nathan et al., 2008), and can account for much of their daily energy expenditure (Chai et al., 1999). Observed patterns of movement are determined not only by evolved behaviors and intrinsic state (age, sex, body condition, etc.), but also by environmental conditions (climate, oceanography, etc), and by prey abundance and distribution (Hays et al., 2016). Understanding how different movement strategies affect foraging success may provide insight into the processes underlying survival and reproductive success and, ultimately, population dynamics. Robust quantification of these patterns and how they change through time is required to meet this goal (Crossin et al., 2014).

Over recent decades, advances in animal tracking and biologging technology have provided an enormous amount of increasingly precise measurements of animal movement paths (Phillips et al., 2007; Block et al., 2011; Hussey et al., 2015; Kays et al., 2015; Hays et al., 2016). Researchers use satellite transmitters or data loggers to collect location information; these provide data at intervals which are often constrained by battery power or memory capacity (Fedak et al., 2002; Edwards et al., 2007; Shillinger et al., 2012). In large birds and mammals, particularly in recent years, the deployment of GPS loggers, and, in marine animals, of salt-water immersion or temperature-depth loggers, have generated a wealth of tracking data at high temporal resolution and at relatively low cost (Mackley et al., 2010; Block et al., 2011; Scales et al., 2016). However, data collected previously using VHF and satellite transmitters (platform terminal transmitters or PTTs), older GPS and immersion loggers, or by human observers are also available, albeit at coarser spatial and/or temporal resolution (e.g. Edwards et al., 2007; Froy et al., 2015). The increasing use of tracking and biologging technology has also been accompanied by initiatives to archive, share, and exchange animal tracking data (Birdlife International, 2004; Kranstauber et al., 2011). This wealth of existing data creates opportunities for informative comparisons between archived and new behavioral data. However, the different recording resolutions add complications both in terms of methodology and interpretation.

With the burgeoning of biologging and other ecological research, detailed observations are now available that span a time period that is relevant to the temporal scales of demographic processes, even for long-lived animals, as well as changes in the Earth’s climate (Hazen et al., 2013; Crossin et al., 2014). New research avenues have therefore opened for using biologging data to study how movement patterns may be changing across time, including in response to environmental variation (Hays et al., 2016). Most studies that deploy tracking devices on animals, such as seabirds, are usually aimed at answering broad ecological questions about habitat use and foraging behavior in one or a few successive years, as opposed to describing patterns of movement across time frames longer than a decade (but see Bogdanova et al., 2014; Carneiro et al., 2016). Consequently, device sampling intervals may be sub-optimal for learning about movement over the longer term in post hoc studies. Assessing whether movement strategies have changed requires robust methods and movement models that allow the synthesis of datasets collected at different temporal scales, with differing accuracy, and often with different research aims at the outset.

There are many ways to describe and quantify movement patterns of animals that depend on the type and quality of available data. Many models of foraging assume that organisms move diffusively, i.e., that animals perform uncorrelated Brownian walks as they search for food (Johnson et al., 1992). However, for most animals the Brownian assumption is clearly inadequate (Turchin, 1998). Super-diffusive descriptions of movement, such as Lévy walks or flights (Shlesinger et al., 1995; Viswanathan, 2010; Watkins et al., 2005) or intermittent search strategies (Bénichou et al., 2006, 2007), which describe movement as small jumps interspersed with occasional longer jumps are popular alternatives to standard diffusion models as they allow for more complex patterns. Lévy walks, which model movements with step lengths determined by a power-law distribution, were first applied in ecology to describe the foraging strategies of wandering albatrosses, Diomedea exulans (Viswanathan et al., 1996), and have since been used to describe search or foraging strategies across many different biological systems (e.g., see references in Edwards et al., 2007). They were also shown theoretically to represent optimal search strategies for revisitable targets when the targets are fractally distributed (Viswanathan et al., 1999). However, the validity of Lévy flights as descriptions of animal movement foraging is hotly debated in the ecological literature (Buchanan, 2008; Humphries et al., 2010; Reynolds, 2012; Travis, 2007; Edwards et al., 2007; Auger-Methe et al., 2011; Viswanathan et al., 2011). For instance, although the initial study on albatrosses indicated a Lévy pattern of foraging (Viswanathan et al., 1996), after correcting and augmenting the original data, and utilizing improved statistical methods, a later study by Edwards et al. (2007) showed that the Lévy flight model was not supported; instead, flight times were more likely to be gamma distributed. Subsequent studies claim new evidence for Lévy like behavior in certain marine predators (Sims et al., 2008; Humphries et al., 2010; Sims et al., 2012; Hays et al., 2012; Focardi and Cecere, 2014; Reynolds et al., 2016), and that humans exhibit more complex behaviors (González et al., 2008). Further studies on albatrosses have concluded that foraging patterns of some (although far from all) individuals are well described by modified Lévy flights or Brownian movement in various contexts, and further concluded that birds utilizing this method are able to consume considerably more prey than they need to satisfy their own energy requirements (Humphries et al., 2012). Thus the evidence on Lévy flights in nature is decidedly mixed.

The controversy surrounding the Levy foraging hypothesis has focused both on the theoretical justification of this process model, as well as the statistical procedures used to distinguish Lévy walks from other random walks (Benhamou, 2007; Auger-Methe et al., 2011; Plank et al., 2013). A substantial body of literature has dealt with different fitting approaches and goodness-of-fit measures (Plank and James, 2008; White et al., 2008; Auger-Methe et al., 2011; Edwards et al., 2012), however the question of parameter identifiability or - for practical purposes - estimability has received considerably less attention (but see Auger-Methe et al., 2011; Auger-Méthé et al., 2016).

In this study, we use foraging data from albatross species collected a decade apart to explore how the changes in logger technology (and hence the scale and mode of sampling), modeled distributions (statistical fitting), and the treatment of both data and distributions, may influence the findings, and our ability to infer and compare behavior over time. Our analyses focus on a particular type of data from loggers which detect and record saltwater immersion, providing information on wet and dry periods (so called immersion loggers) (Edwards et al., 2007; Mackley et al., 2010). Although a geographic location was not available in some of the earlier deployments, PTT or GPS data have been collected concurrently with the immersion data in the past 20 years, providing improved insights into movements and habitat use (Mackley et al., 2010; Scales et al., 2016; Carneiro et al., 2016). Previously, a major consideration was memory capacity, which led to alternative ways of sampling and storing data, which were aggregated at different time-scales on the device during the deployment. The aims of our study were to evaluate model and parameter identifiability for different generations of immersion loggers, using synthetic datasets reflecting different sampling regimes. We further investigated parameter estimation for actual data collected in the wild.

2 Methods

In order to determine to what extent the inference of underlying foraging patterns is influenced by the data collection method, we combine a simulation study with an analysis of two suites of data on flights and water landings at sea collected a decade apart, 1992-1993 and 2002-2004, from wandering Diomedea exulans and black-browed albatrosses (Thalassarche melanophris). Since these data comprise segments of behaviors that have variously been referred to as steps, trips, tracks, flights, etc., a glossary is provided in Table 1.

2.1 Inference Procedure

All flights are assumed to come from one of four possible distributions: (shifted) exponential; (shifted) gamma; (shifted) q-exponential; pareto. Details of the distributions are given in Appendix B. These true flights are then resampled with (real or virtual) data loggers that discretize or aggregate the flights.

We use two approaches to infer the parameters of the underlying process from the data. One is to take a “naive” maximum likelihood approach, i.e., ignore the observational process, and instead assume that the observed data are drawn without noise from the underlying distribution. Another approach is to use a multinomial maximum likelihood approach that explicitly models the observational process (Edwards et al., 2007). In this case, the log-likelihood of the parameters $\boldsymbol{\theta}$, given a record $\mathbf{r}$ (a set of observations, see Table 1) takes the general form

[TABLE]

where $d_{j}$ is the number of recorded flights of length $j\in[1,J]$, and $\mathrm{P}(j|\boldsymbol{\theta})$ is the probability of observing a flight of length $j$ given the underlying flight time distribution and observation process.

One assumption of this multinomial model is that there is some biological lower limit to the possible flight in terms of the length of time spent dry. For instance, if a bird extended its foot out of the water to scratch its head, the logger would record that event as a dry interval; however, ideally, these events would be excluded from any analysis of flights. Following Edwards et al. (2007); Reynolds et al. (2016) and others, we use a lower limit of flight duration (part of an overall trip) of 30 seconds, on the biological assumption that this is not likely to be a flight to a different food patch. This lower limit to flight time is built into the exponential, gamma, and q-exponential distributions as a shift, and into the pareto as the lower set point (See Appendix B for details).

2.2 Simulation Studies

We used a suite of simulations to explore the effect of the different logger sampling schemes (specifically the time scales over which data are aggregated) on our ability to correctly infer parameters values of a known model and to choose the true model if we treat it as unknown.

First we generated a series of “true” flights drawn directly from the known distributions without the observation process. For each of the four distributions, we specified four parameter sets for a total of 16 underlying flight time distributions. When possible, we chose parameters so that the theoretical means between the 4 sets of parameters corresponded between distributions, to ensure that the scales of the processes were comparable. For each of the 16 flight distributions we created 10 simulated data sets of length 3000 (i.e., 10 sets of 3000 flights).

For each of the 160 simulated data sets, we then “observed” the data using our two most extreme sampling regimes, that is intervals of either 1 hour or 30 seconds corresponding to the sampling intervals used in field deployments in 1992 and 2004, respectively (Edwards et al., 2007). More details of the algorithm are in Appendix C.3. This resulted in 320 simulated data sets of length $\leq 3000$ (as the aggregation can result in a subset of flights being labeled as non-flights and thus discarded). These simulated data were used in the following two simulation studies.

Parameter Identifiability

Using the 320 simulated data sets, we attempted to infer the parameters from the underlying flight-time model corresponding to the one that generated the data. We used both a naive maximum likelihood estimate (MLE) (i.e., one that excluded the observational process) and the multinomial with the appropriate observational process. For instance, if the underlying model was an exponential, we fit the exponential model, only. The inferred parameters were then compared to the true parameters that generated the flights.

Model Identifiability

Using the 320 simulated data sets, we fit all four of the possible flight distribution models using of the exact (multi-nomial) likelihood. For each of the 320 data sets, we calculated the Akaike Information Criterion (AIC, Akaike, 1973; Bozdogan, 1987), the Bayesian Information Criterion (BIC, Raftery, 1986; Schwarz et al., 1978), and the approximate model probabilities based on BIC (Burnham and Anderson, 2004, and given by Eqn. (2), below), and used these to select the best model for each data set.

2.3 Immersion Data Analysis

Observational data on flights and water landings were obtained from immersion loggers deployed on the legs of individual wandering albatrosses and black-browed albatrosses from Bird Island, South Georgia (54∘ 00’S, 38∘ 03’W). Multiple types of loggers were deployed between 1992 and 2004. The data are summarized in Table 2.

Flight length calculation from immersion data

Wet/dry records were parsed at the highest temporal resolution for each type of logging device, before flight lengths were calculated by merging consecutive time periods recorded as dry (Supplementary Materials A).

Model fitting and model selection

All four flight distributions under the the multinomial likelihood with appropriate discretization and aggregation parameters, were fit to all sets of data noted in Table 2. Models for each data set were ranked via BIC. Further, approximate model probabilities $P(M_{i})$ (based on BIC, Burnham and Anderson, 2004) were calculated as:

[TABLE]

where $R$ is the total number of models being considered, and $BIC_{min}$ is the minimum BIC value across those models. The second expression is more numerically stable, and so is the one we use in our calculations

Comparing the flight length distributions between years and species

After selecting the best fitting model via BIC, we examine model fit by plotting the theoretical quantiles verses the observed data quantiles. We then compared the estimates of parameters for our three data sets. In the current likelihood framework we obtain point estimates for all parameters. We can then use these parameters to estimate the means/medians and variances between the fitted models.

3 Results

3.1 Simulation Studies

Parameter Identifiability

Using the 320 simulated data sets, we attempted to infer the parameters from the flight-time model corresponding to that which generated the data using both the naive likelihood (excluding observational process) and the exact multinomial likelihood with the appropriate observational process. Results for the exact likelihood are shown in Figure 1 and the naive likelihood in 2.

Overall, parameter estimates were much more precise for the data recorded at high frequency (i.e., 10 sec sampling), regardless of the underlying true distribution or the scale of the true process. This is because the true lengths of dry periods are recorded with high resolution when data are recorded at this high frequency. However, even for the higher resolution data, using the naive likelihood can bias parameter estimates for some cases of the Pareto and q-exponential distributions. This is probably because even the small amount of truncation of dry periods of 30 seconds or less changes the expected ratio of small flights to longer flights, biasing the estimates. Using the exact likelihood helps to account for these shifts.

In contrast, estimating parameters for a 1 hour integration step is difficult even when using the exact likelihood if the mean/median flight times are on the order (or less) of the integration period (i.e. if true flight times are less than $\sim 5$ hours in our simulations). Again, use of the exact likelihood can provide better results, although for the Pareto and q-exponential both approaches perform poorly if the true median flight duration is short and the integration interval is long.

Model Identifiability

Using the 320 simulated data sets, we fit all four of the possible flight distribution models. Focusing on the results from the multinomial approach applied to the appropriate observational process, In Figure 3 we show the proportion of times that the model is identified correctly under the two observation schemes across all 4 flight distribution models. The selected model is indicated when the true model is incorrectly identified.

For the short sampling interval (10 sec), the true model almost always has the highest probability. The exception is the q-exponential model, for which there is not a consistent best model if the true median flight length is one hour. Different patterns are apparent for the long (1 hr) integration interval. In this case, only for data generated from the gamma distribution does the true model have high model probability across all parameter settings. For the q-exponential, the true model is chosen if the true median flight time is ¿1hr (like in the short integration case). If the true underlying models are either Pareto or exponential the true model is never selected. Instead, data from the exponential distribution are always classified as gamma, and the data from the Pareto are either classified as q-exponential or as gamma. These results are congruent with the results from the parameter identifiability simulation experiment: when parameters can be well estimated, that model is likely to be correctly identified as the true model, whereas if the parameter estimates are poor, the patterns in the data that result from the observation process are no longer consistent with the true model, and are better described be one of the alternative models.

3.2 Immersion Data Analysis

Flight length calculation from immersion data

Across both species, and irrespective of the observation regime, the model that is most consistent with the observed data is the gamma distribution (Table 3). This is in line with previous results on a subset of the data (Edwards et al., 2007). Based on gamma Q-Q plots (Figure 4) for all three of the data sets, the fitted gamma distributions appear to be reasonable for data both from black-browed albatrosses in 2002 and wandering albatrosses in 2004, although both exhibit heavier tails than would be expected from the gamma distribution. The fit for the data from wandering albatrosses in 1998 is much poorer, and is under-estimating the number of short flights. Because the fit is relatively poor, directly assessing whether or not the patterns are consistent across time periods, or comparing between species, should be approached with care.

4 Discussion

Quantification of the movements of animals, such as seabirds, provides insights into foraging and migration behavior, the underlying drivers of movement, how movement and behavior may change over time in response to these drivers, and the consequences for individual performance and population dynamics (Crossin et al., 2014; Hays et al., 2016). The continuing development of new biologging technology for monitoring animal movement has greatly increased the resolution and quality of the data available, increased sample sizes, and reduced the effort required in the field, particularly for obtaining long time-series. These data represent invaluable archives for reconstructing historical movement patterns of animals for comparison with more recent observations. They provide a window into the past for understanding animal movements and the influence of changing environmental conditions, including the abundance and distribution of prey (Seco et al., 2016; Pereira et al., 2017). As animal movement databases grow (Birdlife International, 2004; Kranstauber et al., 2011), so do the opportunities for historical comparisons. However, as the resolution and accuracy of tracking devices has changed over time, these comparisons must be made with care to ensure robust interpretation.

In this paper, we used simulated data to assess the ability to quantify and correctly identify patterns of flight lengths under observational regimes that correspond to the range of older and more recent immersion logger technology for recording landings of foraging seabirds at sea. These simulation experiments are the optimal approach for evaluating new statistical methods – if it is impossible to reconstruct the true parameters from a known distribution, then the inference method will almost certainly be unreliable when applied to experimental or observational data. Further, this approach allows us to examine the impact of the observational method on the resulting conclusions, and to identify ways of comparing and combining disparate data sets to maximize their value.

Using simulated data from a set of four underlying flight length distributions that have been hypothesized to describe the flight distributions of seabirds (exponential, gamma, Pareto (corresponding to a Lévy flight), q-exponential) that are then “observed” using a sampling regime typical of immersion loggers deployed in the field, we were able to test the effectiveness of our statistical methods. In particular, we focused on the extent to which incorporating truncation and aggregation, which are part of the sampling procedure, into the likelihood estimation was necessary in order to determine accurate parameter values and correctly identify the underlying model. The results indicated that the inference procedure using the exact likelihood performs as well or better than the naive likelihood in all cases, i.e., accounting for the observation process improves our ability to both identify the model and estimate parameters. This improvement comes with a computational cost, as evaluating the exact likelihood is slower, and some tuning of the maximization procedure is required for individual data sets to achieve convergence. For low and medium resolution data (or if a Pareto distribution is considered, regardless of the resolution) the computational costs are worthwhile. For the very high resolution data (at least every 10 secs.) available from loggers in recent years, it may be sufficient to use the naive likelihood, with no need to incorporate the observational process.

Even when using the likelihood that incorporates the observational process, care must be taken when analyzing data that have been aggregated at time scales that are much longer than the events of interest. In these cases, the parameters can be significantly biased, and a model different from that used to generate the data may be chosen as the best model. In practice, it may therefore be impossible to accurately compare flight patterns from loggers that provide aggregated data at coarse scales from those that provide fine scale data. Instead, the latter may need to be degraded (reaggregated) to the coarser scale to determine if patterns from the two regimes are at least consistent, even if it is not possible to determine if the parameters are the same. This also puts a constraint on the biological questions that may be compared between data taken at different resolutions. For instance, short scale inferences about foraging intervals within a food patch may be unreliable whereas inferences about longer scale movement between patches may be accessable.

Based on the simulated data, we were not able to identify the process model underlying the “observations” with the coarse sampling regime consistently and accurately. Thus, we must be cautious when attempting to infer whether or not this particular aspect of the foraging strategy of the albatross has changed over the past two decades based the type of data at hand. Even for the higher temporal resolution data that we present here, the lack of model fit indicated by the QQ plots (Fig. 4) is concerning. In particular, there are more long flights than would be typical for the best fitting gamma model. The question is why would this be the case? In some cases, where concurrent location data are available, we may be able to determine that some longer flights may not represent foraging behaviors, and can be excluded. This was the case for a proportion of the data for which we had concurrent location data. The longer flights could also indicate individual birds that are not exhibiting foraging behavior (for instance attempting to fly out of a storm). In the current analysis we have treated the behavioral state as known, such that the wet status corresponds to feeding/handling attempts and dry to flying foraging. Thus we have not utilized a more complex statistical approach, such as state-space modeling (Patterson et al., 2008) that can allow concurrent estimation of behavioral state. If the observed longer flight patterns are a result of a separate non-foraging flying, these methods may be useful for identifying them.

Another possibility is that the mismatch between the data and the models is a symptom of inter-individual differences in behavior, or otherwise more complex behavior than the simple models here allow. If this is the case improving the models themselves, as well as the statistical techniques to analyze them, will be a more fruitful way forward. Even in the case of developing new models for the underlying behavior, it may be that direct parameterization of all model components is not possible. Instead, the quantified patterns explored here could be directly compared with model outputs, for instance emergent flight lengths from am optimal foraging or individual based model. Although model parameterization and validation of more complex models can be challenging, they can allow us to better understand why we see particular patterns and to better predict how behavior may change into the future.

Acknowledgments

This project was funded by an NSF grant (PLR-1341649) to L.R.J. and S.J.R. We thank Andrew Edwards for sharing his MLE code.

Author contributions statement

L.R.J., P.H.B, R.A.P., and S.J.R. conceived the ideas; P.H.B. and S.J.R. processed raw data; L.R.J. and P.H.B. performed data simulations and fitting; L.R.J., P.H.B, R.A.P, and S.J.R. analyzed the results. All authors wrote and reviewed the manuscript.

Data access

All observational datasets used in this study are available from the Polar Data Centre at the British Antarctic Survey, Cambridge, UK (polardatacentre@bas.ac.uk). Simulated datasets and associated simulation and analysis code will be deposited on http://zenodo.org and assigned a digital object identifier upon acceptance of the manuscript.

Appendix A Parsing and cleaning of data

Immersion data (wet/dry records) were downloaded in the field and stored in a variety of file formats. We developed a suite of parsing functions to import these records into R and format the data for further analysis. All code for file parsing and data analysis will be deposited in an online repository.

A.1 Flight length calculations

Wet/dry records were parsed at the highest temporal resolution for each type of logging device, before flight lengths were calculated by merging consecutive time periods recorded as dry. In line with previous studies (Edwards et al., 2007), only dry segments with a duration over 30 seconds were counted as flights, to avoid counting preening events and other behaviors that might involve the leg with the logger extended out of the water.

Appendix B Distributions for step lengths

We compare 4 distributional models for flight times: shifted exponential, pareto, shifted gamma (which is related to an exponentially truncated pareto), and the shifted q-exponential. Below we give brief details on each of these distributions. Implementations of these distributions in R are included in the supplementary materials.

Shifted Exponential Distribution

The shifted exponential is a generalization of the exponential where the support has been shifted some positive amount $x_{0}$. That is the pdf, $f(x)$ is given by

[TABLE]

where $\lambda>0$ is the rate parameter. Similarly was can say that $X-x_{0}\sim\mathrm{Exp}(\lambda)$. The theoretical mean of the exponential is $1/\lambda$. If the data were observed directly, the maximum likelihood estimator (MLE) is given by

[TABLE]

where $\bar{x}$ is the sample mean.

Pareto Distribution

The Pareto (Type 1) distribution is is power law probability distribution defined above a lower limit, $x_{0}$. The pdf, $f(x)$ is given by:

[TABLE]

where $x_{0}>0$ is the minimum possible value (that also determines the scale of the process), and $\alpha>0$ is a unitless shape parameter. This formulation is the same as the typical power law used in Levy studies, but with the shape parameter redefined so that $\alpha=\mu-1$ so that the parameter is constrained to be positive instead of $>1$. The mean and variance of the Pareto are only defined for a subset of parameter values: for $\alpha<1$ the mean approaches infinity, and the variance is also infinite for $\alpha\leq 2$.

In this paper, we assume that we have a biologically defined lower limit, $x_{0}$, so we are only concerned with estimating the shape parameter, $\alpha$. If the data were observed without error, the MLE for $\alpha$ is given by

[TABLE]

where $n$ is the number of data points. If we were also estimating $x_{0}$, the MLE is simply $\hat{x_{0}}=\min_{i}x_{i}$, and this estimator would be plugged into Eqn. 6.

Shifted Gamma Distribution

The shifted gamma distribution is a generalization of the gamma where the support has been shifted some positive amount $x_{0}$. That is the pdf, $f(x)$ is given by

[TABLE]

where $\alpha>0$ is the shape parameter and $\beta>0$ is the rate parameter. Similarly was can say that $X-x_{0}\sim\mathrm{Gamma}(\alpha,\beta)$. The theoretical mean of the standard Gamma distribution in this case is given by $\alpha/\beta$ and the variance by $\alpha/\beta^{2}$. There are no closed form MLEs for both parameters of the Gamma distribution. Instead we find these numerically.

Shifted Q-Exponential Distribution

The q-exponential distribution is a generalization of the exponential distribution with heavy (possibly power-law) tails. The pdf, $f(x)$ of the shifted distribution is given by

[TABLE]

where $1\leq q<2$ is the shape parameter, $\lambda>0$ is the rate parameter and

[TABLE]

When $q=1$ we regain the exponential distribution. Similarly was can say that $X-x_{0}\sim\mathrm{Pareto}(q,\lambda)$. Here we restrict ourselves to the case where $1\leq q<2$ as this is the range of $q$ for which the support of the distribution is on $[0;\infty)$. As with the gamma, we determine the MLE for the parameters analytically.

Appendix C Multinomial likelihoods for the the biologger data

As in Edwards et al. (2007), we use a multinomial maximum likelihood approach to estimate the parameters of the underlying flight process while taking into account the observational process of the biologgers. Here we briefly present the likelihoods described in detail in the Supplementary Materials 1 and 2 from Edwards et al. (2007), with the equations generalized slightly to account for differences in sampling protocols.

Most generally, the log-likelihood of the PDF parameters $\boldsymbol{\theta}$, given a set of observations $\mathbf{r}$ take the general form

[TABLE]

where $J$ is the number of recorded flights of length. The form of the multinomial probability $P(j|\theta)$ depends on both the underlying continuous distribution as well as the observation protocol. We recognize and implement likelihood functions for two classes of observation protocols: 1) discretized only (2004-type likelihood from Edwards et al. (2007)) and 2) discretized and aggregated (1992-type likelihood from Edwards et al. (2007)). Both forms of the likelihood, as well as functions for creating data, are implemented in R and included as part of the supplemental materials.

C.1 Discretized only data (2004-type)

Data of these types consist of sequences of wet/dry indicators within short intervals of length $s$ seconds (10 sec in the wandering albatross data from 2004). The record, $R$, is defined to be the number, $j$, of consecutive dry readings in between two wet readings (see Section A for further details). We define $m$ as the minimum interval, in seconds, that constitutes a flight (e.g., 30 sec for the data from 2004). The minimum record length that we include as a flight in our data set is $m/s+1$, as records shorter than this can include flights that are shorter than $m$. Following Edwards et al. (2007)), we can write the probability of observing flights of length $j$ as:

[TABLE]

where $f(x;\theta)$ is the pdf of the underlying flight time distribution.

Following Edwards et al. (2007)), we exclude records of length $j_{\mathrm{min}}=m/s$ as these can correspond to flights of less than $m$. However, we must account for the dry intervals of $>m$ that we miss by excluding the records of length $j_{\mathrm{min}}$. The probability of obtaining a record of length $j_{\mathrm{min}}$ is given by

[TABLE]

We can then obtain the full likelihood for our data by inserting Eqn. 10 into Eqn 9 and subtracting $n$ (the number of records) times the log of 1 minus Eqn. 11:

[TABLE]

The pdf, $f(x;\theta)$, from any of the probability distributions described in Appendix B can be used in this log likelihood. This function is then minimized numerically to obtain estimates of the parameters of the distribution.

C.2 Discretized and Aggregated data (1992-type)

Given memory limitations in the older types of immersion loggers, data were aggregated. The loggers deployed on wandering albatrosses in 1992 tested for saltwater immersion every three seconds, and if at least half of the tests were positive in segments of length $\epsilon$ (15 sec for the wandering albatross data from 1992), the segment was counted as wet. The segments of length $\epsilon$ were then aggregated over a larger interval $s$ (1 hour), where only the number of $\epsilon$-length segments that were wet in $s$ are recorded. Thus for each interval of length $s$ a number of segments that are wet will be an integer number on $[0;s/\epsilon]$ ($s/\epsilon=240$ for the 1992 WALB data). Since it is impossible to discern the exact number or pattern of immersion events within the interval, the exact flight times cannot be distinguished. Thus, the record $R$ for this case is defined to be the number, $j$, of consecutive completely dry intervals in between two intervals with at least one immersion. As before we define $m$ as the minimum interval, in seconds, that constitutes a flight (e.g., 30 sec for the data from 1992). Thus the minimum flight length is shorter than the minimum record length, which is by definition 1.

Analogous to the previous section, and following Edwards et al. (2007)) we can write the probability of obtaining a record of length $j$:

[TABLE]

Similarly to the in the previous section, this equation is included in the expression for the multinomial likelihood probability and subtracted from the portion relating to $\mathrm{P}(R=0|\theta)$ to get the full log-likelihood:

[TABLE]

As for the previous section, and distribution of interest can be inserted, and the log likelihood is then maximized to obtain estimates of the parameters for the underlying flight time distributions.

C.3 Sampling from the multinomial data distribution

Obtaining samples from the multinomial data distribution is straightforward, as for both observation methods (discretized or discretized and aggregated) a flight is the number of intervals that are fully dry.

Take a “true” flight time draw from a known distribution, $\tau$. 2. 2.

Randomly choose the start time of the flight uniformly within the observation interval. That is, the start time is $t_{s}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mathrm{U}(0,s)$ where $s$ is the length of the interval, for instance 10 sec for the 2004 WALB data or 1 hour for the 1992 WALB data. 3. 3.

The length of the flight, in segments, is then calculated as $\tau_{\mathrm{obs}}=\mathrm{floor}((\tau-t_{s})/s)$. 4. 4.

If $\tau_{\mathrm{obs}}>m$ (i.e., of a minimum length to be counted as a flight), that flight is added to the record, $R$.