A Bayesian Markov model with P\'olya-Gamma sampling for estimating individual behavior transition probabilities from accelerometer classifications

TL;DR

This paper introduces a Bayesian Markov model with Pólya-Gamma sampling to estimate animal behavior transition probabilities from accelerometer data, accounting for environmental factors and providing detailed behavioral insights.

Contribution

The study presents a novel Bayesian Markov model incorporating Pólya-Gamma sampling for efficient inference of behavior transitions from accelerometer data, including covariates like habitat and weather.

Findings

Significant differences in transition rates among habitats.

Diurnal behavior patterns identified.

Weather influences on behavior transitions.

Abstract

The use of accelerometers in wildlife tracking provides a fine-scale data source for understanding animal behavior and decision-making. Current methods in movement ecology focus on behavior as a driver of movement mechanisms. Our Markov model is a flexible and efficient method for inference related to effects on behavior that considers dependence between current and past behaviors. We applied this model to behavior data from six greater white-fronted geese (Anser albifrons frontalis) during spring migration in mid-continent North America and considered likely drivers of behavior, including habitat, weather and time of day effects. We modeled the transitions between flying, feeding, stationary and walking behavior states using a first-order Bayesian Markov model. We introduced P\'olya-Gamma latent variables for automatic sampling of the covariate coefficients from the posterior distribution and we calculated the odds ratios from the posterior samples. Our model provides a unifying framework for including both acceleration and Global Positioning System data. We found significant differences in behavioral transition rates among habitat types, diurnal behavior and behavioral changes due to weather. Our model provides straightforward inference of behavioral time allocation across used habitats, which is not amenable in activity budget or resource selection frameworks.

Tables

Table 1. Table 1 : Means and 95% credible intervals for odds ratios for habitat coefficients for the feeding to feeding transition arranged by increasing mean estimated in a Bayesian framework for six greater white-fronted geese in March 2018. A significant difference between habitats is indicated by non-overlapping credible intervals which can be visualized in Fig. 4 . The values correspond to exponentiating the estimates of the coefficients depicted in the from feeding to feeding panel of Fig. 2

Covariate	Mean	95% CI
Woody Wetlands	0.80	0.38	1.48
Open Water	0.91	0.78	1.07
Herbaceous Wetlands	1.13	0.81	1.53
Grasses	1.40	1.20	1.62
Soybeans	1.42	1.19	1.69
Fallow/Idle Cropland	1.49	0.84	2.46
Corn	1.52	1.37	1.70
Other Crops	1.72	1.06	2.64
Wheat	1.75	1.01	2.79
Developed	1.82	0.83	3.59

Table 2. Table 2 : A selection of means, 95% credible intervals and proportion of samples with an estimate greater than 1 for the odds ratios for quantitative covariate coefficients estimated in a Bayesian framework for six greater white-fronted geese in March 2018. The quantitative variable has a significant effect on the transition probability if the credible interval does not overlap with 1 which corresponds to a credible interval in Fig. 3 not overlapping 0.

Transition	Covariate	mean	95% CI		Proportion of Samples > 1
	Maximum Temperature	1.26	1.14	1.39	1.00
	Minimum Temperature	0.83	0.75	0.91	0.00
Flight to Flight	Wind Speed	0.84	0.76	0.92	0.00
Flight to Stationary		1.23	1.04	1.45	0.99
Feeding to Stationary		1.18	1.06	1.31	1.00
Stationary to Stationary		1.20	1.13	1.27	1.00
Walking to Stationary	cos(time)	1.22	1.13	1.32	1.00
Stationary to Stationary	sin(time)	1.17	1.11	1.24	1.00

Equations

$$\begin{array}{c}{\mathbf{P}}_{nt}\equiv \left[\begin{array}{cccc}{p}_{n11t}& p_{n12t}& \cdots & p_{n1Jt}\\ p_{n21t}& p_{n22t}& \cdots & p_{n2Jt}\\ ⋮& ⋮& \ddots & ⋮\\ p_{nJ1t}& p_{nJ2t}& \cdots & p_{nJJt}\end{array}\right],\end{array}$$

$$\begin{array}{c}{\psi }_{nijt}\equiv log\left((p_{nijt})/(p_{niJt})\right)={\text{x}}_{nt}^{\prime }{\mathbf{\beta }}_{ij},\text{and}\\ p_{nijt}=(exp({\psi }_{nijt}))/(\sum _{k=1}^{J}exp({\psi }_{nikt})),\end{array}$$

$${\mathbf{\beta }}_{ij}={\left[\begin{array}{ccc}({\mathbf{\alpha }}_{ij}{)}^{\prime }& ({\mathbf{\zeta }}_{ij}{)}^{\prime }& ({\mathbf{\theta }}_{ij}{)}^{\prime }\end{array}\right]}^{\prime }.$$

