Imputation Approaches for Animal Movement Modeling

TL;DR

This paper reviews imputation methods for animal movement modeling, evaluates their accuracy through simulations, and demonstrates their application on telemetry data from northern fur seals, highlighting their strengths and limitations.

Contribution

It introduces imputation approaches for animal movement models, assesses their performance via simulations, and applies them to real telemetry data, filling a gap in understanding their effectiveness.

Findings

Location-related parameters are estimated more accurately than higher-order movement parameters.

Imputation methods can be useful but vary in accuracy depending on the parameter of interest.

Simulation results highlight the importance of model choice in movement inference.

Abstract

The analysis of telemetry data is common in animal ecological studies. While the collection of telemetry data for individual animals has improved dramatically, the methods to properly account for inherent uncertainties (e.g., measurement error, dependence, barriers to movement) have lagged behind. Still, many new statistical approaches have been developed to infer unknown quantities affecting animal movement or predict movement based on telemetry data. Hierarchical statistical models are useful to account for some of the aforementioned uncertainties, as well as provide population-level inference, but they often come with an increased computational burden. For certain types of statistical models, it is straightforward to provide inference if the latent true animal trajectory is known, but challenging otherwise. In these cases, approaches related to multiple imputation have been employed to account for the uncertainty associated with our knowledge of the latent trajectory. Despite the increasing use of imputation approaches for modeling animal movement, the general sensitivity and accuracy of these methods have not been explored in detail. We provide an introduction to animal movement modeling and describe how imputation approaches may be helpful for certain types of models. We also assess the performance of imputation approaches in a simulation study. Our simulation study suggests that inference for model parameters directly related to the location of an individual may be more accurate than inference for parameters associated with higher-order processes such as velocity or acceleration. Finally, we apply these methods to analyze a telemetry data set involving northern fur seals (Callorhinus ursinus) in the Bering Sea.

Figures

Tables

Table 1. Table 1: Posterior medians, 95% credible intervals, and prior distributions for the application to northern fur seal movement.

	posterior		prior
parameter	median	(2.5%, 97.5%)	density
${\sigma }_s$ (kms)	1.54	(1.50, 1.58)	$\text{Unif}(0,\infty )$
${\sigma }_v$ (hour-1)	0.269	(0.262, 0.275)	$\text{Unif}(0,\infty )$

Equations

$$\left[\mathbf{\theta }|\mathbf{s}\right]$$

$$=\int \left[\mathbf{\theta }|\mathbf{s},{\mathbf{s}}_m\right]\left[{\mathbf{s}}_m|\mathbf{s}\right]d{\mathbf{s}}_m,$$

$$\mathbb{E}\left(\mathbf{\theta }|\mathbf{s}\right)=\mathbb{E}\left(\mathbb{E}\left(\mathbf{\theta }|\mathbf{s},{\mathbf{s}}_m\right)|\mathbf{s}\right),$$

$$\text{Var}\left(\mathbf{\theta }|\mathbf{s}\right)=\mathbb{E}\left(\text{Var}\left(\mathbf{\theta }|\mathbf{s},{\mathbf{s}}_m\right)|\mathbf{s}\right)+\text{Var}\left(\mathbb{E}\left(\mathbf{\theta }|\mathbf{s},{\mathbf{s}}_m\right)|\mathbf{s}\right).$$

$$\begin{array}{rl}d\mathbf{\mu }(t)& =\mathbf{v}(t)dt\\ d\mathbf{v}(t)& =-\nabla H\left(\mathbf{\mu }(t),t\right)dt-{\sigma }_v\mathbf{v}(t)dt+{\sigma }_{v}d\mathbf{b}(t),\end{array}$$

$$H(\mathbf{\mu }(t),t)=\beta (t)\Vert \mathbf{\mu }(t)-\mathbf{c}{\Vert }_2,$$

$$\mathbf{s}(t)|\mathbf{\mu }(t)$$

$${\sigma }_s^2$$

$${\sigma }_v^2$$

$$\mathbf{\beta }$$

$$\mathbf{\alpha }$$

$$\left[\mathbf{\theta }|\mathbf{s}\right]$$

$$d\mathbf{\mu }(t)$$

$$\mathbf{\mu }(0)|{\sigma }_0^2$$

$$\left[\mathbf{\mu }(t_i)|\mathbf{\mu }(t_{i-1}),\beta \right]$$

$$=\text{N}(\mathbf{\mu }(t_{i-1})+\beta \frac{\mathbf{c}-\mathbf{\mu }(t_{i-1})}{\Vert \mathbf{c}-\mathbf{\mu }(t_{i-1}){\Vert }_2}\Delta t_i,\Delta t_i{\mathbf{I}}_2)$$

$$\left[\mathbf{s}(t_i)|\mathbf{\mu }(t_i),{\sigma }_s^2\right]=\text{N}\left(\mathbf{\mu }(t_i),{\sigma }_s^2{\mathbf{I}}_2\right).$$

$$\left[\beta |{\sigma }_{\beta }^2\right]$$

$$\left[{\sigma }_s^2|a_s,b_s\right]$$

$${\sigma }_0^2$$

$$\left[\mathbf{\mu }(t_1=0)|\cdot \right]$$

$$\begin{array}{rl}\left[\mathbf{\mu }(t_i)|\cdot \right]& \propto \left[\mathbf{\mu }(t_{i+1})|\mathbf{\mu }(t_i),\beta \right]\left[\mathbf{\mu }(t_i)|\mathbf{\mu }(t_{i-1})\right]\\ & \times {\left[\mathbf{s}(t_i)|\mathbf{\mu }(t_i),{\sigma }_s^2\right]}^{1_{\left\{\text{observation occurred at }{t}_i\right\}}}\end{array}1<i<m$$

$$\left[\mathbf{\mu }(t_m)|\cdot \right]$$

