Process convolution approaches for modeling interacting trajectories

TL;DR

This paper introduces a flexible multi-stage process convolution framework for modeling complex interacting trajectories, demonstrated through social whale movement, improving inference and reducing uncertainty in dependent path modeling.

Contribution

It generalizes process convolutions into a chain structure, enabling interpretable modeling of complex dependencies in spatio-temporal data.

Findings

Model captures social interactions in whale movement.

Reduces uncertainty in trajectory estimation.

Provides interpretable insights into dynamic social structures.

Abstract

Gaussian processes are a fundamental statistical tool used in a wide range of applications. In the spatio-temporal setting, several families of covariance functions exist to accommodate a wide variety of dependence structures arising in different applications. These parametric families can be restrictive and are insufficient in some situations. In contrast, process convolutions represent a flexible, interpretable approach to defining the covariance of a Gaussian process and have modest requirements to ensure validity. We introduce a generalization of the process convolution approach that employs multiple convolutions sequentially to form a "process convolution chain." In our proposed multi-stage framework, complex dependencies that arise from a combination of different interacting mechanisms are decomposed into a series of interpretable kernel smoothers. We demonstrate an application of process convolution chains to model killer whale movement, in which the paths taken by multiple individuals are not independent, but reflect dynamic social interactions within the population. Our proposed model for dependent movement provides inference for the latent dynamic social structure in the study population. Additionally, by leveraging the positive dependence among individual paths, we achieve a reduction in uncertainty for the estimated locations of the whales, compared to a model that treats paths as independent.

Tables

Table 1. Table 1: True values and prior distibutions used in simulation study.