$$\begin{array}{c}{\mathbf{\alpha }}_{ij}={\left[\begin{array}{ccc}{\alpha }_{1ij}& \dots & {\alpha }_{(N-1)ij}\end{array}\right]}^{\prime }\sim N(0,{\sigma }_{{\alpha }_{ij}}^2\mathbf{I}),\text{such that}\sum _{n=1}^N{\alpha }_{nij}=0,\\ {\mathbf{\zeta }}_{ij}={\left[\begin{array}{ccc}{\zeta }_{1ij}& \dots & {\zeta }_{Hij}\end{array}\right]}^{\prime }\sim N({\mu }_{ij}1,{\sigma }_{{\zeta }_{ij}}^2\mathbf{I}),p({\mu }_{ij})\propto 1,\\ {\mathbf{\theta }}_{ij}={\left[\begin{array}{ccc}{\theta }_{1ij}& \dots & {\theta }_{B-(N-1+H)ij}\end{array}\right]}^{\prime }\sim N(0,{\sigma }_{{\theta }_{ij}}^2\mathbf{I}).\end{array}$$

$$\begin{array}{c}L({\mathbf{\beta }}_{ij}|{\mathbf{\beta }}_{i,-j},{\mathcal{S}}_1,...,{\mathcal{S}}_N)\propto \prod _{n=1}^N\prod _{t=1}^{T_n}\prod _{j=1}^{J-1}{\left[e^{{\psi }_{nijt}}/(\sum _{k=1}^J{e}^{{\psi }_{nikt}})\right]}^{y_{nijt}},\\ \propto \prod _{n=1}^N\prod _{t=1}^{T_n}{\left[e^{{\eta }_{nijt}}/(1+e^{{\eta }_{nijt}})\right]}^{y_{nijt}}{\left[1/(1+e^{{\eta }_{nijt}})\right]}^{1-y_{nijt}},\end{array}$$

$$\begin{array}{c}\left[\mathbf{\beta }|\mathcal{A}\right]=\int \left[\mathbf{\beta }|\mathcal{S},\mathcal{A}\right]\left[\mathcal{S}|\mathcal{A}\right]d\mathcal{S},\\ =\int \left[\mathbf{\beta }|\mathcal{S}\right]\left[\mathcal{S}|\mathcal{A}\right]d\mathcal{S},\end{array}$$

Full text

Running Title: Markov model for ACC behavior transitions

1 Key-words

animal behavior, auxiliary variables, hierarchical models, multinomial logistic, multiple imputation

2 Introduction

Animals make decisions daily that can result in differential fitness (i.e., survival or reproductive success; Brown et al., 2004; Breed and Moore, 2015). The knowledge of an animal’s behavior provides insight into its decision-making process. Historically, behavioral studies of animal populations were conducted using methods such as direct observation of focal individuals or instantaneous scan sampling of the group (Altmann, 1974). However, data collection by direct observation prevents a comprehensive understanding of the decision-making process because it is limited to specific times and places when the animal is observable. The introduction of wildlife tracking devices has largely mitigated these constraints and provided observations over greater time and space. Recent improvements in tracking technology have allowed researchers to gather high frequency data over extended periods of time (e.g., >1 year), which has led to unprecedented insights into animal decision-making (Nathan et al., 2012; Leos-Barajas et al., 2017).

Behavioral inference from animal movement models has long been a research goal. Early work by Morales et al. (2004) showed the accuracy of movement state-space models based on location data only could be improved by including latent behavioral states allowing the model to switch between movement modes such as exploratory and encampment. Specifically, an exploratory state was assumed to be characterized by longer step lengths and smaller turning angles, whereas an encamped state was characterized by shorter step lengths and larger turning angles. State-space models considering more than two states typically require more information. Michelot et al. (2017) included constraints in the temporal sequence of latent behavioral states based on expert knowledge of the annual cycle of the study species in order to distinguish among four states. McClintock et al. (2017) used a three state model, and states were distinguished by the addition of an auxiliary data stream. State estimation in this context typically is accomplished via hidden Markov models (McClintock, 2017), which have proven to be an efficient method for estimating an unobserved sequence of a categorical variable, such as behavior, that is associated with the values of observed quantities (e.g., step length and turning angle; Zucchini et al., 2016). In general, the estimated states are not interpretable as true behaviors, but rather proxies as they arise from a mixture model clustering procedure (Leos-Barajas et al., 2017; Patterson et al., 2017; Michelot and Blackwell, 2019). Therefore, when “ground truth” data are available and inference on specific behavior states is the primary goal, behavioral inference is better achieved by building a classifier and then applying a behavior model to the classified states as we propose in this paper.

The inclusion of accelerometers in wildlife tracking devices has become increasingly common. An accelerometer is a tool for measuring an object’s acceleration (ACC), and when placed on animals, can be used to derive energy expenditure and behavior of the tagged individual independent of location information (Nathan et al., 2012). Data collected from accelerometers are substantially different from Global Positioning System (GPS) observations in quantity, resolution and quality. That is, ACC data can be collected at a high frequency throughout the life of the tracking device, which results in richer and relatively larger data sets. The frequency of ACC collection can range from nearly continuous to more widely spaced intervals and is typically set to the highest frequency possible before affecting battery performance. Importantly, unlike location tracking devices, accelerometers do not typically miss “fixes” because the accelerometer instrument does not require linkages with external equipment such as satellites or radio trackers. In addition, ACC data are collected in two or three axes of movement relative to the device position, which provides the ability to discriminate between behavioral states with similar trajectory profiles.