$$\left[\beta |\cdot \right]$$

$$\left[{\sigma }_s^2|\cdot \right]$$

Full text

Keywords: animal movement models, hierarchical models, telemetry data, multiple imputation

1 Introduction

An increasing diversity of telemetry data types and their use in a variety of scientific studies has resulted in rapid development of models for animal trajectories (McClintock et al., 2014; Kays et al., 2015). Statistical models for animal movement processes mostly assume one of three forms: point process models, discrete-time models, and continuous-time models (Hooten et al., 2017). Each of these three types of statistical models are aimed at answering certain characteristic questions. For example, point process models are often used to infer individual- and/or population-level “resource selection,” (e.g., Brost et al., 2015), whereas discrete-time models are commonly used to temporally cluster trajectories into discrete components that may relate to animal behavior (e.g., Langrock et al., 2012; McClintock et al., 2013). Continuous-time models for animal movement have linked trajectories to mechanistic stochastic processes (e.g., Brillinger, 2010; Hooten and Johnson, In Press) and are also used to interpolate telemetry data; that is, to predict an individual’s position at times of interest (Johnson et al., 2008; Hooten et al., 2017).

Hierarchical approaches to model specification facilitate population-level inference (Hooten et al., 2016) and account for telemetry measurement error (e.g., Brost et al., 2015; Buderman et al., 2016). Some forms of telemetry measurement error have unique features due to the data collection procedure (e.g., the satellite orbit of Service Argos; McClintock et al. 2015) that must be explicitly accounted for in statistical models for animal movement. Hierarchical models provide a rigorous framework that both allows for feedback between the data- and process-level components of a complete movement model, and properly accounts for uncertainty due to measurement error.

While the advantages of hierarchical models are clear, they can often be challenging to fit to real data sets due to their inherent complexities. In some cases, a full hierarchical model may be impossible to implement when the process is specified for a functional of the true trajectory that is insufficient (i.e., carries less information content than the position process). If an animal movement model is specified for a certain type of latent process that does not reconcile well with the telemetry data, then traditional hierarchical modeling approaches may not be feasible. For example, the discrete-space models for movement in Hooten et al. (2010), Hanks et al. (2015), and Hanks and Hughes (2016) involve a transformation of the continuous underlying path into residence times associated with a finite set of areal regions. It is straightforward to calculate the areal residence times based on the true trajectory, but impossible to recover the true trajectory based on the areal residence times alone. This noninvertible functional of the trajectory precludes the ability to fit the full hierarchical model using standard methods (e.g., maximum likelihood or Bayesian algorithms; Hooten et al. 2010).

In other settings, it may be possible to specify and fit a fully hierarchical model, but doing so would involve excessive computational resources (e.g., Scharf et al., 2016). For example, a transformation of the trajectory into one or more of its temporal derivatives results in a consistent representation of the trajectory (if the initial conditions are known); however, if the temporal derivative itself is treated as a latent stochastic process in a hierarchical model, we may not be able to analytically marginalize over the true trajectory to obtain an integrated likelihood. Thus, an algorithm to fit the model, such as Markov chain Monte Carlo (MCMC), must be constructed to sample the latent continuous-time trajectory in addition to the other model parameters. It may not be feasible to quickly obtain samples of the true latent trajectory process in such a model, thus, alternative approaches based on imputation have been considered.

Figure 1 shows the observed locations of a northern fur seal (Callorhinus ursinus) as it makes several foraging excursions from the Pribilof islands in Alaska. By imputing several samples from an approximation to the distribution of the true, continuous-time movement process, it is possible to obtain approximate inference from a movement model that would otherwise be too computationally expensive to fit (Section 4).

Imputation strategies have been used in statistics for decades to address missing data issues. For example, a simple, yet ad hoc, imputation procedure for telemetry data involves reconstructing the trajectory using linear interpolation (i.e., connecting the dots). The imputed trajectory is then treated as data in a statistical model that requires knowledge of the true unknown continuous path. This procedure is usually inadequate for several reasons: (1) the linearly interpolated path is not based on an understanding of the dynamics in the trajectory and, thus, provides a biased understanding of the individual’s position between observed time points, (2) the telemetry data are considered to be perfectly accurate observations of the trajectory at the observation times, and (3) the uncertainty in our understanding of the trajectory between (or at) observation times is not accommodated.

Multiple imputation for random trajectories was developed to help alleviate some of these issues (Rubin, 1996; Hooten et al., 2010; Hanks et al., 2011). By imputing several, model-based paths, it is possible to account for underlying movement dynamics and measurement error. In what follows, we review imputation strategies from a Bayesian perspective. We then show how to use imputation for fitting continuous-time continuous-animal space movement models (see Scharf et al., 2016; McClintock, 2017, for the discrete-time setting). We examine the strengths and weaknesses of imputation approaches empirically, using simulation. Finally, we demonstrate imputation while modeling real telemetry data for a northern fur seal.

2 Methods

2.1 Multiple imputation

Conventionally, statistical methods for imputation were developed to address issues with missing data. Prior to more formal multiple imputation approaches (e.g., Rubin, 1996), a suite of ad hoc methods were proposed (e.g., complete-case analysis, available-case analysis, fill-in with means; Little and Rubin 1987). In the Bayesian context, the basic idea behind imputation is that inference for the model parameters $\boldsymbol{\theta}$ under the observed data can be obtained using the posterior

[TABLE]

where, $\mathbf{s}_{m}$ represents the missing data, and, following Gelfand and Smith (1990), we use brackets to denote a probability distribution. The distribution $\left[\boldsymbol{\theta}|\mathbf{s},\mathbf{s}_{m}\right]$ is often referred to as the “complete-data” posterior and the distribution $\left[\mathbf{s}_{m}|\mathbf{s}\right]$ is the distribution of the missing data, conditioned only on the observed data (i.e., the imputation distribution). Thus, heuristically, the desired posterior distribution $\left[\boldsymbol{\theta}|\mathbf{s}\right]$ is the complete-data posterior distribution averaged over the imputation distribution (Rubin, 2004).