parameter	true	prior density
${\varphi }_{inl}$ (“low tortuosity")	0.04	$\text{Gamma}(2,100)$
${\varphi }_{inl}$ (“high tortuosity")	0.04/3	$\text{Gamma}(2,100)$
${\sigma }_0^2$	1	$\text{IG}({10}^{-3},{10}^{-3})$
${\sigma }_{\mu /s}^2$	800	$\text{IG}({10}^{-3},{10}^{-3})$
${\sigma }_{\mu }^2$	10
${\sigma }_s^2$	0.0125	$\text{IG}({10}^{-3},{10}^{-3})$

Table 2. Table 2: Posterior credible intervals and prior distributions for the IP-DEP model in killer whale application.

	posterior (IP-DEP)		posterior (IP-IND)		prior
parameter	median	(2.5%, 97.5%)	median	(2.5%, 97.5%)	density
${\varphi }_{inl}$	0.00848	(0.00756, 0.00954)	0.00850	(0.00762, 0.00951)	$\text{Gamma}(2,100)$
${\sigma }_0^2$	705	(292, 2040)	1440	(627, 4510)	$\text{IG}(1,10)$
${\sigma }_{\mu /s}^2$	3710	(2600, 5490)	1700	(1380, 2400)	$\text{IG}({10}^{-3},{10}^{-3})$
${\sigma }_{\mu }^2$	4.94	(3.52, 7.21)	2.39	(1.88, 3.14)	NA
${\sigma }_s^2$	0.00133	(0.00122, 0.00145)	0.00133	(0.00123, 0.00145)	$\text{IG}({10}^{-3},{10}^{-3})$
${\sigma }_w^2$	0.491	(0.370, 0.660)	NA	NA	$\text{IG}(52,10)$

Equations

$$\mu (t)$$

$$\text{Cov}\left(\mu (t),\mu (t^{\ast })\right)$$

$$\mu (t)$$

$$\begin{array}{rl}{\mu }^{(1)}({\tau }_1)& ={\int }_{\mathcal{T}}{h}^{(1)}({\tau }_1,{\tau }_0)dB({\tau }_0)\\ {\mu }^{(\text{uline}2)}({\tau }_2)& ={\int }_{\mathcal{T}}{h}^{(2)}({\tau }_2,{\tau }_1)d{\mu }^{(1)}({\tau }_1)\\ ⋮& \\ {\mu }^{(\text{uline}L)}({\tau }_l)& ={\int }_{\mathcal{T}}{h}^{(L)}({\tau }_L,{\tau }_{L-1})d{\mu }^{(\text{uline}L-1)}({\tau }_{L-1}).\end{array}$$

$${\mu }^{(\text{uline}l)}({\tau }_l)$$

$$h^{(\text{uline}l)}({\tau }_l,{\tau }_0)$$

$${\int }_{\mathcal{T}}{\int }_{\mathcal{T}}{h}^{(l_1)}(t,\tau )h^{(l_2)}(\tau ,{\tau }_0)dB({\tau }_0)\ne {\int }_{\mathcal{T}}{\int }_{\mathcal{T}}{h}^{(l_2)}(t,\tau )h^{(l_1)}(\tau ,{\tau }_0)dB({\tau }_0).$$

$$\tilde{\mathbf{\mu }}\sim \mathcal{N}\left(0,\Delta \tau \tilde{\mathbf{H}}{\tilde{\mathbf{H}}}^{\prime }\right).$$

$$h_{ij}^{(soc)}({\tau }_{soc},\tau )$$

$$|w_{i\cdot }({\tau }_{soc})|$$

$${\mu }_1^{(soc)}({\tau }_{soc})$$

$$=h_{11}({\tau }_{soc}){\mu }_1({\tau }_{soc})+h_{12}({\tau }_{soc}){\mu }_2({\tau }_{soc})$$

$$=\frac{{\mu }_1({\tau }_{soc})}{1+|w_{12}({\tau }_{soc})|}+\frac{w_{12}({\tau }_{soc}){\mu }_2({\tau }_{soc})}{1+|w_{12}({\tau }_{soc})|}.$$

$${\mu }_i^{(bm)}({\tau }_{bm})$$

$${\mu }_i^{(\text{uline}soc)}({\tau }_{soc})$$

$${\tilde{\mu }}_i({\tau }_{inl})={\mu }_i^{(\text{uline}inl)}({\tau }_{inl})$$

$$s_i(t)$$

$$w_{ij}({\tau }_{soc})$$

$${\mathbf{\mu }}_i^{(w)}({\tau }_w)$$

$$h^{(w)}({\tau }_w,\tau )$$

$${\mathbf{\mu }}_i^{(w)}\sim \mathcal{N}\left(0,{\sigma }_w^2\Delta t{\tilde{\mathbf{H}}}_w^{\prime }{\tilde{\mathbf{H}}}_w\right).$$

$${\mathbf{s}}_i(t)|{\mathbf{\epsilon }}_w,{\varphi }_{inl},{\sigma }_s^2,{\sigma }_{\mu }^2,{\sigma }_0^2\sim \mathcal{N}\left(0,\mathbf{\Sigma }\right)$$

$$\mathbf{\Sigma }\equiv {\sigma }_s^2\mathbf{I}+{\sigma }_{\mu }^2\Delta \tau {\tilde{\mathbf{H}}}^{\prime }\tilde{\mathbf{H}},$$

$$\text{SPE}(\hat{\mathbf{\mu }};{\mathcal{T}}_g,{\tilde{\mathbf{\mu }}}_{true})=\frac{1}{|{\mathcal{T}}_g|}{\int }_{{\mathcal{T}}_g}\Vert {\tilde{\mathbf{\mu }}}_{true}(t_g)-\hat{\mathbf{\mu }}(t_g){\Vert }_2^{2}d{t}_g.$$

$$\mathbf{s}$$

$$\tilde{\mathbf{H}}$$

$${\mathbf{\Sigma }}_{dB}(t,\tau ,i,j)$$

$${\mathbf{H}}^{(bm)}(t,\tau ,i,j)$$

$${\mathbf{H}}^{(soc)}(t,\tau ,i,j)$$

$${\mathbf{H}}^{(inl)}(t,{\tau }_{soc},i,j)$$

$$\mathbf{\tau }$$

$$\Delta {\tau }_i$$

$$dB({\tau }_i)$$

$$d\mathbf{B}$$

$$\tilde{\mathbf{\mu }}$$

$${\varphi }_{inl}$$

$${\sigma }_{\mu /s}$$

$${\sigma }_w^2$$

Full text

\jnlcitation\cname

, , , and , (\cyear2017), \ctitleProcess convolution approaches for modeling interacting trajectories, \cjournalEnvironmetrics.

1 Introduction

Models for continuous random processes using kernel convolutions are known in much of the spatial and spatio-temporal statistical literature as “process convolutions” (e.g., Higdon 23; Calder 6; Bolin \BBA Lindgren 3). Also referred to as spatial moving averages e.g., 10, process convolutions first arose as a means for constructing valid covariance matrices for Gaussian processes (GPs), while relaxing the typical assumptions of stationarity and isotropy e.g., 2. One application area for process convolutions that has driven much of their recent theoretical development is for random spatial processes on stream networks 42. Stream networks occupy sub-manifolds of physical space on which specialized measures of distance often complicate enforcement of non-negative definiteness in covariance functions and have an inherent need for anisotropic dependence structures to accommodate the effects of direction in stream flow. Process convolutions offer an approach to modeling dependence that addresses both of these issues. We employ process convolutions to model the dependence observed in the trajectories of interacting animals.

A wealth of data arising from measurements of trajectories has spurred the development of many new statistical models for individual-based movement processes. These models have been used to study the individual movement patterns of a wide variety of animals, where scientific questions of interest are often focused on which exogenous environmental factors drive movement 35. A fundamental commonality of research related to animal movement is a need for accurate and precise knowledge about the paths traversed by each individual, leading to the development of many models used for reconstructing the true underlying movement processes from telemetry data e.g., 16. Several excellent models for telemetry data have been proposed in both discrete- and continuous-time formulations that have sought to address the significant challenges making prediction of the true paths difficult (for a review see McClintock \BOthers. 31; Hooten \BBA Johnson 26). Telemetry devices are often subject to environmental conditions, and technological limitations that restrict their ability to deliver precise measurements of location at dense temporal resolutions. Thus, many analyses of animal movement must account for both measurement error and irregular temporal observations when estimating the true path traversed by an individual. In some cases, it may be reasonable to assume that the data include independent Gaussian noise; however, for many telemetry devices, it is necessary to account for complex, non-Gaussian forms of measurement error 4, 5. Irregularity in the frequency of telemetry observations can complicate the implementation of discrete-time models although see 39 for one solution using a multiple-imputation procedure and potentially lead to large uncertainty during time intervals with few or no observations.

Ecologists are also interested in the extent to which individuals in a population move in direct response to each other. Particularly in the case of large mammals, complex and dynamic social connections within a population can play a critical role in the movement behavior of individuals (e.g., Williams \BBA Lusseau 46; S\BPBIZ. Goldenberg \BOthers. 19; Scharf \BOthers. 38). For example, two individuals with a strong social connection may exhibit similar movement patterns and visit the same locations, resulting in positive dependence between their respective paths. Inference about the social ties within a population, and how those ties change over time, can provide valuable information about the behavioral ecology of a population. However, there are few existing methods available that explicitly account for interactions among individuals (e.g., Haydon \BOthers. 21; Codling \BBA Bode 8; Langrock, Roland and Hopcraft, J. Grant C. and Blackwell, Paul G. and Goodall, Victoria and King, Ruth and Niu, Mu and Patterson, Toby A. and Pedersen, Martin W. and Skarin, Anna and Schick, Robert S. 30; Russell \BOthers. 37; Scharf \BOthers. 38). Especially when direct observation of a species is infeasible, careful analysis of data gathered from telemetry devices can reveal useful information about a population’s social structure 38.

We propose a joint model for the movement of multiple individuals that expands the suite of statistical methodology available for studying animal movement. Our approach allows researchers to study dependence among the paths of a population of individuals that arises as the result of social interactions. Also, by taking into account the dependence among individuals in a population, our approach has the potential to reduce both bias and uncertainty in reconstructing trajectories in the study population compared to models that treat individuals independently. We obtain improvements in path reconstruction by modeling the true movement processes of the population conditioned on an unobserved, dynamic social network. The inferred network provides a description of the social ties within a population and how those ties change in time.

In what follows, we present a generalized approach to process convolutions in which a random process is decomposed into a sequence of one or more smoothing kernels that are convolved with a white noise process. Multiple stages of smoothing allow us to model dependence in time and between individuals separately, resulting in a multivariate Gaussian process that captures the combination of effects. It is important to note that, while the kernels responsible for inducing temporal and path-wise dependence are constructed separately, the resulting covariance function is not separable in the geostatistical sense.

We provide a detailed introduction to our generalized approach to process convolutions in Sections 2.1 and 2.2 and construct our joint model for movement in Sections 2.3 and 2.4. After discussing the implementation details in Section 3, we demonstrate path reconstruction in a simulation study in Section 4. We apply our model to the study of telemetry data arising from the movement of killer whales in Section 5 and close with a summary of our findings and future directions.

2 Methods

2.1 A multiple-kernel convolution framework

We first outline a new flexible framework for the development of Gaussian process models based on a synthesis of ideas from geostatistics, multivariate time series, and trajectory modeling. Our hierarchical framework relies on the kernel convolution approach for modeling random processes, known, in the spatial and spatio-temporal statistical literature, as “process convolutions” (e.g., Higdon 23; Calder 6; Bolin \BBA Lindgren 3). A mean zero random process $\mu(\cdot)$ is called a process convolution if it is constructed by convolving a continuous random process, $d\!B(\cdot)$, with a kernel function, $h$, over a domain, $\mathcal{T}$, so that

[TABLE]

If the random process $d\!B(\cdot)$ is Gaussian, the resulting process $\mu(\cdot)$ will also be Gaussian with covariance function

[TABLE]

It is common to define the process $d\!B(\cdot)$ to be Gaussian white noise, which yields the simplified covariance function

[TABLE]

(although see Nychka \BOthers. 33 for an example of a spatial process convolution where $d\!B(\cdot)$ has covariance specified through a Gaussian Markov random field). A process convolution with kernel function $h(\cdot,\cdot)$, as in (1), represents a smoothing of the process $d\!B(\cdot)$, and the kernel is therefore often referred to as a smoother. In what follows, we refer to $h$ as both a kernel function and a smoother, interchangeably.

The function defined in (3) is guaranteed to be non-negative definite, and therefore a valid covariance function, if $\int_{\mathcal{T}}h(t,\tau)d\tau<\infty$ and $\int_{\mathcal{T}}h^{2}(t,\tau)d\tau<\infty$ for all $t$ 23. In many applications, it may be easier to specify the proper form of dependence in a GP through the form of a kernel smoother rather than that of a covariance function. For example, by constructing process convolutions with asymmetrical kernels, Ver Hoef \BOthers. 42 accounted for directional dependence in spatial processes arising on stream networks, where directional flow plays a crucial role. For certain choices of $h$, it is possible to express the covariance function analytically (e.g., Higdon 23; Paciorek 34; Ver Hoef \BBA Peterson 43), although this is not necessary for non-negative definiteness. In the cases where an analytic solution to (3) is not available, one can evaluate the integral numerically. The class of covariance functions constructed using kernel convolutions is general, containing many of the parametric families commonly used in geostatistical settings (e.g., exponential, Gaussian, spherical), as well as more flexible dependence structures. For example, Higdon 22 and Paciorek 34 showed that process convolutions allow researchers to model both anisotropic and non-stationary dependence by letting $h(t,\tau)$ vary with both $|t-\tau|$ and $t$. In principle, any covariance function evaluated over a discretized domain has a process convolution representation on the same discretized grid, where the kernel function’s values over the grid may be defined using a decomposition of the covariance matrix (e.g., Cholesky).

Hooten \BBA Johnson 26 used process convolutions to develop new models for trajectories by convolving a Wiener process, rather than white noise, with different kernel functions. Using Wiener process convolutions, they specified realistic models for animal movement without the characteristic “roughness” found in Brownian motion and define a framework that incorporated other existing models for movement, such as that in 27. Brownian motion itself can be thought of as a process convolution where the kernel is a step function given by $h^{(bm)}(t,\tau)\equiv\mathbf{1}_{\left\{\tau\leq t\right\}}$ and the process is Gaussian white noise. Thus, the Hooten \BBA Johnson 26 framework can be viewed as a nested smoothing procedure, where the position at time $t$ is

[TABLE]

To account for Brownian motion with an initial position far from the origin, one may define the process $d\!B(\cdot)$ to be a continuous white noise process for all $t>0$, with arbitrary initial variance so that $d\!B(0)\sim\mathcal{N}\left(0,\sigma_{0}^{2}\right)$. Hooten \BBA Johnson 26 provided several examples of possible kernel functions and demonstrated a computationally efficient procedure to incorporate temporal non-stationarity into the kernel, allowing them to fit highly flexible models to telemetry data arising from animal movement trajectories.

It is possible to write the two-stage process convolution in (4) in terms of a single effective smoothing kernel $\tilde{h}(t,\tilde{\tau})=\int_{\mathcal{T}}h(t,\tau)h^{(bm)}(\tau,\tilde{\tau})d\tau$ using Fubini’s theorem to change the order of integration. Moreover, the equivalent, single-stage representation for $\mu(\cdot)$ is not limited to two-stage smoothing processes. One may specify any number of kernel functions and the resulting multi-stage kernel convolution can be written in terms of a single effective kernel. This nested structure motivates a novel model-building approach in which an arbitrarily large collection of convolution kernels is used to specify a valid GP. We refer to a stochastic process resulting from the iterative convolution of a chain of kernels $h^{(1)},\dots,h^{(L)}$ as a “process convolution chain” (PCC) and write it in its expanded form as

[TABLE]

The underline in stages $2\leq l\leq L$ (i.e., $\mu^{(\uline{l})}(\tau_{l})$) distinguishes the PCC composed of all kernels $\left\{h^{(k)}:k\leq l\right\}$ from that of convolving $h^{(l)}$ with the process $d\!B(\cdot)$ directly (written $\mu^{(l)}(\tau_{l})$). We use matching superscripts to denote the covariance function for a given process (e.g., $\text{Cov}(\mu^{(\uline{l})}(t),\mu^{(\uline{l})}(t^{*}))=C^{(\uline{l})}(t,t^{*})$). A collapsed version of a PCC using a single effective kernel and the white noise process $d\!B(\cdot)$ can be written as

[TABLE]

To simplify notation, we write $\tilde{h}=h^{(\uline{L})}$ for the effective kernel constructed from the entire chain, $\tilde{\mu}$ for the full process $\mu^{(\uline{L})}$, and $\widetilde{C}$ for the corresponding covariance function. In general, the kernels $h^{(k)}$ are not commutative, in the sense that

[TABLE]

Thus, both the forms of the kernels and their order are important when specifying a PCC.

Specifying a GP through an ordered chain of smoothers allows for considerable flexibility. For example, if we specify a chain of size $L=1$ with kernel function $h^{(1)}(t,\tau_{1})\equiv\mathbf{1}_{\left\{\tau_{1}<t\right\}}$, we recover simple Brownian motion. If we increase the length of the chain by including a second kernel, $h^{(2)}(t,\tau_{2})=(1+\frac{\tau_{2}-t}{\phi})\mathbf{1}_{\left\{-\phi<\tau_{2}-t\leq 0\right\}}$, we recover a structure used by Ver Hoef \BBA Peterson 43 for modeling dependence in stream networks also mentioned in 26 Figure 2, row 3. Another specific case of PCCs are so-called $(k-1)$-fold integrated Wiener processes 40, 44, 36, which use $L=k-1$ kernels of the form $h^{(l)}(t,\tau)=(t-\tau)/(l-1)$ and allow one to model GPs with exactly $m$ continuous derivatives. The multi-stage decomposition is similar in spirit to the way random processes with multiple dependence scales are decomposed additively in multi-resolution processes (Higdon 23; Nychka \BOthers. 33; Katzfuss 29); the difference is the convolution chain approach combines model components using convolution rather than addition.

In what follows, we construct convolution kernels that allow us to specify multivariate GPs; these provide a powerful tool for modeling multiple trajectories arising from interacting individuals. We return to this specific application of PCCs in Section 2.3.

2.2 Finite representation

The joint distribution of $\tilde{\boldsymbol{\mu}}\equiv\left(\tilde{\mu}(t_{1}),\dots,\tilde{\mu}(t_{n})\right)^{\prime}$ for $\mathbf{t}\equiv\left(t_{1},\dots,t_{n}\right)^{\prime}\subseteq\mathcal{T}$ is mean-zero Gaussian with covariance function $\widetilde{C}(t,t^{*})=\int_{\mathcal{T}}\tilde{h}(t,\tau)\tilde{h}(t^{*},\tau)d\tau$ for any $t,t^{*}\in\mathbf{t}$ (if the integrals of $\tilde{h}$ and $\tilde{h}^{2}$ exist). When this integral cannot be computed analytically, one can perform the integration numerically by selecting a fine grid of values $\boldsymbol{\tau}\equiv\left(\tau_{0},\dots,\tau_{m}\right)^{\prime}$ from the domain $\mathcal{T}$, and evaluating the discrete sum $\sum_{i=1}^{m}\tilde{h}(t,\tau_{i})\tilde{h}(t^{*},\tau_{i})\Delta\tau_{i}$, where $\Delta\tau_{i}=\tau_{i}-\tau_{i-1}$. In the discretized setting, the process convolution may be written as a matrix product $\tilde{\boldsymbol{\mu}}=\widetilde{\mathbf{H}}\boldsymbol{\varepsilon}_{0}$, where the $j^{th}$ row of $\widetilde{\mathbf{H}}_{n\times m}$ is the function $h(t_{j},\tau)$ evaluated at all $\tau\in\boldsymbol{\tau}$, and $\boldsymbol{\varepsilon}_{0}\equiv\left(\varepsilon_{0}(\tau_{1}),\dots,\varepsilon_{0}(\tau_{m})\right)^{\prime}$ with each $\varepsilon_{0}(\tau_{i})\sim\mathcal{N}(0,\Delta\tau_{i})$. Assuming the grid times $\boldsymbol{\tau}$ are equally spaced with intervals of size $\Delta\tau$, the joint distribution of $\tilde{\boldsymbol{\mu}}$ can be expressed using matrix notation as

[TABLE]

As the density of the grid grows to infinity ($m\rightarrow\infty$), the outer product in (9) approaches the covariance defined by the continuous integral over $\mathcal{T}$. The granularity of the grid (i.e., the size of $m$) required to adequately approximate the integral will generally depend on the characteristics of $\tilde{h}$. Alternatively, one can follow the approach of Higdon 23, who used finite process convolutions also called discrete process convolutions in 7 as an approximation, and choose $m<n$ to yield a fixed-rank model for the continuous process. Fixed-rank models can offer considerable computational efficiency and provide an adjustable level of implicit regularization 45.

2.3 A process convolution chain for dependent movement

2.3.1 Social smoothing

We construct a novel model for animal movement using the PCC approach within a Bayesian hierarchical modeling framework. To describe the basic procedure, we first consider paths in one dimension. The model we describe can be readily extended to movement in two or more dimensions, and we demonstrate that in the example that follows. We use a three-stage ($L=3$) PCC that includes a specialized kernel function constructed to induce dependence among the paths of $p$ different individuals. The inter-path dependence that arises is based on a weighted, undirected, latent social network 18.

Let the random variable $\mu_{i}(t)$ represent the position of individual $i$ at time $t$, and let $w_{ij}(\tau)\in\left[0,1\right]$ denote the connection weight between individuals $i$ and $j$ at time $\tau$, with $w_{ii}(\tau)\equiv 1$ and $i,j\in\left\{1,\dots,p\right\}$. We propose a “social” smoothing kernel of the form

[TABLE]

We write the kernel in (10) with explicit dependence on arbitrary time $\tau$ for completeness, but will hereafter suppress the second argument for brevity. To provide intuition for the effect of convolving with the social smoothing kernel, we first consider the case of a multivariate process with two individuals. We assume that, for the process $\boldsymbol{\mu}(\cdot)\equiv\left(\mu_{1}(\cdot),\mu_{2}(\cdot)\right)^{\prime}$, the variables $\mu_{1}(\tau)$ and $\mu_{2}(\tau^{*})$ are independent a priori for all $\tau$ and $\tau^{*}$. After social smoothing, the resulting process for individual 1 is defined by

[TABLE]

From (13) and (14), it is clear that $\mu_{1}^{(soc)}(\tau_{soc})$ is a weighted average of the independent variables $\mu_{1}(\tau_{soc})$ and $\mu_{2}(\tau_{soc})$, where the weights are equal to $h_{11}(\tau_{soc})$ and $h_{12}(\tau_{soc})$, respectively. When the network connection is strong (i.e., $w_{12}(\tau_{soc})\approx 1$), the two weights $h_{11}(\tau_{soc})$ and $h_{12}(\tau_{soc})$ will be approximately equal, and the effect of the kernel is to “squeeze” the previously independent paths toward a scaled version of their mutual mean. When the network connection is weak ($w_{12}(\tau_{soc})\approx 0$), the weights will be close to 1 and 0 respectively, and the processes will be unaffected by the kernel, retaining the a priori independence between individuals. These effects are visible in the simulated process shown in the top plot of Figure 1.

The social smoothing kernel we propose is motivated by an interest in describing a mechanistic, interpretable driver of dependent movement. Specifically, it is useful for modeling the paths of individuals who have a tendency to move toward, and alongside, other individuals with whom they share connections. This type of dependence is ubiquitous in animal movement, but is by no means the only meaningful type of interaction. As an alternative, consider the movement of two highly territorial animals with overlapping territories in which observed paths would appear to “avoid” each other. To model mutual avoidance, we require a social kernel that “spreads apart” two otherwise independent paths when a strong social connection is present. In what follows, we focus on the particular kernel defined in (16), but we emphasize that PCCs, with suitably constructed kernels, offer a flexible way to specify models for a wide range of possible behavioral mechanisms. Additionally, our choice of social smoothing kernel is most appropriate when social connections among individuals vary slowly relative to the movement processes. When relationships among individuals are allowed to vary too rapidly, it may become difficult for the model to distinguish between brief encouters that arise due to social effects, and those that arise as part of the stochasticity inherent in Brownian motion.

Equipped with the social kernel given in (16), we outline the first two stages of our full three-stage PCC model for the dependent movement of $p$ individuals. We begin with Gaussian white noise and smooth at the first stage convolving $h_{ij}^{(bm)}(t,\tau)\equiv\mathbf{1}_{\left\{\tau\leq t\right\}}\mathbf{1}_{\left\{i=j\right\}}$ with $d\!B(\cdot)$ to generate $p$ independent instances of Brownian motion, denoted $\mu_{i}^{(bm)}(\cdot)$, each with its own initial position. At the second stage, we apply the social smoothing kernel defined in (10) to the collection of all $p$ processes, $\boldsymbol{\mu}^{(bm)}(t)\equiv\left(\mu_{1}(t),\dots,\mu_{p}(t)\right)^{\prime}$. Smoothing with the social kernel returns weighted averages of the Brownian processes, where the weights are proportional to the latent social weights $w_{ij}(t)$. Thus, two individuals $i$ and $j$ for whom $w_{ij}(t)$ is close to 1 will tend to have smoothed locations $\mu^{(\uline{soc})}_{i}(t)$ and $\mu^{(\uline{soc})}_{j}(t)$ that are close together in space. A third and final stage in the PCC, introduced in the following section, ensures that the random process has the proper temporal smoothness required for modeling animal movement.

2.3.2 Inertial smoothing

Marginally, the individual paths $\boldsymbol{\mu}_{i}^{\uline{soc}}$ generated from the two-stage PCC constructed from $h^{(bm)}$ and $h^{(soc)}$ are each an instance of Brownian motion. However, Brownain motion is unsuitable for direct modeling of animal movement because the instantaneous velocity of a particle traversing a Brownian path is discontinuous, and therefore the acceleration is not well defined. Discontinuities in the first derivative of the path processes imply that individuals are capable of instantaneous changes in their velocities and is inconsistent with the physical laws governing the mechanics of massive bodies e.g., 14.

Several process convolutions have been explored that impart specific smoothness properties on a GP. For example, Shepp 40, Wecker \BBA Ansley 44, and Rue \BBA Held 36 discuss applications of $(k-1)$-fold integrated Wiener processes, which ensure the existence of $k-1$ continuous derivatives. Johnson \BOthers. 27 constructed a model for the movement of individual harbor and northern fur seals by specifying an Ornstein-Uhlenbeck process for the velocity, rather than the position. Hooten \BBA Johnson 26 modeled the true locations of an individual animal using a two-stage PCC in which the second kernel function is Gaussian, and 5 used a similar model with cubic splines as their kernels. In all of these examples, smoothness is imparted to a process through convolution with a kernel, either implicitly or explicitly. We introduce a final stage of smoothing to yield paths that are guaranteed to be continuously differentiable.

Enforcing the existence of a continuous derivative for the movement processes ensures that individuals have an interpretable acceleration at every time $t$. Therefore, we refer to the final kernel function as an “inertial” smoother because its purpose is to generate paths with the physical properties required to obey the laws of classical mechanics. We specify a kernel from the Matérn family of correlation functions 9 with an unknown range parameter, $\phi_{inl}$, and smoothness $\nu=1$ as the inertial smoother. The effect of this final kernel, $h^{(inl)}_{ij}(t,\tau)\equiv\frac{|\tau_{inl}-\tau_{soc}|}{\phi_{inl}}K_{1}\left(|\tau_{inl}-\tau_{soc}|/\phi_{inl}\right)\mathbf{1}_{\left\{i=j\right\}}$, where $K_{1}$ is a modified Bessel function, is visible in the bottom plot of Figure 1. As we did for the social smoothing kernel, in what follows we suppress the superfluous index and write $h^{(inl)}(\cdot)$.

2.3.3 The full three-stage process convolution chain

The expanded specification of the full three-stage PCC is given by

[TABLE]

One of the primary benefits of using a PCC framework is that it allows one to decompose complex dependencies into an iterative sequence of relatively simple mechanisms. Positions of multiple individuals can, in general, exhibit dependencies in both time and among individuals, and the characteristics of these dependencies may be dynamic themselves. We construct a PCC using kernels that compartmentalize the mechanisms generating the complex dependencies in the collective movement of interacting individuals. The first stage of smoothing is in time, and not across individuals, so that $\mu_{i}^{(bm)}(\tau_{bm})$ and $\mu_{j}^{(bm)}(\tau_{bm})$ are independent for $i\neq j$. The second stage of smoothing is across individuals, and not in time, in the sense that $h^{(soc)}_{ij}(t,t^{*})=0$ for $t\neq t^{*}$. The final stage of smoothing is in time only.

Restricting individual smoothing components to operate in only one “dimension” at a time (in this case, either temporal or social) is not necessary. However, a primary motivation for constructing GP models in the PCC framework is to decompose a complex mechanism into components that are easier to understand separately. Therefore, it is natural to form PCC models with kernel functions that operate in a limited number of dimensions simultaneously.

Let $s_{i}(t)$ denote the observed position of individual $i\in\left\{1,\dots,p\right\}$ at time $t\in\mathcal{T}$. We model the true, unobserved position of individual $i$ at time $t$ proportional to $\tilde{\mu}_{i}(t)$ and add a measurement error process, $\varepsilon_{i}(t)$, to yield the data model

[TABLE]

The parameter $\sigma_{\mu}$ scales the random process $\tilde{\boldsymbol{\mu}}(\cdot)$ to account for the overall spatial domain of the locations. We model the measurement error $\varepsilon_{i}(t)$ as i.i.d. mean-zero Gaussian random variables with variance $\sigma_{s}^{2}$, although other measurement error processes could also be used (e.g., 4; Buderman \BOthers. 5).

2.4 Process prior for the dynamic social network

The general model for dependent movement is highly parameterized. To evaluate the covariance $\Delta\tau\widetilde{\mathbf{H}}\widetilde{\mathbf{H}}^{\prime}$ for a grid of $m$ times, one must estimate ${p\choose 2}m$ parameters associated with the social network that drives the second stage of the PCC. As we describe in Section 3, we fit the proposed model to data within a Markov chain Monte Carlo (MCMC) paradigm, which requires repeated inversion of a large, dense covariance matrix. Without further structure for $w_{ij}(\cdot)$, estimating the underlying network would be computationally infeasible. We increase computational efficiency by leveraging reasonable assumptions about the underlying network to constrain the space of possible latent social networks, thereby reducing the effective number of parameters in the model. We impose constraints on the network’s marginal complexity at each fixed time as well as the smoothness of its temporal evolution.

Our approach makes use of recently developed methods for modeling dynamic stochastic networks. To reduce the dimensionality of the latent network at a fixed time, we follow the approaches of Hoff \BOthers. 25, Hoff 24, and Durante \BOthers. 11, who modeled social connections using a latent space. Latent space models for stochastic networks proceed by defining network connections through functions of the locations of particles in a latent $d$-dimensional “social” space. As an example, consider the following latent space model for the edges in a weighted, undirected, dynamic social network. Suppose the latent variables $\boldsymbol{\mu}_{i}^{(w)}(t)$, $\boldsymbol{\mu}_{j}^{(w)}(t)$, and $\boldsymbol{\mu}_{k}^{(w)}(t)$ are positions in $\mathbb{R}^{2}$, and edge weights are defined by $w_{ij}(t)=e^{-d_{ij}(t)}$, where $d_{ij}(t)$ is the Euclidean distance between $\boldsymbol{\mu}^{(w)}_{i}(t)$ and $\boldsymbol{\mu}^{(w)}_{j}(t)$. Defining edge weights based on the distance between particles in a latent space induces positive dependence among the connections $w_{ij}(t)$, because the latent distances are constrained by the characteristics of the latent space and the measure of distance. That is, if $d_{ij}(t)$ and $d_{jk}(t)$ are small, resulting in edge weights $w_{ij}(t)$ and $w_{jk}(t)$ that are close to 1, the triangle inequality implies $d_{ik}(t)$ must also be small, and $w_{ik}(t)$ will also be close to 1. Therefore, latent space approaches for modeling social networks can be understood as a way to induce a tendency to complete triangles in a network. This phenomenon is often observed in human social networks (e.g., friends of friends tend to themselves be friends; Hansell 20), and is also reasonable for many applications to animal social networks. A latent space approach to constraining the underlying social network has the computational advantage of reducing the number of parameters to estimate by a factor of $p$, because we only need to estimate $p$ latent paths instead of ${p\choose 2}$ pairwise relationships.

We define the latent social connection between individuals $i$ and $j$ at time $\tau_{soc}$ using the latent positions and an appropriately chosen functional $g$ by

[TABLE]

In general, the type of network desired for the application (e.g., binary, weighted) may inform the particular choice for $g$. As discussed in Section 2.3, we construct the kernel $h_{ij}^{(soc)}$ under the assumption that the weights in the social network have support $\left[0,1\right]$. The compact support of $w_{ij}(\cdot)$ motivates our choice of $g$, defined by $g(\mathbf{x},\mathbf{y})\equiv e^{-\|\mathbf{x}-\mathbf{y}\|_{2}^{2}}$. The functional $g$ maps two vectors in the latent space to the unit interval $\left[0,1\right]$ and follows the methods used for latent space network modeling by Hoff 24 who showed that this construction induces positive dependence in the edges $w_{ij}(\cdot)$.

[TABLE]

and $d\mathbf{B}^{(w)}_{i}(\cdot)$ are instances of two-dimensional Brownian motion, now used to define the positions in the latent social space. The parameter $\sigma_{w}$ controls the dispersion of the latent paths $\boldsymbol{\mu}_{i}^{(w)}$, and is therefore related to the overall density of the social network, with smaller $\sigma_{w}$ corresponding to higher connectivity. The parameter $\phi^{2}_{w}$ controls the tortuosity of the trajectories through the latent social space, and is related to the temporal stability of the network over time, with larger values of $\phi_{w}^{2}$ corresponding to a more stable network. As described in Section 2.2, we approximate the continuous stochastic processes $\boldsymbol{\mu}_{i}^{(w)}(\cdot)$ using a set of independent normal random variables anchored at a finite number of knots yielding

[TABLE]

3 Model implementation

We obtain realizations from the posterior distribution of the model parameters using a Markov chain Monte Carlo (MCMC) algorithm. In this section, we briefly discuss the most relevant features of model implementation. A more detailed description, including our choice of priors and hyperparameters, is available in Appendicies B.1 and B.2.

We specify independent normal priors for the initial-location parameters and a conjugate hyperprior on the population-level variance so that $\mu_{0i}\sim\mathcal{N}(0,\sigma_{0}^{2})$ and $\sigma_{0}^{2}\sim\text{IG}(a_{0},b_{0})$. Additionally, we integrate the true location process $\tilde{\boldsymbol{\mu}}$ out of the likelihood, which allows us to avoid sampling the continuous movement process directly. This is a common technique used in spatial statistics (e.g., Gelfand \BOthers. 17; Finley \BOthers. 15) to improve mixing for the remaining parameters in the model. The underlying movement process $\tilde{\boldsymbol{\mu}}$ can be recovered post hoc using composition sampling e.g., 15. The resulting integrated model formulation is

[TABLE]

where the parameters $\phi_{inl}$, $\sigma_{0}^{2}$, and $\boldsymbol{\varepsilon}_{w}$ enter the density through the definition of $\widetilde{\mathbf{H}}$ in (15) - (17). To update each latent-space path, $\boldsymbol{\mu}_{i}^{(w)}$, we employ a Metropolis-Hastings algorithm.

4 Simulation

As the dependence among individuals in a population weakens (i.e., $w_{ij}(t)\rightarrow 0$ for all $t$, $i\neq j$), our proposed model simplifies to one that considers each path independently, but shares information about all non-social parameters across individuals. The PCC model constructed from a fixed, empty social network presents a natural baseline for comparison with our full model for dependent movement. Estimating the latent social network that gives rise to the dependence among paths represents the majority of the computational effort to fit our proposed model, and grows rapidly with the number of individuals under study. Thus, there is a natural incentive to model movement under the assumption of path-independence unless it can be shown to be deficient. We compare the performances of the full model for dependent movement (IP-DEP) with the special case of the model under the assumption of inter-path independence (IP-IND). In many cases, we find that accounting for the dependence among individuals results in significantly improved reconstruction of the true underlying paths.

Telemetry devices are subject to a wide range of environmental conditions that frequently result in large intervals of time during which no observed locations are recorded. Estimates of the true locations, $\tilde{\boldsymbol{\mu}}$, during long gaps between observations are often accompanied by large amounts of uncertainty that present a challenge for researchers studying animal behavior. In simulation, we show that the presence of moderate to strong dependence among observed paths provides an opportunity for improvement in path reconstruction when we take into account the joint distribution of all observed individuals. To illustrate the potential gains, we consider a simple setting of two individuals in which the connection status, $w_{12}$, is constant over the period of observation. We assume regular, uninterrupted observations for individual 1, but a gap in time, $\mathcal{T}_{g}$, exists in the observed sequence of telemetry locations for individual 2.

We evaluate the quality of a given model for path reconstruction in terms of both accuracy and precision. To assess the accuracy of a path reconstruction, we define a loss function termed the “squared path error” (SPE) that quantifies the agreement between a specific path reconstruction $\hat{\boldsymbol{\mu}}$ and the true path $\tilde{\boldsymbol{\mu}}_{true}$ as

[TABLE]

We define the reconstruction $\hat{\boldsymbol{\mu}}$ to be the posterior mean for $\tilde{\boldsymbol{\mu}}$ so that (25) represents a measure of how accurately the center of the posterior distribution matches the true underlying path. To assess the precision associated with a model, we compute the radii for circular 95% credible regions surrounding each point in the path $\hat{\boldsymbol{\mu}}(t_{g})$ and average across the entire temporal gap in observations, $\mathcal{T}_{g}$, to yield an overall summary of precision we term the average circular 95% credible region radius ($\text{ACCRR}_{0.95}$).

To provide a general sense of when accounting for path dependence results in the greatest gains in path reconstruction, we vary the strength of dependence, $w_{12}$, the proportion of the study interval made up by the observation gap, $|\mathcal{T}_{g}|/|\mathcal{T}|$, and the tortuosity of the paths, and use SPE in tandem with $\text{ACCRR}_{0.95}$ to compare the performance of IP-DEP and IP-IND. We vary tortuosity through the range parameter $\phi_{inl}$ that appears in the inertial smoother, with smaller values corresponding to more tortuous paths on average. For each combination of $w_{12}$, $\mathcal{T}_{g}$, and tortuosity, we used the five step procedure:

Simulate a realization of $\tilde{\boldsymbol{\mu}}_{true}$ and $\mathbf{s}$ from our proposed model (see Figure 2). 2. 2.

Fit both the dependent and independent-paths models to $\mathbf{s}$. 3. 3.

Sample 1000 paths from the posterior distribution of $\boldsymbol{\mu}$ using composition sampling. 4. 4.

Compute the SPE defined in (25) for the posterior mean of $\boldsymbol{\mu}_{2}$. 5. 5.

Use the same 1000 draws from the posterior to compute the $\text{ACCRR}_{0.95}$ for $\boldsymbol{\mu}_{2}$.

We repeat the procedure 20 times for each combination of parameters to obtain an estimate of the variability among realizations of the simulated paths. Values for the parameters used to simulate each path as well as all prior distributions and hyper parameters are provided in Table 1, and an example simulation is shown in Figure 2.

To facilitate direct comparison between the two models under consideration, we examined the ratios of SPE and $\text{ACCRR}_{0.95}$ for IP-DEP and IP-IND, defining ratios such that the relevant value for IP-DEP appears in the denominator and values greater than 1 show support for IP-DEP. Figure 3 displays the ratios of SPE (top row) and $\text{ACCRR}_{0.95}$ (bottom row) under all combinations of values for the parameters $w_{12}$, $|\mathcal{T}_{g}|/|\mathcal{T}|$, and $\phi_{inl}$ (tortuosity). Individual plots show the median ratio across all simulations for high (solid) and low (dashed) tortuosities, with an associated polygon delineating the 25% and 75% quantiles across simulations. Columns organize the plots by gap size, and the strength of the social connection $w_{12}$ increases along the x-axis within each individual plot.

Several general observations can be made about the circumstances under which fitting the full model for dependent movement offers the greatest improvements in path reconstruction. First, both precision and accuracy improve near monotonically with increasing $w_{12}$. Second, the greatest gains come when the gap size is moderate to large (columns 3 and 4). When the gap is brief (column 1), the difference between the reconstructed paths for each model is modest. Finally, the gains in performance for the full model are greater for the case of high tortuosity. In our simulation study, more tortuous paths are characterized by a shorter range of dependence in time; thus, increasing tortuosity is similar to increasing the size of the observation gap.

5 Killer whales

We analyzed telemetry data for four killer whales near the Antarctic Peninsula (see Figure 4) over the course of five days in February 2014. Geographic positions were measured using Argos telemetry tags (for a complete description of the tags and study area see Andrews \BOthers. 1; Durban \BBA Pitman 13). Although multiple types of killer whales have been described in this area, all four tags were deployed on individuals from the same population of the most common Type B2 killer whales 12.

It is immediately apparent from the data that individuals 1 and 2 (Figure 4, bottom left) show potential evidence of a close connection. In addition, there is some ambiguity about the relationship between individuals 3 and 4 (Figure 4 top right), because they occupy approximately the same spatial region during the study period. In contrast, there is little reason to suspect dependence between the two pairs of individuals. We analyzed the movement of all four individuals jointly, which allowed for pooling of information about measurement error across the entire group, and provided an opportunity for basic model validation. By fitting a model for the joint movement of all four individuals with a fully flexible latent network structure, we were able to check for the presence of potentially spurious network connections, because existing knowledge about the system suggests that the underlying network should exclude connections between the two subgroups.

The top plot of Figure 5 shows the observation times for the killer whales in our study, with darker regions corresponding to a denser rate of telemetry measurements in time. Two different day-long gaps in observation occur for individual 1 (Figure 5), a consequence of the original study design used in the deployment of the telemetry tags which sought to balance the need for temporally dense observations with limitations in the battery life of the tags by collecting measurement only every other day for select individuals. The observation times for the other three individuals cover these gaps, suggesting that modeling the four paths jointly may allow for more precise and/or accurate estimates of the true path taken by individual 1 if moderate to strong dependence exists between it and any of the other three whales.

As in the simulation study (Section 4), we investigated path reconstructions generated by our proposed full model for dependent movement (IP-DEP), as well as the analogous model under an assumption of inter-path independence (IP-IND). We found that the full model generates path reconstructions in which the uncertainty about the position of individual 1 is dramatically reduced, compared to the reconstructions generated under the assumption of inter-path independence as measured by 95% circular credible interval radii (bottom plot of Figure 5). Circular credible regions for the true position of the individual during the two large gaps in observation occuring on February 12th and 14th for the IP-IND model have large radii, sometimes exceeding 70km. In contrast, circular credible regions for the true position based on the IP-DEP model are less than a third this size, and similar in magnitude to the uncertainties for individuals for which we have dense observations. Despite the uncertainty about the latent social network, the IP-DEP model offers a substantial reduction in uncertainty about the true path taken by individual 1.

The posterior distributions for the network relationships are shown in Figure 6. We can see strong evidence for a significant relationship between individuals 1 and 2, however there is weak evidence of a meaningful social connection between individuals 3 and 4. Similarly, there is no evidence of connections existing between any other pairs of individuals. Credible intervals for other model parameters, as well as the specified prior distributions and hyperparameter values are provided in Table 2.

6 Discussion

Appropriate models for dependent, multivariate data are application-dependent. Gaussian processes offer a flexible, parsimonious tool for analyzing complex data, such as those that arise from measuring the movement of animals, if one can specify a valid form for the covariance. The need for covariances that are mathematically sound, realistic, and interpretable has motivated decades of research with a vast array of applications because satisfying all three of these goals simulataneously is challenging. Our proposed PCC framework represents a novel perspective in the construction of covariance functions. Through an application involving the movement of killer whales, we demonstrated that the PCC framework can be used to create sophisticated covariance functions that account for several important mechanisms, without resorting to unrealistic assumptions such as separability.

The joint movement of interacting individuals can be viewed as a multivariate temporal process. Methods for interpolating multivariate spatial and spatio-temporal processes, also called “cokriging,” have been studied for several decades (e.g., Myers 32; Cressie 9). The primary challenge has been to develop models for multivariate processes that accurately capture the dynamics within each process, as well as across processes, while ensuring the resulting cross-covariance structure is valid (symmetric and non-negative definite). While not explicitly framed as a PCC, 41 proposed a model for cokriging that shares important connections with our proposed model for the joint movement of interacting animals. 41 also approached the problem of specifying a valid, mechanistically motivated cross-covariance structure for multivariate processes through the use of kernel convolutions, relaxing unrealistic assumptions about the covariance structure. Similar to the way dependence among individuals arises through the application of a social smoothing kernel (Section 2.3.1), 41 induced dependence in a multivariate process by smoothing across variables. The cross-covariance described by 41 also accommodates spatio-temporal lags in the dependence among variables. This generalization may also be a useful feature in future models for animal movement, where a temporal lag may allow researchers to capture the effect of animals following one another.

The particular model we constructed for animal movement cannot be used to understand all forms of dependence that may exist among interacting individuals. Rather, it has been tailored to the case in which interactions among animals manifest themselves as movement along proximally close trajectories. Particles may also exhibit other forms of dependence, such as a tendency to repel in the case of strong territorial behavior, that would not be well-described by the same covariance function. The PCC framework allows researchers to define and order kernels as necessary to appropriately model the mechanisms under study. Under a different ordering, the same three kernels we employed result in another useful covariance function for the study of particle movement. We briefly discuss this alternative to highlight the flexibility offered by the PCC framework and provide some intuition for the role certain kernel functions play in the characteristics of the resulting random process.

As mentioned in Section 2.1, choosing the ordering of the kernels in a PCC represents a meaningful modeling decision when kernel convolutions do not commute. In our proposed model for dependent movement, the second and third stages of smoothing commute, while the first and second do not. We can interpret the commutativity between the social and inertial smoothers mechanistically, by noting that it does not matter whether the effect of the latent social connections operates on the raw, Brownian paths, or the temporally smoothed paths. In contrast, if we socially smoothed white noise and then used the kernel associated with Brownian motion, a very different form of dependence results. This alternative ordering ($h^{(soc)}\rightarrow h^{(bm)}\rightarrow h^{(inl)}$) presents another plausible mechanism through which dependence might arise in paths taken by multiple particles.

When applied to the velocity of a particle, smoothing with the step function kernel $h^{(bm)}$ returns the associated position process (after taking into consideration the appropriate initial location of the particle). Therefore, Brownian motion can be thought of as a random position process in which the velocity of the particle during each infinitesimal span of time is a realization from a Gaussian white noise process. In our proposed PCC model for dependent movement, the social smoother has the effect of “shrinking” the positions of strongly-connected particles toward each other. If we instead employed the social smoother before the step function kernel, the result would be to “shrink” the velocities of the particles together, rather than their positions. Thus, we would be inducing a tendency for particles to move in similar directions, though not necessarily a tendency to be in similar locations. Such an effect would be visible in the particles as movement in parallel, perhaps with a considerable distance between connected particles.

A PCC approach to constructing covariance functions for GPs allows for broad flexibility in model development, offering researchers a highly customizable framework that can be used in a wide variety of applications. We demonstrated the value of this approach with an application to animal movement, however PCCs can be used to model a broad range of random processes. Rather than relying on parametric families of covariance functions, PCCs encourage the use of interpretable, intuitive, and problem-specific convolution kernels, allowing for direct incorporation of scientific knowledge.

acknowledgements

Killer whale tagging was conducted under permit #14097 from the National Marine Fisheries Service and Antarctic Conservation Act permit #2009-013. Shipboard tagging operations were supported by Lindblad Expeditions and the National Geographic Society, and by an NSF rapid grant to Ari Friedlaender. Robert Pitman helped with tag deployments and identification of killer whale types in the field.

This research was supported by 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. The findings and conclusions of the NOAA authors in the paper are their own and do not necessarily represent the views of the National Marine Fisheries Service, NOAA.

Appendix A Model specification

Data model

Each orthogonal spatial direction (longitude and latitude) is modeled independently using the same specification.

[TABLE]

Process model (integrated out for model fitting)

[TABLE]

Prior model

[TABLE]

Application only:

[TABLE]

Fixed parameters

Simulation study:

[TABLE]

Application:

[TABLE]

Appendix B Model fitting details

B.1 Simulation Study

We simulated approximately continuous true paths using a grid of 500 equally spaced time points on the unit interval. We simulated 100 observation times drawn uniformly from the same unit interval. For estimation, we fixed the hyperparameters associated with the network to $\phi_{w}=4/15$ and $\sigma_{w}^{2}=10$, and used the priors given in Table 1. For each fit, we obtained 4,000 iterations and discard the first 3,000 as burnin. The entire simulation study was parallelized across 6 3GHz cores and required approximately one week of computation time.

B.2 Killer whales

We acquired 100,000 iterations on a single computing node, and discared the first 50,000 as burnin to yield a sample size of 50,000. We employ diffuse priors for all parameters. We set the hyperparameters associated with the network to $\phi_{w}=0.3$ and $\sigma_{w}^{2}=10$. Model fitting was performed using a processor speed of 3 GHz and required approximately 100 hours of computing time.