The high frequency of data collected by accelerometers often allows researchers to “ground-truth" acceleration profiles with visual observations of the animals. When “ground-truth” observations of behavior are available, behavioral classification of ACC data can be conducted using machine learning algorithms (Resheff et al., 2014; Chakravarty et al., 2019). These approaches have greater reliability at identifying more than two states when compared to methods that consider only location data (Resheff et al., 2014). The richness of the ACC data allows for the derivation of many features to build the classification model. The fitted classification model can be used to assign behavior labels to ACC data. Then the labeled data can be used to make inference on animal behavior in a second stage. Currently, researchers summarize the labeled data within a time scale of interest (e.g., days or hours) and use a traditional activity budget framework, such as linking proportions of observed behaviors to covariates through generalized linear mixed models (Broekhuis et al., 2014; Heurich et al., 2014). However, by analyzing the aggregate summarized data, inherent temporal structure and small scale processes in the data may not be fully utilized. The analysis of the activity budgets generally ignores temporal dependence in the data by modeling the proportions of behaviors separately. Rugg and Buech (1990) showed that a Markov model for behavior resulted in improved estimates of time allocation compared to traditional activity budget analyses. Additionally, the use of the labeled data ignores uncertainty associated with the chosen classification method.

We propose a two-stage framework in which we first build a behavioral state classifier and then explicitly model temporal dependence in the high frequency behavioral observations via a Markov model and multiple imputation. Multiple imputation is often used in missing data scenarios and averages inference across a suite of potential true data sets (Scharf et al., 2017; McClintock, 2017). We build a classifier using features of the accelerometer data then build potential data sets based on the classification probabilities or proportions. For a given behavior data set, our Markov model works directly on the behavioral state transitions and inherently assumes that current behavior depends on recent behavior (i.e., the Markov assumption states that given the most recent past, the current behavior is independent of the long-term past). The transition probabilities are modeled with a logistic function to link covariates to the probability of transition. Therefore, covariate effects retain the useful odds ratio interpretation as in multinomial logistic regression, and provide insight into how covariates affect the tradeoff between time spent in each behavior. By using additional habitat information from location data, we can infer habitat use differences among behaviors. Resource selection models determine selection of habitats by animals based on differences in frequency of use and availability. As compared to traditional resource selection frameworks, we can make these inferences without having to define an "availability distribution" (Hooten et al., 2017). Although these approaches are important for identifying frequently used habitat, they do not define use and availability with respect to behavior. Our Markov model framework allows one to consider behavior profiles for different used habitats. Our approach provides the ability to answer questions about differential use of selected habitats that may be apparent from the focal animal’s behavior. For instance, a resource selection study may identify two different agricultural crops as preferentially used, but these results do not identify the behavior associated with the crops or the differences in behavioral rates between the crops.

The relationship between energy expenditure and fitness is especially critical for migratory animals because the time allocation of behaviors during migration, such as feeding, can impact survival during migration and subsequent reproductive success (Harrison et al., 2011). We consider a long-distance migrant bird, the greater white-fronted goose (Anser albifrons frontalis), to demonstrate the inference capabilities of our Markov model. The migration route of white-fronted geese spans a wide range of habitats in mid-continent North America. It is known that geese use different habitat types for different behaviors (e.g., roosting in water and feeding in crops) and thus, we expect differences in behavior transition probabilities or rates depending on habitat (Krapu et al., 1995). Our Markov model has the ability to provide inference on potential differences in behavior transition rates among habitats that serve a similar purpose. For example, feeding may primarily occur in agricultural fields, but rates of feeding may differ by crop type. Variability in rates of behavior transitions may also be related to environmental factors such as weather. Although activity budget analyses during winter have not found strong effects of weather (Ely, 1992), studies during spring have found important relationships between weather and the timing of migratory movements (Fox et al., 2003), and previous work has demonstrated heterogeneity in the movement of the geese (Hooten et al., 2018). We include weather variables comprising temperature and wind because we anticipate that these features can explain variation in behavior transitions during spring migration (e.g., favorable weather conditions increase rates of flight).