The Bayesian imputation concept yields useful properties of posterior moments. In particular, the posterior mean and variance of the parameters can be written using the iterated expectation and variance formulations:

[TABLE]

and

[TABLE]

In practice, for a given observed data set, if $K$ samples of the missing data can be obtained from the imputation distribution, then the posterior expectation can be approximated by fitting the complete-data models to each complete-data set and then averaging the posterior means. Similarly, by averaging the posterior variances from the complete-data models and adding the sample variance of the posterior means, we arrive at an approximation of the correct posterior variance $\text{Var}\left(\boldsymbol{\theta}|\mathbf{s}\right)$. These approximations improve as the number of missing data samples from the imputation distribution ($K$) increases.

Multiple imputation depends on the ability to evaluate the complete-data posterior, and the ability to sample missing data sets from the imputation distribution. In the context of animal movement, telemetry data are measurements of the true positions with error and models for movement dynamics are most commonly specified for the true position process. By treating the true position process as missing data and taking a Bayesian multiple imputation approach, it is possible to obtain approximate inference under models that we would otherwise be unable to use.

2.2 Animal movement

To introduce the concept of imputation in a continuous-time animal movement setting, we denote the set of telemetry data as $\mathbf{s}\equiv(\mathbf{s}(t_{1})^{\prime},\ldots,\mathbf{s}(t_{n})^{\prime})^{\prime}$, a $2n\times 1$ vector, for observation times $t_{1},\ldots,t_{n}$. In this setting, each $\mathbf{s}(t_{i})$ is a $2\times 1$ vector representing an observed telemetry position in two-dimensional space at time $t_{i}$. Similarly, we represent the true latent position process as $\boldsymbol{\mu}\equiv\left(\boldsymbol{\mu}(t_{1})^{\prime},\ldots,\boldsymbol{\mu}(t_{m})^{\prime}\right)^{\prime}$, a $2m\times 1$ vector where each $\boldsymbol{\mu}(t_{j})$ represents a true, but unknown, position at time $t_{j}$.

The true position process is typically modeled as arising from a distribution conditioned on a set of process model parameters $\boldsymbol{\theta}$, such that $\boldsymbol{\mu}\sim\left[\boldsymbol{\mu}|\boldsymbol{\theta}\right]$. In a hierarchical framework, the telemetry data are often modeled as conditional on the true latent trajectory and observation parameters, $\boldsymbol{\psi}$, as $\mathbf{s}\sim\left[\mathbf{s}|\boldsymbol{\mu},\boldsymbol{\psi}\right]$. The hierarchical framework is useful for analyzing telemetry data when it is reasonable to assume that the processes that gives rise to measurement error is conditionally independent from the processes that describe the behavior of the animal.

An important class of models for animal movement that fits into our hierarchical framework can be described using sets of stochastic differential equations (SDEs). Brillinger and Stewart (1998) introduced one such model to study the migration of northern elephant seals that explicitly accounted for the curvature of the Earth. A general model that describes the motion of a Brownian particle in an external, conservative force field is given by the SDEs

[TABLE]

where $d\mathbf{b}(t)\sim\text{N}\left(\mathbf{0},\;\mathbf{I}_{2}dt\right)$ and $H(\cdot,\cdot)$ is a time-varying potential function, the negative gradient of which yields the force acting on the Brownian particle (Nelson, 1967). In practice, a dense grid of times $t_{1},\dots,t_{m}$ are used to fit the model, and each location, $\boldsymbol{\mu}(t_{j})$, and velocity, $\mathbf{v}(t_{j})$, have an associated $dt_{j}=t_{j}-t_{j-1}$.

Modeling animal movement using (5) allows researchers to study the impacts of a wide variety of natural phenomena on an individual’s behavior by specifying an appropriate potential function. For instance, if we define $H$ as

[TABLE]

where $\|\boldsymbol{\mu}(t)-\mathbf{c}\|_{2}$ is the Euclidean distance between $\boldsymbol{\mu}(t)$ and $\mathbf{c}$, and $\mathbf{c}$ is a feature of interest in the environment, then the force on the particle is a pull toward $\mathbf{c}$ for $\beta(t)>0$, and repulsion from $\mathbf{c}$ for $\beta(t)<0$. The model in (5) – (6) provides an excellent description of the behavior of many marine mammals, such as northern fur seals. In summer months, female northern fur seals spend several days at sea foraging before returning to a specific rookery to nurse their pups (Johnson et al., 2013). By defining $\mathbf{c}$ to be the location of the rookery, we can accommodate the time-varying effect it has on a nursing female’s movement (see Section 4).

To form a complete hierarchical model, we specify independent, circular Gaussian measurement error, and conjugate priors for $\sigma_{s}^{2}$ and the process parameters in $\boldsymbol{\theta}$ as

[TABLE]

where $\mathbf{W}$ is a matrix of basis functions used to model the dynamic process $\boldsymbol{\beta}=\left(\beta(t_{1}),\dots,\beta(t_{m})\right)^{\prime}$ (e.g., splines, Hefley et al. 2017). Fitting the SDE-based model to telemetry data using an MCMC algorithm is computationally intensive due to the non-linear relationship between $\boldsymbol{\mu}$ and $\boldsymbol{\theta}\equiv\left(\boldsymbol{\beta}^{\prime},\sigma_{v}^{2}\right)^{\prime}$. Standard methods require an MCMC algorithm in which $\boldsymbol{\mu}(t)$ must be updated over a fine grid of times, and for many modern telemetry data sets, the effort required to fit the model can exceed the computational resources available on even the most advanced high performance computing platforms. However, conditioned on the true path, $\boldsymbol{\mu}$, inference for $\boldsymbol{\theta}$ is straightforward.

2.3 Process Imputation

Hooten et al. (2010) and Hanks et al. (2011) proposed an approximate-model fitting procedure that is similar in spirit to multiple imputation, which we term “process imputation.” The motivation for the process imputation approach arises from a decomposition of the posterior distribution of the process parameters. Using standard properties of conditional probability, we write the posterior distribution as

[TABLE]

The form of the integral in (12) is similar to that used in conventional data imputation approaches (2); however, in this case, the process $\boldsymbol{\mu}$ takes the place of the missing data $\mathbf{s}_{m}$.

As with conventional imputation approaches, if we can obtain the posterior distribution $\left[\boldsymbol{\theta}|\boldsymbol{\mu},\mathbf{s}\right]$ and samples from the process imputation distribution $\left[\boldsymbol{\mu}|\mathbf{s}\right]$, then we can take the same approach of fitting the process model to each of the imputed process samples and approximating posterior moments for $\boldsymbol{\theta}$ using the relevant iterated expectation and variance equations.

For many movement models, including both the SDE-based model introduced in (5), and the discrete-space models presented by Hooten et al. (2010), Hanks et al. (2015), and Hanks and Hughes (2016), evaluating $\left[\boldsymbol{\theta}|\boldsymbol{\mu},\mathbf{s}\right]$ up to proportionality is straightforward. The process model parameters, $\boldsymbol{\theta}$, are often conditionally independent of data, $\mathbf{s}$, and the observation parameters, $\boldsymbol{\psi}$, such that the complete-data posterior simplifies to $\left[\boldsymbol{\theta}|\boldsymbol{\mu}\right]\propto\left[\boldsymbol{\mu}|\boldsymbol{\theta}\right]\left[\boldsymbol{\theta}\right]$. In practice, the challenge to implementing a process imputation approach is in sampling from the imputation distribution, $\left[\boldsymbol{\mu}|\mathbf{s}\right]$. The difficulty arises for the same reasons we are unable to fit the full hierarchical model. Sampling from the imputation distribution, $\left[\boldsymbol{\mu}|\mathbf{s}\right]\propto\int\left[\mathbf{s}|\boldsymbol{\mu},\boldsymbol{\psi}\right]\left[\boldsymbol{\mu}\right]\left[\boldsymbol{\psi}\right]d\boldsymbol{\psi}$, requires that we have access to the marginal distribution $\left[\boldsymbol{\mu}\right]$, which, in turn, requires we evaluate the integral $\left[\boldsymbol{\mu}\right]=\int\left[\boldsymbol{\mu}|\boldsymbol{\theta}\right]\left[\boldsymbol{\theta}\right]d\boldsymbol{\theta}$. For the models we consider, this integral is intractable, either for reasons involving noninvertibility or computation.

In what follows, we assume that direct sampling from $\left[\boldsymbol{\mu}|\mathbf{s}\right]$ is impossible. To circumvent this obstacle, we assume the existence of another conditional random variable $\boldsymbol{\mu}^{*}$. The distribution of $\boldsymbol{\mu}^{*}$ is sufficiently similar to $\boldsymbol{\mu}$ that draws from the distribution $\left[\boldsymbol{\mu}^{*}|\mathbf{s}\right]$, which may or may not depend on $\boldsymbol{\psi}$, can be used to approximately evaluate the integral in (12). We term the distribution given by $\left[\boldsymbol{\mu}^{*}|\mathbf{s}\right]$ the “approximate imputation distribution” (AID), and denote the resulting approximation to the true posterior by $\left[\boldsymbol{\theta}^{*}|\mathbf{s}\right]$. When implementing process imputation, selection of the AID is performed by first specifying a model for $\left[\boldsymbol{\mu}^{*}|\mathbf{s},\boldsymbol{\phi}\right]$ parameterized by $\boldsymbol{\phi}$. Estimates for the parameters $\boldsymbol{\phi}$ are obtained in a preliminary stage by fitting $\left[\boldsymbol{\mu}^{*}|\mathbf{s},\boldsymbol{\phi}\right]$ to the data. Common models for the AID used in the literature include the velocity-based Ornstein-Uhlenbeck model of Johnson et al. (2008) and the spline-based functional model of Buderman et al. (2016).

As motivation for substituting draws for the true imputation distribution with draws from an approximate distribution, we offer the following heuristic argument based on asymptotically dense observations with negligible measurement error. Additionally, we discuss the results of two simulation studies in Section 3 that provide practical guidance for using a process imputation approach.

When observations with small measurement error are available, the distribution $\left[\boldsymbol{\mu}|\mathbf{s}\right]$ has probability mass concentrated on paths passing through the immediate vicinity of the telemetry locations, $\mathbf{s}$. For the extreme case in which an observation $\mathbf{s}(t)$ exists for all possible times $t$, and measurements are made without error, we have $\mathbf{s}(t)=\boldsymbol{\mu}(t)$ for all $t$, and the conditional distribution $\left[\boldsymbol{\mu}|\mathbf{s}\right]$ collapses to a point mass on the true path. At this limit, the distribution $\left[\boldsymbol{\mu}\right]$ could, in principle, be replaced by any probability measure with the correct support without affecting $\left[\boldsymbol{\mu}|\mathbf{s}\right]$. Thus, in the limiting case of infinitely dense observation times and zero measurement error, any AID that concentrates probability mass around paths through the data is sufficient. Following this argument, it is reasonable to expect that the degree to which substituting $\left[\boldsymbol{\mu}^{*}|\mathbf{s}\right]$ in equation (12) is an adequate approximation will depend on the density of observation times and the severity of measurement error.