We implemented our Markov model in a Bayesian framework. At the core of our Markov model is multinomial logistic regression. There is a rich literature on Bayesian estimation of logistic models because sampling from the posterior distribution requires Metropolis-Hastings Markov Chain Monte Carlo algorithms that can be difficult to implement (tune) for moderate numbers of covariates or low observed frequencies of categories (Albert and Chib, 1993; Holmes and Held, 2006; Frühwirth-Schnatter and Frühwirth, 2010; Polson et al., 2013). To mitigate this challenge, automatic sampling strategies have been developed that rely on augmentation of data with latent variables. For example, Albert and Chib (1993) introduced the use of truncated normal random variables for category membership using probit link regression. Hooten et al. (2010) used the Albert and Chib sampling for a multinomial movement model. Holmes and Held (2006) expressed a multinomial logistic model as a product of binary models and used truncated scale mixture normals to define category membership. Alternatively, Polson et al. (2013) developed a class of Pólya-Gamma distributions for automatic sampling of Bayesian logistic regression models that is faster and provides exact inference when compared to the scale mixture latent variables. The introduction of Pólya-Gamma latent variables induces simple Gibbs updating steps for the logistic regression coefficients when using a normal (Gaussian) prior. Additionally, Polson et al. (2013) showed the Pólya-Gamma scheme is significantly faster than other data augmentation methods and more efficient than Metropolis Hastings for mixed effects logistic models. Note that Pólya-Gamma latent variables have been used to estimate coefficients of a transition probability matrix of a hidden Markov model for rainfall data (Holsclaw et al., 2017) and to learn dependence structure in topic models for text mining applications (Chen et al., 2013; Linderman et al., 2015; Glynn et al., 2019). Pólya-Gamma latent variables have yet to be used for analysis of animal behavior data. Moreover, to the best of our knowledge, this is the first comprehensive example of Pólya-Gamma latent variables and multiple imputation within a Markov model.

The presented framework provides rich inference about covariate effects including information related to the location from GPS data on the sequence of observed behaviors derived from ACC data. We expected to estimate behavior transition probabilities for greater white-fronted geese during spring migration that confirm prior knowledge of the diurnal pattern of behavior and habitat use. The coefficients should reflect an increase in stationary behavior overnight in open water and wetland habitats while also suggesting more movement during the day. Additionally, we hypothesized that the estimates would exhibit variability among habitat types, suggesting variability in time allocation. Furthermore, we expected behavioral transitions associated with flight and feeding to be influenced by inclement weather more than walking and stationary behaviors. Our method provides a fine-scale picture of the behavioral decision-making process that is driven by both ACC and GPS data rather than mechanistic drivers of movement.

3 Materials and Methods

3.1 Dataset