An analogy to variational approximations (Ormerod and Wand, 2010) provides additional motivation for the process imputation approach. Variational approximations are often employed in a Bayesian setting to find approximations to posterior densities when sampling from the true densities is computationally infeasible. Optimal approximations are derived by minimizing the Kullback-Leibler divergence between a family of proposal distributions (e.g., distributions for which certain model parameters are assumed to be independent) and the true distribution. In the process imputation setting, the model for the AID takes on the role of the family of proposal distributions. Instead of minimizing the Kullback-Leibler divergence, we select the optimal approximation by fitting the AID model to the data. The accuracy of a variational approximation to the true posterior density depends, in part, upon the suitability of the family of proposal distributions. Our simulation study suggests that the choice of AID model in the process imputation setting has an analogous relationship with the approximate posterior distribution, $\left[\boldsymbol{\theta}^{*}|\mathbf{s}\right]$.

In practice, when measurement errors are non-zero and the number of observations is finite, the model for the AID can have a significant impact on inference. When implementing this procedure, a finite number of draws from the AID are used to compute the integral in (12); using too few draws can fail to adequately represent the AID, biasing inference. We investigated the effects of each of these factors on inference obtained for the SDE-based movement model using process imputation with a simulation study (Section 3).

3 Simulation studies

We conducted two simulation studies of the process imputation approach in the context of SDE-based models for movement. The results revealed important information about when the approach may be useful, and offer insight into the viability of process imputation for other movement models.

3.1 Study 1: Brownian particle in an external force field

The first simulation study was designed to explore the impacts of four different factors on the validity of inference obtained for the movement model specified in (5) – (11) using process imputation: 1) the severity of measurement error, 2) the temporal density of observations, 3) the model for the AID, and 4) the number of draws, $K$, from the AID used in the analysis. Figure 2 shows examples of simulated data (black points) along with imputed paths (gray curves) for the different combinations of measurement error severity and observation frequency. For all four combinations of measurement error (“large” and “small”) and observation frequency (“sparse” and “dense”), we simulated 24 paths from the SDE-based movement model (see Appendix A for details). We then fit two different AID models to the simulated data: the integrated Ornstein-Uhlenbeck model of Johnson et al. (2008) using the R package crawl (“OU”), and a Gaussian process (e.g., Fleming et al., 2015; Hooten and Johnson, In Press) with Gaussian covariance function (“GP”) (see Appendix A for further details). Finally, for each combination of measurement error, observation frequency, and AID model, we drew $K$ realizations from the AID to approximate the integral (12). We investigated three different values for $K$ (8, 32, and 128), and also implemented the process imputation approach using the posterior mean of the AID. To provide a baseline against which to compare the performance of process imputation, we also fit the SDE-based model conditioned on the true, underlying path $\boldsymbol{\mu}$.

3.1.1 Implementation

After a suitable AID has been identified, the process imputation approach can be implemented using an MCMC algorithm. We obtained samples from the approximation to the true posterior distribution, $\left[\boldsymbol{\theta}^{*}|\mathbf{s}\right]$, by drawing realizations form the complete-data posterior, $\left[\boldsymbol{\theta}|\boldsymbol{\mu}\right]$, and integrating with the respect to the AID. The two steps can be combined into a single MCMC algorithm using the following procedure:

Draw realizations from the AID, ${\boldsymbol{\mu}^{*}}^{(k)}\sim\left[\boldsymbol{\mu}^{*}|\mathbf{s}\right]$, for $k=1,\dots,K$. 2. 2.

MCMC procedure:

1. (a)

Randomly select one of the $K$ samples, ${\boldsymbol{\mu}^{*}}^{(k)}$, with probability $1/K$. 2. (b)

Update model parameters $\boldsymbol{\theta}|\boldsymbol{\mu}={\boldsymbol{\mu}^{*}}^{(k)}$ conditioned on the imputed path, ${\boldsymbol{\mu}^{*}}^{(k)}$. 3. (c)

Repeat steps 2.(a) - 2.(b) at each iteration of the MCMC algorithm.

3.1.2 Results

To evaluate the performance of the process imputation approach at estimating the effect of the potential function, we compared the posterior distributions of $\boldsymbol{\beta}$ to the values used to simulate the data through two summary measures. We assessed the accuracy of process imputation by computing the proportion of the study period for which the pointwise, equal-tailed, 95% credible interval contains the true value of $\beta(t)$, which we term the “coverage” region. We assessed the precision of process imputation by computing the proportion of the study period during which the pointwise credible interval both contains the truth, and does not contain 0, which we term the “detection” region. Thus, evidence of strong performance by the process imputation approach is characterized by a combination of large coverage and detection regions. Figure 3 shows an example of coverage and detection regions. For further details concerning coverage and detection regions, see Appendix A. To evaluate the performance of the process imputation approach at estimating $\sigma_{v}^{2}$ and $\sigma_{s}^{2}$, we computed the proportion of simulations in which equal-tailed 95% credible intervals contained the true values of the parameters.

Figure 4 shows the sizes of coverage and detection regions under each process imputation regime, repectively. The vertical range of the bars shows the 12.5% to 87.5% quantiles across each set of 24 simulations. Lighter shades correspond to cases for which the AID model is the Ornstein-Uhlenbeck model of Johnson et al. (2008), and darker shades correspond to cases for which the AID model is a Gaussian process. Rows and columns organize cases by measurement error and observation frequency, respectively. Finally, within each plot, cases are organized by the number of draws used from the imputation distribution in fitting the approximate hierarchical model. Figure 5 shows the proportion of the 24 simulated paths for which the equal-tailed 95% credible interval contained the true value of $\sigma_{v}^{2}$, organized analogously. The far left of each subplot in Figures 4 and 5 corresponds to the case of a single imputed path that is the posterior mean of the AID, and the far right corresponds to the case of a single imputation of the true, simulated path. Thus, the far right values in each subplot represent a baseline against which we can compare the performance of the other regimes. Figure 8 in Appendix A shows the coverage proportions for $\sigma_{s}^{2}$.

From Figure 4, it is clear that merely imputing the posterior mean of the AID, rather than several samples, does not adequately represent the variability of the imputation distribution, leading to bias and poor coverage of the true function $\boldsymbol{\beta}$. However, using a small number of imputed paths is generally sufficient to achieve a large coverage region.

Second, in all cases, imputing from the OU-based AID yields larger coverage regions as compared to the GP-based AID. A large coverage region is not, by itself, evidence of strong performance. However, the OU-based AID also generally yields larger detection regions, suggesting that the OU-based AID better represents the true imputation distribution than the GP-based AID, at least for the particular case of the SDE-based movement model defined by (5) – (11).

Third, the sizes of the coverage and detection regions for cases corresponding to an OU AID model are generally insensitive to measurement error and observation frequency (Figure 4), with the only exception being the case of high measurement error and sparse observations. The same cannot be said of the GP AID model, for which coverage and detection of $\boldsymbol{\beta}$ deteriorates when measurement errors are small, especially when observations are dense in time. The weakness of the GP AID model for the case of small measurement error is apparent in Figure 5 where coverage drops to zero.

Finally, inference for the function $\boldsymbol{\beta}$ using the process imputation approach appears to be more robust across different values of measurement error, observation frequency, and AID models than inference for $\sigma_{v}^{2}$. Neither AID model resulted in adequate coverage for $\sigma_{v}^{2}$, with the possible exception of the OU model for the cases corresponding to small measurement error. Therefore, the results of the study suggest that process imputation is a reasonable approach for obtaining inference about the effect of the potential surface, $H$, which depends directly on the position of the individual. However, there is also evidence that process imputation may not be a valid approach for obtaining inference about the parameter $\sigma_{v}^{2}$, which relates to the velocity of the individual.

3.2 Study 2: Comparing exact and two-stage approaches

The simulation study discussed in 3.1 is useful because it demonstrates the viability of the two-stage approach to approximate inference for the case when exact inference is computationally infeasible. By the same token, we cannot compare the approximate inference from the two-stage procedure to exact inference. Any model for which we could obtain both forms of inference must necessarily be somewhat contrived, because the very accessibility of exact inference negates the need for an approximate procedure. Still, it is instructive to consider a simplified setting when both exact and approximate inference is available to see what differences exist.

We investigate the model employed by Brillinger and Stewart (1998) for the continuous movement of a particle under the influence of an external force. As earlier, let $\boldsymbol{\mu}(t)$ denote the location of a particle in two dimensions at time $t$. Let $\mathbf{H}(\boldsymbol{\mu}(t),\beta)$ be a potential surface whose negative gradient represents the force operating on the particle, the value of which may depend upon the parameter $\beta$. In contrast to the model used in the first simulation study, we assume the parameter $\beta$ is constant in time. The continuous path of the particle is modeled with the stochastic differential equation

[TABLE]

where $d\mathbf{b}(t)$ are infinitesimal increments of Brownian motion. As before, we model the measured locations, $\mathbf{s}(t)$, of the particle as conditionally Gaussian with mean equal to the true location, $\boldsymbol{\mu}(t)$, and variance $\sigma_{s}^{2}$.

An important distinction between this model and that presented in Section 3 is that the “force” defined by the potential function now operates directly on the position process, rather than the velocity. Therefore, interpretation of the effect in the context of classical Newtonian physics is no longer possible, and continuing to refer to the function $\mathbf{H}$ as a potential function and its gradient as a force may be misleading. Additionally, paths generated from this model do not have well-defined velocities at all time points, and may not be realistic models for the movement of massive bodies.

Nevertheless, the stochastic differential equation is sometimes sufficient for approximating the movement of particles under the influence of an external force. We consider the case when $\mathbf{H}(\boldsymbol{\mu}(t),\beta)=\beta\|\boldsymbol{\mu}(t)-\mathbf{c}\|_{2}$, as in Section 2.2, such that the particle feels a constant “push” toward the location $\mathbf{c}$.

Similar to the simulation study in Section 3.1, we fit the model to simulated data using the proposed two-stage approximation. For the purposes of this simulation study, we used the Ornstein-Uhlenbeck model of Johnson et al. (2008) for the AID and obtained inference using the posterior mean of the fitted AID, as well as $K=8,32,$ and $128$ draws from the AID. Additionally, we obtained exact inference for the model in an MCMC framework by incrementally updating the true, unobserved path $\boldsymbol{\mu}$ over a dense grid of time points $t_{1},\dots,t_{m}$ (see Appendix A for further details). Settings with both “sparse” and “dense” rates of observation (100, 500 observation times, respectively), as well as “large” and “small” measurement errors ($\sigma_{s}^{2}=10^{-2},10^{-4}$ respectively) were considered (see Figure 9 in Appendix A for examples).

We evaluated the performance of each method by examining coverage and detection rates (see Section 3) across 24 different simulated paths. The left half of Figure 6 shows proportion of simulated paths for which the equal-tailed 95% posterior credible interval overlapped the true value (“o”), and the proportion of simulated paths for which the credible intervals both overlapped the true value, and failed to contain 0 (“+”). Coverage was similar for both methods, across all choices of $K$, however detection was low for the process imputation approach when only the posterior mean of the AID was used. When 8 or more draws from the AID were used, detection rates for the two-stage process imputation method was similar to those obtained from exact inference, although for the case of “sparse” observation times, detection rates for the process imputation approach were slightly lower than those obtained for exact inference, suggesting that the two-stage procedure may result in inflated uncertainty about $\beta$ when the density of observations is low.