Greater white-fronted geese, hereafter white-fronted geese, migrate from wintering areas in the southern US (i.e., Arkansas, Louisiana, Mississippi and Texas) to breeding areas in Alaska and northern Canada (Fig.1; Baldassarre, 2014). Individual female white-fronted geese were captured using rocket nets between December 2017 and January 2018 in Texas. Captured individuals were each fitted with a U.S. Geological Survey metal leg band and solar-powered Ornitela (http://www.ornitela.com) neck collar, comprising GPS, ACC and Global System for Mobile communications (i.e., for daily data upload) technology. Age, sex and morphometric measurements (normal wing cord [cm], head, culmen, tarsus, middle toe lengths [mm] and mass [nearest 0.1 kg]) were also collected. We set tracking devices to obtain ACC values at 10 Hz for 3 seconds at 6 minute intervals and GPS locations at 30 minute intervals with data upload every 24 hours. To demonstrate our modeling approach, we used GPS and ACC data from six white-fronted geese subset to 1-31 March 2018, which comprises a portion of the spring migration period (Fig.1). For simplicity, we chose a time period when all of the geese were migrating within the U.S. rather than defining the specific dates associated with the geese leaving wintering habitat and arriving at breeding areas. The average number of ACC fixes per individual was 7,257 with a range of 7,054 to 7,319.

The location data were used to determine habitat and weather factors experienced by geese during spring migration. The U.S. Department of Agriculture maintains a raster of habitat and crop data throughout the contiguous lower 48 states (CropScape; https://nassgeodata.gmu.edu/CropScape/), which is updated annually. We assumed the habitat at the GPS fix was best represented by the CropScape grid cell containing the observed point. The original CropScape categories were combined into fewer groups (supplementary materials S.1). The National Centers for Environmental Prediction’s North American Regional Reanalysis (NARR; https://rda.ucar.edu/datasets/ds608.0/) data provide high resolution historical weather data in space and time. NARR is available eight times daily (i.e., every three hours) on an approximately 32 km grid. We assumed that the weather at the time of GPS fix was best represented by the NARR value corresponding to the nearest time and grid cell containing the observed point. We further assumed the weather and habitat at the time of an ACC fix were best represented by the values assigned to the most recent GPS location. The weather variables obtained were temperature [K], wind direction [ °] and wind speed [ $m/s^{2}$ ]. Temperature was further summarized to daily minima and maxima. Weather variables were standardized to have mean 0 and standard deviation 1. We assumed diurnal variability in behavior patterns because geese move most predictably near dawn and dusk between roosting and feeding areas; thus, we included such variability by using the local solar time of GPS fixes. We captured the diurnal behavior of the geese through two continuous covariates with a 24 hour period calculated from the local solar time of day (in seconds) by $cos(2\pi(seconds)/86400)$ and $sin(2\pi(seconds)/86400)$, referred to as $\cos(\mbox{time})$ and $\sin(\mbox{time})$. That is, $\cos(\mbox{time})$ is a representation of night (high values) and day (low values), while $\sin(\mbox{time})$ represents the first half of the day (high values) and the second half of the day (low values).

3.1.1 Behavior Classification

We summarized the ACC fixes into the 52 features described in the appendix of Resheff et al. (2014). The features were used to classify the ACC observations into four behavioral categories (i.e., flying, feeding, stationary and walking) using a random forest model. We chose a random forest model for its efficient handling of the large feature space. By using the random forest model, we assume the behavior labels are conditionally independent in time given the ACC data. The number of variables used to build the classification trees was varied from 1 to 15 and chosen by repeated 10-fold cross validation. The final model used 4 variables and had 96.5% accuracy on the training data set acquired by video recording. The random forest model was used to predict classification probabilities for the ACC fixes from the time period of interest. From the classification probabilities, we constructed $M=200$ possible data sets for multiple imputation. We discuss the details of imputation in the Model Fitting section. For more details on fitting the classification model, see supplementary materials (S.2).

3.2 Model

3.2.1 Behavior Transition Model

Each imputation data set is a regularly spaced time series of categorical behaviors, $\mathcal{S}_{nT_{n}}=\{s_{n0},...,s_{nT_{n}}:s_{nt}\in\{1,...,J\}\}$, where $s_{nt}$ is the observed behavior category for individual $n$, $n=1,...,N$, and time point $t$, $t=0,1,...,T_{n}$, from the set of $J$ behavioral categories. A Markov model for categorical time series data is defined by a transition probability matrix, $\boldsymbol{P}_{nt}$, describing the time-varying transition probabilities between observed states $s_{n,t-1}$ and $s_{nt}$ as follows:

[TABLE]

where $p_{nijt}\equiv P(s_{nt}=j|s_{n,t-1}=i)$, and the row probabilities sum to one, $\sum\limits_{j=1}^{J}p_{nijt}=1$. Therefore, given the behavioral state at the previous time step ($t-1$), the current observed category (at time $t$) is modeled as a multinomial trial with probabilities from the corresponding row of the transition matrix (1). We let $y_{nijt}$ be an indicator for individual $n$’s transition from state $i$ at time $t-1$ to state $j$ at time $t$, in other words, $y_{nijt}=1(s_{nt}=j|s_{n,t-1}=i)$ is defined to be $1$ for the case when the $n^{th}$ individual is in state $i$ at time $t-1$ and state $j$ at time $t$, and [math] otherwise. Then the aggregated transition indices vector $\boldsymbol{y}_{nit}=\begin{bmatrix}y_{ni1t}&\ldots&y_{niJt}\end{bmatrix}^{\prime}$, along with the corresponding transition probabilities, $\boldsymbol{p}_{nit}=\begin{bmatrix}p_{ni1t}&\cdots&p_{niJt}\\ \end{bmatrix}^{\prime}$, describe a multinomial trial.

Covariates are introduced by the multinomial logistic link function on the elements of the vector of transition probabilities, $\boldsymbol{p}_{nit}$, (e.g., see Sung et al., 2007; Holsclaw et al., 2017):

[TABLE]

where we assume the vector of parameter coefficients $\boldsymbol{\beta}_{iJ}=0$ for all $i=1,2,...,J$ for identifiability. In our application, $J$ refers to the walking behavior. Therefore, coefficients for transitions to walking are 0 and do not appear in subsequent figures 2 and 3. Note that covariates may vary by individual and/or time, which induces non-homogeneous transition probabilities.

3.2.2 Parameter Model

The $B$-dimensional parameter vectors in the data model multinomial logistic functions, $\boldsymbol{\beta}_{ij}$, for $i=1,2,...,J$ and $j=1,2,...,J-1$ are partitioned into three components:

[TABLE]

where $\boldsymbol{\alpha}_{ij}$ is an $(N-1)$-dimensional vector of random individual effects, $\boldsymbol{\zeta}_{ij}$ is an $H$-dimensional vector of habitat intercepts and $\boldsymbol{\theta}_{ij}$ are $B-(N-1+H)$-dimensional vectors of fixed quantitative covariate effects. The random individual effects are subject to a “sum to zero” constraint and therefore we only sample $(N-1)$ coefficients. We assume independent normal distribution priors for each of the sampled coefficients as given below. The habitat coefficients for each transition are assumed to have the same mean, $\mu_{ij}$. This mean is equivalent to the average intercept and, therefore, each habitat coefficient can be interpreted as the mean plus the habitat effect, similar to a cell means model in ANOVA. The common mean is assigned a flat prior. The hierarchical centering of these habitat parameters leads to less correlation between parameters compared to estimating $H$ habitat effects and an intercept in the regression model (Gilks et al., 1995). In addition, centering provides more interpretability compared to including a baseline intercept and $H-1$ coefficients. In summary, the prior distributions for these coefficients are given by

[TABLE]

where $\mathbf{0}$ and $\mathbf{1}$ are vectors of zeroes and ones, respectively. Note, all prior variances ($\sigma_{\alpha_{ij}}^{2},\sigma_{\zeta_{ij}}^{2},\sigma_{\theta_{ij}}^{2}$) are fixed at 100 to induce a vague prior and reduce number of parameters needing to be estimated.

3.3 Model Fitting

The likelihood of each coefficient vector, $\boldsymbol{\beta}_{ij}$, is the product of multinomial logistic functions. As the dimension of $\boldsymbol{\beta}_{ij}$ increases, the tuning of a Metropolis-Hastings algorithm becomes increasingly difficult. Latent variable schemes provide ways to automatically sample from the posterior without tuning and are generalized from binomial logistic regression to multinomial logistic regression (Albert and Chib, 1993; Holmes and Held, 2006; Frühwirth-Schnatter and Frühwirth, 2007; Polson et al., 2013). Specifically, the latent variable schemes express the likelihood as a product of binary logistic functions by the following transformation:

[TABLE]

where $\eta_{nijt}=\psi_{nijt}-C_{nijt}$ and $C_{nijt}=log\sum_{k\neq j}exp\{\psi_{nikt}\}$. In this formulation, $\eta_{nijt}$ is the log odds for a binomial random variable, which indicates whether or not the transition was from state $i$ to $j$ for individual $n$ at time point $t$ (Holmes and Held, 2006; Polson et al., 2013). From Theorem 1 of Polson et al. (2013), the product of logistic functions in the likelihood (5) is proportional to the product of of an exponential of $\eta_{nijt}$ and a Pólya-Gamma kernel. Therefore, by introducing the Pólya-Gamma latent variables, the sampling is done by conditional Gibbs updates alternating between the coefficients, $\boldsymbol{\beta}_{ij}$, and Pólya-Gamma latent variables. The details for the derivation of the full conditionals are in supplementary material (see supplementary materials S.3).

3.3.1 Multiple Imputation

The uncertainty associated with classifying the ACC data into categorical data by random forest is incorporated in the posterior distribution of the parameters via multiple imputation. Multiple imputation is used in situations involving missing data (e.g., the true behavior category) to provide approximate inference for the parameters based on the observed data (e.g., ACC fixes) (Rubin, 2004; Scharf et al., 2017; McClintock, 2017).

Let $\boldsymbol{\beta}$ be the collection of regression coefficients for all transitions, $\mathcal{S}$ be the collection of the state transitions for all individuals, and $\mathcal{A}$ be the collection of ACC data for all individuals. Our target distribution is the posterior of the unknown parameters, $\boldsymbol{\beta}$, given the behavior labels, $\mathcal{S}$. However, we do not observe the true behavior labels directly, but we predicted the behavior labels by classifying the observed auxiliary data, $\mathcal{A}$. Therefore, our target distribution is

[TABLE]

where given the behavior labels, $\mathcal{S}$, the parameters are conditionally independent of the ACC data. We assumed the distribution of the behavior labels given the ACC data, $\left[\mathcal{S}|\mathcal{A}\right]$, is the prediction from the supervised classification random forest.

We follow the multiple imputation MCMC algorithm outlined by Scharf et al. (2017) which numerically marginalizes over $\left[\mathcal{S}|\mathcal{A}\right]$ by randomly sampling from a set of $M$ potential data sets where $M=200$:

Draw realizations from the imputation distribution, $\mathcal{S}^{(m)}_{1T_{1}},...,\mathcal{S}^{(m)}_{NT_{N}}$, for $m=1,...,M$. 2. 2.

For each iteration of the MCMC repeat:

• •

Select an imputation data set with probability $1/M$.

• •

Use the Gibbs updating steps in supplementary materials (S.3) with $(\mathcal{S}_{1T_{1}},...,\mathcal{S}_{NT_{N}})=(\mathcal{S}^{(m)}_{1T_{1}},...,\mathcal{S}^{(m)}_{NT_{N}})$.

Throughout model fitting, we used “walking” as our reference behavioral category, $J$, in the multinomial logistic function (2). The coefficients for the reference were set to [math] and therefore are not displayed in the results. We chose walking as the reference category because we did not have specific hypotheses for this behavior and there was little variation in acceleration due to walking among individuals. Our method is amenable to choosing a different reference category specific to the study species. If it is of interest to more easily interpret effects on state duration, the reference category would depend on the previous state (i.e., for transitions from flight, the reference category would be flight).

We assessed parameter convergence by monitoring trace plots and setting different random starting values (Brooks and Gelman, 1998). For inference, we sampled 15000 iterations from the model posterior and the first 5000 were discarded as burn in to ensure summaries were not influenced by the starting values. We generated behavior sequences from the current values of the parameters each iteration to investigate the goodness of fit of the model (supplementary web material S.4).

The posteriors for odds ratios were obtained by exponentiating the coefficients each iteration. The posterior samples for coefficients and odds ratios were summarized to the posterior mean and 95% credible intervals. We determined significance by whether the 95% credible intervals included zero for the quantitative covariates for weather and time of day. The significance of the quantitative coefficients corresponded to the 95% credible interval for the odds ratio not containing one and a proportion of >0.95 (positive effect on transition probability) or <0.05 (negative effect on transition probability) of iterations in which the odds ratio was >1. For each transition probability, we investigated the differences in behavioral transition probabilities among habitats by calculating the pairwise proportion of iterations in which habitat coefficients differed in magnitude, $\zeta_{aij}>\zeta_{bij},\;a,b\in\{1,...,H\},\;and\;a\neq b$. Two habitat coefficients for a transition were considered significantly different if the proportion of iterations with a difference in magnitude was >0.95 (the probability of transition from behavior $i$ to $j$ is greater in habitat $a$ than habitat $b$) or <0.05 (the probability of transition from behavior $i$ to $j$ is lesser in habitat $a$ than habitat $b$). We did not adjust the proportion cut off or widen credible intervals to account for inherent multiplicity in this case because we modeled the habitat coefficients with a common mean, which pooled the coefficient estimates (4). The hierarchical centering shrinks estimates toward the common mean, $\mu_{ij}$, which makes it harder for significant pairwise differences to occur, thus eliminating the need to make additional post hoc adjustments for multiple comparisons (Gelman et al., 2012).

The assumption of a discrete time Markov process implicitly accommodates inference on the original sampling scale of the ACC data schedule (e.g., every 6 minutes). We interpreted all coefficient estimates as effects on behavioral transition probabilities from state $i$ to state $j$ at a 6 minute interval, relative to a base behavior state of walking. For brevity, we did not explicitly restate “relative to walking at a 6 minute interval" for each interpretation. Similarly, the odds ratios are interpreted as the multiplicative change in odds of transitions to state $j$ from state $i$ versus transitions to walking from state $i$.

4 Results

Our modeling framework establishes flexibility and efficiency in estimating covariate effects in the behavioral decision-making process by specifying transition-specific coefficients. The variability in coefficients across from states is indicative of the complexity in the decision-making process in migrating white-fronted geese. The estimated effects of different habitats on behavior transitions rates did not follow the same patterns across the from states (Fig. 2). By contrast, for weather and time of day covariates, the pattern in the coefficient estimates was similar across from states (Fig. 3).

There were significant differences between habitats for every behavior transition probability except flight to stationary (see supplementary materials S.5). We estimated a significantly higher probability of transitioning to feeding from feeding (i.e., continued to feed), stationary and walking for food crops such as corn and soybeans and a significantly lower probability when the birds were in open water habitat (Fig. 4). The positive coefficient estimates for corn and soybeans indicated an odds ratio >1 and an increase in the probability for transitioning to feeding from flight and feeding (Fig. 2 to feeding subplots). For example, the mean effect of corn on the odds of continuing to feed was 0.42 with corresponding mean odds ratio 1.52, indicating the odds of continuing to feed at the next time step was 52% greater than the odds of transitioning to walking from feeding in corn habitat (Table 1). Transitions to stationary from feeding and walking were more probable in open water habitats than in corn (see supplementary materials S.5). In general, the widths of the credible intervals in Figure 2 have an inverse relationship with the observed frequency of habitats (see supplementary materials S.1 for frequencies). For example, the most frequently used habitats were corn and open water, which tended to have the narrowest credible intervals.

The probability of remaining in flight (i.e., transitioning to flight from flight) increased significantly during the first half of the day as indicated by the left most credible interval for sin(time) in the top left subplot of Fig. 3(a). There were few ACC fixes classified as flight at night, which corresponds to the significant negative effects of cos(time) on all transition probabilities to flight (Fig. 3(a)). ACC fixes were classified as stationary most frequently and in greater proportions at night. Therefore, it is not surprising that many diurnal coefficients corresponding to transition probabilities to the stationary state were significant. Table 2 shows that 6 out of 8 diurnal coefficients for transitions to stationary were significant. The transition probabilities to feeding and to flight from any behavior decreased overnight (negative effect of cos(time) across all from states to feeding and flight; Fig. 3(a)).

Daily minimum and maximum temperatures did not affect the probability of remaining in flight the same way. Warmer than average daily minimum temperatures decreased the probability of remaining in flight (negative credible interval in top minimum temperature subplot of Fig. 3(b)), while the warmer than average daily maximum temperatures increased the probability (positive credible interval in top maximum temperature subplot of Fig. 3(b)). Increased wind speeds decreased the probability of remaining in flight and transitioning to feeding from feeding and walking. We found no evidence of the weather variables effecting transition probabilities from stationary behavior. The probability of remaining in flight was significantly affected by minimum daily temperature, maximum daily temperature and wind speed, but not wind direction (Table 2).

5 Discussion

We provide a unified framework to connect variation in animal behavior with variation in habitat use and weather by propagating uncertainty in transitions using multiple imputation within a Bayesian Markov model with data from a long-distance migratory bird. Our approach is broadly applicable to other focal species and study systems across ecology. By analyzing data on the scale of frequency of collection, inferences are more intuitive and appropriate than aggregating to proportions. Importantly, our approach allows analysis of the behavior sequence with inherent temporal dependence and inference about covariate effects on behavior transition probabilities.

Transition matrices are well studied and can provide a wealth of inference and prediction beyond what is presented in this study including simulation of behavior sequences in different settings. Further, transition matrix models are used in ecological and evolutionary research beyond behavioral applications and Markov models. We extended our Markov model with a Pólya-Gamma sampling scheme, which will be useful for fast and automatic estimation of other complex ecological models that utilize the logistic link function.

We also implemented a unique approach to analyzing both ACC and GPS data from tracking devices. When we know the specific location of animals on the landscape, quantifying the effects of habitat on behavior transitions provides unprecedented information regarding the differential rates of behavior in specific habitats. We made simple assumptions when we constructed our covariates by assigning values to ACC fixes from the most recent GPS fix, but there is the potential to link to existing animal movement models for interpolation of location data at times of ACC observations. McClintock et al. (2017) used a continuous time correlated random walk model to predict locations at a regular time interval (Johnson et al., 2008). The locations could be predicted to the time points of the ACC fixes and covariate values could reflect the prediction location or be imputed from the prediction distribution (Hooten et al., 2010; Scharf et al., 2017).

In the case of white-fronted geese during spring migration, different rates of behavior can be attributed to different habitats. Larger relative effects of habitats associated with food sources such as corn and soybeans on transition probabilities to feeding compared to effects of wetland habitats and open water aligns with previous knowledge of white-fronted goose ecology (Ely, 1992; Krapu et al., 1995). Open water consistently had lower transition probabilities than food habitats for transitions to feeding and flight. Our finding that white-fronted geese were less likely to transition to feeding or flight in open water compared to food habitats was consistent with our expectations. Also, we found significant effects of weather on behavior transition probabilities of white-fronted geese during spring migration. Higher winds decreased the odds of remaining feeding or transitioning to feeding from stationary or walking behavior. Duration of flight behavior (flight to flight transitions) was the most influenced by weather. The opposite effects of minimum and maximum daily temperature on flight durations may be indicative of more complex decision-making processes by white-fronted geese. An increase in flight duration with an increase in maximum daily temperature aligns with our a priori assumptions that these birds do not often migrate beyond the snow line during spring because food is relatively inaccessible under snow.

Furthermore, we were able to control for diurnal patterns in activity by using a continuous transformation rather than discretizing the time of day. The effects of time of day were consistent across from states. Specifically, the coefficients for cos(time) all appeared to be significant in the same direction and magnitude in Figure 3(a). In the future, it may be beneficial to relax the assumption that all covariate effects are transition specific and instead estimate a mix of transition-specific and behavior-specific coefficients. For example, Holsclaw et al. (2017) estimated a transition matrix for a Hidden Markov model with a transition-specific intercept and state specific coefficients for all covariates. They estimated weather covariates specific to rainfall states at weather stations in India to identify global effects of weather systems. For animal migration and behavior, the occurrence of certain weather systems or the time of day may always result in a specific decision regardless of the current behavior which suggests a model with behavior specific coefficients may be a better fit. Within the Bayesian framework, models with different parameter formulations can be compared using tools such as Bayes factors, information criteria, or prediction of hold out samples to test hypotheses about the behavior process (Hooten and Hobbs, 2015).

Our analysis was limited to a month-long subset of the spring 2018 migration and this pattern may become more clear with the inclusion of more years. Although much of the inference was verification of previous knowledge about white-fronted geese, the methodology allows us to infer about habitat use and behavior simultaneously using both the GPS and ACC information. Most importantly, we developed a detailed picture of time allocation in reference to the specific habitat types used by white-fronted geese during spring migration which is not addressed by traditional models used in previous activity budget analysis, movement trajectory prediction, or resource selection frameworks.

There is diminishing return on estimation and inference after a sufficiently large number of imputation data sets are used (Scharf et al., 2017; McClintock, 2017). In order to investigate our choice of 200 data sets, we compared posterior inference among 200 imputation data sets, 100 imputation data sets and 1 data set corresponding to the most likely behavior classification. There was a general consensus among the three scenarios on the direction of the effects. Compared to the estimates based on the most likely classification, the estimates for the imputation scenarios were shrunk towards zero. Credible interval widths for the imputation scenarios appeared consistent or larger than the widths for coefficients estimated from the most likely behaviors (see supplementary web material (S.6) for side-by-side comparisons).

Advances in animal tracking technologies continue to provide more frequently collected data for a greater duration of time. Thus, rich data sets are emerging as never before for ecologists and evolutionary biologists. Hence, there is an increasing need for development of models that appropriately handle the structure and volume of collected information for improved inference. Our Markov model framework provides much more capability for directly interpreting behavior patterns. In addition, the Pólya-Gamma latent variables facilitate for more efficient sampling and have yet to be used in the animal behavior and movement literature. If classification of behaviors is not feasible, the Pólya-Gamma sampling scheme can be incorporated into Bayesian estimation of transition probability matrices in a hidden Markov model framework (Holsclaw et al., 2017). The model directly handles temporal dependence in ACC data and learns about the behavior process from both ACC and GPS data. Our results suggest that new data sources coupled with appropriate modeling have unprecedented potential to provide a comprehensive understanding of complex ecological and evolutionary processes in animal movement.