The right half of Figure 6 shows the coverage rates for the measurement error variance, $\sigma_{s}^{2}$. All methods performed similarly. The findings of the simulation study suggest that the two-stage approach is a viable option when multiple draws from the AID are used.

4 Application

We analyzed telemetry data gathered for a single female northern fur seal during the fall of 2008 (Figure 1). Observations were collected using high precision GPS devices at intervals of approximately 15 minutes (for more details see Johnson et al. (2013)) over the course of 6 weeks. The individual under study makes several long foraging excursions, returning approximately every 7 days to nurse pups at a rookery on Saint Paul Island. A vignette giving a detailed description of the model fitting procedure is provided in Appendix C.

We used process imputation to fit the SDE-based model given by (5) – (11), defining the center of attraction, $\mathbf{c}$, to be the location of the rookery, $(57.13^{\circ}\text{N},170.27^{\circ}\text{W})$. Informed by the results of the simulation study (Section 3), we selected the Ornstein-Uhlenbeck process of Johnson et al. (2008) as the AID model, and used $K=128$ imputed paths.

Figure 7 shows a pointwise equal-tailed 95% credible interval for the dynamic effect of the rookery on the movement of the fur seal. The dashed line shows the distance the animal was from the rookery over time. The periodic pattern of repulsion ($\beta(t)<0$) and attraction ($\beta(t)>0$) that the fur seal experiences with respect to the rookery coincide with foraging excursions. A subtler pattern of behavior is visible in the fine scale oscillations that occur during the individual’s travel away from the rookery (e.g., from approximately September 7 – 10). During this phase of the foraging excursions, the force of repulsion away from the rookery is varying, indicating movement that is less directed and more exploratory. In contrast, during the second half of each foraging excursion (e.g., approximately September 12 – 16), the force of attraction steadily increases until just before arrival, indicating strongly directed movement.

Several other relevant questions of scientific interest related to the patterns observed in $\boldsymbol{\beta}$ can be addressed using approximate inference from the SDE-based model. Examples of how this model could be employed include: (1) the study of patterns in $\boldsymbol{\beta}$ over the course of a long time period to learn about the physical development of northern fur seal pups as demonstrated by their mother’s changing behavior; (2) investigation of models with additional potential surface terms for other factors from which researchers could learn about the effects of important environmental covariates on movement, while controlling for the effect of the rookery; and (3) comparison of the function $\boldsymbol{\beta}$ across multiple animals to learn about individual-level variability among nursing females. Without using process imputation, mechanistically-motivated models such as this one are inaccessible to researchers studying animal movement.

5 Discussion

A key result from our simulation study for the case for SDE-based movement models is that process imputation is an effective approach for studying some types of model parameters, but not necessarily others. The approximate posterior distribution for the function $\boldsymbol{\beta}$ is reasonably robust to the severity of measurement error, observation frequency, and AID model, as evidenced by the close agreement in coverage and detection between inference using process imputation and inference conditioned on the true movement process $\boldsymbol{\mu}$ (Figure 4). In contrast, approximate inference for the parameter $\sigma_{v}^{2}$ varies substantially with the choice of AID model, and appears to be strongly dependent on the severity of measurement error and observation frequency (Figure 5).

A practical interest for researchers is understanding the types of variables for which approximate inference using process imputation will be reliable. We suggest one possible indicator for reliability based on a distinction between “position-order” and “higher-order” effects on movement.

Generally, we define position-order effects on movement to be those directly related to the position process of an individual. We define higher-order effects to be those related to derivative processes such as velocity or acceleration. For the SDE-based movement model, the magnitude and direction of the force induced by the potential surface is defined by the the particle’s location in space and is therefore primarily a position-order effect. Alternatively, the influence of the parameter $\sigma_{v}^{2}$ is through the particle’s velocity rather than its location, and is therefore a higher-order effect. Telemetry data typically contain large amounts of information about an individual’s position. However, information in the data about a particle’s higher-order processes can be far less precise and deteriorates much more rapidly with increasing measurement error and gaps between observations. Therefore, we can expect parameters such as $\sigma_{v}^{2}$ to be more difficult to estimate and more strongly influenced by the choice of AID model in a process imputation framework. Currently, process imputation-based inference for higher-order parameters may be unreliable. However, many modern telemetry devices now collect direct measurements of an individual’s velocity and acceleration, which may allow for a broader use of process imputation in the future. Additionally, the primary parameter of scientific interest in the SDE-based movement model is $\boldsymbol{\beta}$; bias in the estimation of higher-order parameters such as $\sigma_{v}^{2}$ may be ignorable when it does not affect ecological learning.

In selecting an AID, the goal is to represent the distribution of the true path conditioned on the data as faithfully as possible. Thus, a first step is to ensure that candidate AIDs generate paths that are physically proximate to those arising from the true imputation distribution. Researchers should verify that the paths from candidate AIDs lie near the observed data. A greater challenge is to verify that the form of temporal dependence present in the AID is similar to that of the true imputation distribution, because the structure of dependence can vary considerably without affecting the central measure. For example, the OU and GP AIDs implemented in the simulation study in Section 3 generate paths with very different levels of smoothness, in the sense that the former has no continuous derivatives, while the latter has an infinite number. As a general rule, one should select a model for the AID whose dependence structure is similar to that of the process model. In practice, the family of models available for the AID will be limited by various constraints, including computational ones, and the researcher will be required to balance several competing priorities. For many applications, it may be helpful to consider multiple potential AIDs and investigate the sensitivity of all parameters of interest.

A common family of movement models that make use of process imputation are discrete-space models (Hooten et al., 2010; Hanks et al., 2015; Hooten et al., 2016; Hanks and Hughes, 2016), which model an individual’s movement as a series of transitions between areal regions. Each transition is modeled using a multinomial distribution, where the probability of an individual moving to an adjacent areal region is a function of the associated environmental covariates. The effects of the environmental covariates on these transition probabilities are examples of position-order effects because they depend explicitly on the location of the individual (i.e., the areal region currently occupied). While the transition probabilities also depend on the residence time of an individual within each grid cell, and thus the average velocity of the animal during its residence, they do not depend explicitly on the instantaneous velocity at each time point. Our simulation study did not directly investigate the performance of process imputation for discrete-space movement models, but our findings suggest that process imputation may provide reliable inference about the relationship between environmental covariates and the inter-cell transition probabilities.

Although imputation-based approaches are an increasingly common tool for fitting models of animal movement to data, we are not aware of any other studies of their impact on inference. Our findings provide useful insight about the circumstances for which process imputation delivers valid inference on parameters of interest.

Acknowledgments

The authors thank Ephraim Hanks for early insights and discussions about the research. Funding for this research was provided by NOAA (RWO 103), CPW (TO 1304), and NSF (DMS 1614392). Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Appendix

Appendix A Simulation studies

A.1 Coverage and detection regions

To evaluate the performance of the process imputation approach, we compared the posterior distributions of $\boldsymbol{\beta}$ to the values used to simulate the data through two summary measures. First, we computed the proportion of the time interval during which a pointwise, equal-tailed 95% credible interval contained the true value of the function $\boldsymbol{\beta}$, which we term the “coverage region.” Small coverage regions suggest that the associated method is biased or under-represents uncertainty about model parameters and may be unreliable in practice. In contrast, large coverage regions do not necessarily suggest reliability, because uncertainty about model parameters may be inflated, resulting in high coverage rates without yielding precise inference. A desirable property of posterior credible regions for model parameters is that they have high coverage, while remaining narrow enough to allow researchers to say something meaningful about relevant parameters when an effect on movement is present. We therefore computed a second summary measure defined as the proportion of the time interval during which the credible interval did not contain zero, which we term the “detection region,” because it represents regions for which there is substantial evidence that $\beta(t)\neq 0$. Figure 3 in the main body of the manuscript shows one example of coverage and detection regions. Our definitions of “coverage” and “detection” regions are not intended to be interpreted from a hypothesis testing perspective. Rather, they are intended to provide insight into the amount of bias and uncertainty that occur under the process imputation inferential procedure.

A.2 Study 1: Brownian particle in an external force field

A.2.1 MCMC details

The Gaussian process AID model is specified as $\boldsymbol{\mu}^{*}\sim\text{N}\left(\mathbf{0},\mathbf{C}^{*}\right)$, where the covariance of $\boldsymbol{\mu}^{*}(t_{j})$ and $\boldsymbol{\mu}^{*}(t_{j})$ is defined as $\mathbf{C}^{*}(t_{j},t_{k})=e^{-\frac{(t_{j}-t_{k})^{2}}{2\phi^{2}}}$. The Ornstein-Uhlenbeck AID model is specified through the velocity, and was first introduced in Johnson et al. (2008). We used the R package crawl (Johnson, 2016). The Ornstein-Uhlenbeck model has been used in several previous instances of process imputation (e.g., Hanks et al., 2011, 2015; Scharf et al., 2016).

A.2.2 Results

Inference for the variance in the measurement error process was consistent across almost all simulation scenarios. The exception was for the case of small measurement errors and dense observation times when using the Gaussian process AID (bottom right, Figure 8). Under this combination of circumstances, the GP AID proved too inflexible to accurately represent the true imputation distribution.

A.2.3 Timing

On a personal laptop with 8 GB of random access memory, and a 2.6 GHz Intel Core i5 processor, the combined two-stage fitting procedure required approximately 10 minutes of computation time per simulated data set.

A.3 Study 2: Comparing exact and two-stage approaches

A.3.1 Exact inference

We estimated the unobserved, true location of the particle over a grid of 500 time points in addition to the observation times. The conditional distribution of $\boldsymbol{\mu}(t_{i})$ is given by

[TABLE]

where $\Delta t_{i}\equiv t_{i}-t_{i-1}$. The conditional distributions for the observed data are given by

[TABLE]

Priors/Hyperparameters

[TABLE]

Full Conditionals

True path:

[TABLE]

Strength of attraction:

[TABLE]

Measurement error variance:

[TABLE]

A.3.2 Timing

On a personal laptop with 8 GB of random access memory, and a 2.6 GHz Intel Core i5 processor, the combined two-stage fitting procedure required approximately 12 seconds of computation time per simulated data set, while the exact-inferential procedure required approximately 10 minutes of computation time. Caution should be exercised in generalizing the relative computational costs of the two inferential procedures explored here to other process imputation applications.

A.3.3 Sample paths

Appendix B Application

Inference is based on a chain of $10^{4}$ iterations, with the first half discarded as burn-in. Evidence of convergence was determined by visual inspection of trace plots, and the Gelman-Rubin diagnostic (Gelman et al., 2014). We fit the model using a grid of values ($10^{i/4}$ for $i\in\left\{-16,\dots,12\right\}$) for the hyper-parameter $\sigma^{2}_{\alpha}$, and found that $10^{-1.75}$ provided an optimal amount of tuning as measured by DIC. As an alternative to score-based selection of the tuning parameter, one could instead specify a prior distribution on $\sigma^{2}_{\alpha}$. Our approach has the benefit of regularizing $\beta(t)$, which can provide improved estimation (e.g., Hastie et al., 2009).

B.1 Timing

On a personal laptop with 8 GB of random access memory, and a 2.6 GHz Intel Core i5 processor, the combined two-stage fitting procedure required approximately 80 minutes of computation time.

Appendix C Application vignette

This vignette shows how the two-stage process imputation procedure was implemented for the application to the movement of a Norther fur seal.

http://www.stat.colostate.edu/~scharfh/supplemental_multiple_imputation_JABES/