## 1 Introduction

For decades, study of animal movement has been an important source of scientific understanding and discovery in ecology (Hooten et al., 2017). Observations of animal trajectories using telemetry devices, such as radio collars, have provided researchers with information about the way animals interact with their environment (e.g., Manly et al. 2002; Johnson et al. 2008) and each other (e.g., Niu et al., 2016; Scharf et al., 2016; Scharf et al., 2018). In the case of the former, one of the most common approaches used by ecologists is to analyze the rates at which animals use certain types of habitats relative to the distribution of habitat types available to them (Manly et al. 2002; Lele and Keim 2006; Hooten et al. 2017

). Such analyses typically make use of a so-called “use-availability” framework in which the probability of an individual using a particular location is modeled as a weighted combination of all available locations

(e.g., Johnson et al., 2006; Northrup et al., 2013). The weights are then referred to as the resource selection function (RSF), and provide insight into which portions of a landscape are most valuable to the study species (Boyce and Mcdonald, 1999). What constitutes an available location depends on the characteristics of the particular species under study and the rate at which telemetry observations are gathered. Typically though, the degree to which a location is available is a function of how far an animal can reasonably be expected to move between observation times, and/or the size of the home range of the individual (e.g., Christ et al. 2008; Brost et al. 2015).For many species, specific geographical features in the landscape can have a strong effect on where individuals choose to move. Such features are sometimes well-summarized by a single point (e.g., dens or kill sites), but may also correspond to higher-dimensional subspaces (e.g., rivers or lakes). Their locations may be relatively static in time (e.g., coastlines or home range centers (Brost et al., 2016)), or may be dynamic (e.g., sea ice extent or areas of high-quality forage for herbivores). While not always framed in the context of resource selection, these landscape features can nevertheless be thought of as resources, and the behavior of animals may demonstrate selection for (or against) points near the feature. We introduce a novel suite of models for the analysis of animal movement that incorporates active selection for features in a landscape that may have complex and dynamic shapes.

Our modeling framework is motivated by the study of polar bear (Ursus maritimus) movement. Polar bears spend much of their time on sea ice over shallow, biologically productive water where they hunt seals. During the sea ice melt season, the part of the year when sea ice first breaks up and contracts toward the pole, then freezes and expands southward again, the most valuable habitat is near the interface between sea ice and the ocean (Durner et al. 2009; Rode et al. 2015; Atwood et al. 2016). The changing distribution and characteristics of sea ice throughout the late spring through early fall means that the location of valuable ice-edge habitat is constantly shifting. As climate change alters the rate at which sea ice thaws and freezes, as well as the size of its minimum and maximum extents, there is increasing concern about how polar bears are responding to these dramatic shifts in their environment (Rode et al., 2014). Our goal is to develop a model for the movement of polar bears that explicitly accounts for the effect of the changing sea ice and can be incorporated into a wide variety of hierarchical models used to better understand polar bear ecology. In Section 3, we use our proposed model to answer a particular question posed by wildlife managers who seek to cluster polar bears from the Beaufort and Chukchi seas into two sub-populations.

## 2 Model Development

### 2.1 Feature preference

To account for an individual’s preference for areas in a landscape near (or far from) a particular feature of interest, we take a use-availability approach. We define the selection weight of a particular location to be a parametric function of the distance from the location to a feature of interest on the landscape. Estimates of the relevant parameters provide a summary of how strong an effect the feature has on the behavior of the observed individuals. We model availability similar to

Hjermann (2000), Christ et al. (2008), Johnson et al. (2008), and Brost et al. (2015) who used radial distributions centered on the most recently-observed location to define the continuously-valued availability at each point in time.Let be the location of an individual at time , where we will assume for now that observation times are equally spaced with no missing values. Define the conditional probability density for as

(1) |

where we use square brackets to denote a probability density. For

, the conditional distribution is proportional to the product of two components, the first of which is the density of a bivariate Gaussian distribution centered on the previous location of the individual, and defines the availability of each point on the landscape as in

Christ et al. (2008). The availability component induces positive auto-correlation in the joint process , with larger values of resulting in processes with greater distances between consecutive locations and faster, more erratic movement.The function is a RSF and controls the effect a particular feature in the landscape has on an individual’s movement. Let denote the set of points that comprise the feature of interest (e.g., the interface between sea ice and ocean). We define the function as

(2) |

(where is the or Euclidean norm) so that the value of is highest near , and drops to zero as moves away from . The value of controls the range at which effectively reduces to zero. We show in Section 2.2 why this particular parametric form for the RSF can be used to leverage computational efficiencies in parameter estimation.

In practice, it will often be useful to define such that it achieves its largest values at locations near so that the conditional density given in (1) has probability mass concentrated near . Specifying in this way provides a method for modeling movement that exhibits preference for the region of the landscape near the feature of interest. Alternative specifications could also be used to model the movement of individuals displaying preference for portions of the landscape far from the feature .

The model for the discrete-time process provides a useful tool for modeling the movement of an individual responding to a one-dimensional feature on a landscape. In Section 3, we apply the model to the movement of polar bears with the ultimate goal of clustering individuals into disjoint sub-populations based on space use. By including availability and resource selection as part of a larger hierarchical structure, we are able to account for polar bears preference for habitats that facilitate the depredation of seals, which, if ignored, might result in biased inference about sub-population membership.

The conditional density in (1) is only defined up to a constant of proportionality that must be computed as part of any likelihood-based estimation procedure. In Section 3

, we employ a Bayesian hierarchical methodology and fit our model for polar bear movement using Markov chain Monte Carlo (MCMC), which requires computation of the normalizing constant several times at each iteration of the algorithm. For a general feature,

, the normalization constant is not analytically tractable. Thus, some form of numerical integration is required to fit the model to data, the computational cost of which precludes such an approach in our application. In the RSF literature, the normalization constant has typically been approximated using either a coarse spatial discretization (e.g., Warton and Shepherd 2010; Brost et al. 2015), or a randomized scheme based on an “availability sample” (e.g., Northrup et al., 2013). In the next section, we introduce a novel approach for computing the necessary normalization constant in which we approximate the RSF in a way that induces conjugacy in the distributional form for and greatly reduces the computational cost of model fitting.### 2.2 Linearization approximation

We implement a novel approximation technique that assumes locally linear structure in the shape of , allowing for efficient approximation of the true conditional density of . To motivate our approximation, we note that, for the special case when is a straight line, the RSF as defined in (2) can be written in a form similar to that of a bivariate Gaussian density function with a rank-deficient covariance matrix.

#### 2.2.1 RSF for straight lines

First, consider the case of a vertical line, , in the real plane so that . For this case, we have

for all real-valued .

To allow for that are straight, but not necessarily vertically oriented, we rotate the coordinate system through an angle . Let be the rotation matrix defined as

and let be defined such that for some real-valued . Note that the inverse of a rotation matrix, , is also equal to its transpose, . The RSF defined in (2) is invariant under rigid transformations such as rotations, therefore

#### 2.2.2 Linearizing complex landscape features

The resulting form for is proportional to the product of two Gaussian distributions, one of which is improper. Provided , the product is a proper bivariate Gaussian distribution with mean, , and covariance, , given by

where is any point in . The distributional form of implicitly defines the appropriate normalization constant in (1). Thus, if there exists some straight line that represents a close approximation to near , may provide a reasonable approximation for that alleviates the computational burden of repeatedly calculating the necessary normalization constant. A natural candidate for is the line that is tangent to at the point on closest to , because this is the portion of the feature that contributes most to the conditional distribution of .

Let be the point in nearest to , and let be the set of points that lie on the line tangent to at . Figure 1 shows a schematic illustrating the way the product of the two Gaussian components in (1) combine to result in a third Gaussian distribution. The dashed line represents an example of the edge of a one-dimensional feature of interest, and the solid line represents the linearized edge.

For to result in an adequate approximation of the RSF, the linearized feature need only resemble the true feature in the vicinity of . Outside of the immediate neighborhood, the availability distribution will be essentially zero, reducing the impact of errors in the RSF approximation. Therefore, the linearization approximation is justifiable if it is reasonable to assume that, at the scale at which locations are measured, the edge of the feature of interest can be considered approximately linear.

## 3 Application

### 3.1 Goals and previous work

There are a total of 19 recognized polar bear sub-populations in the circumpolar Arctic (Figure 2, Obbard et al., 2010). However, the boundaries that delineate the sub-populations are challenging to precisely define because there are few barriers to movement for polar bears, and the changing extent and drift of the sea ice leads to periods of the year when individuals from different sub-populations may use overlapping portions of the landscape. Nevertheless, there are important reasons to determine a clear delineation of the sub-population boundaries. For example, wildlife management agencies such as the U.S. Fish and Wildlife Service (USFWS) use sub-population boundaries to help guide management decisions for polar bears, which are currently listed as ‘threatened’ under the Endangered Species Act (U.S. Fish and Wildlife, 2016). There is also evidence that polar bears from different sub-populations are responding to climate change with differing degrees of success (Rode et al., 2014; Ware et al., 2017). In what follows, we focus on estimating a sub-population boundary between the Chukchi Sea (CS) and Southern Beaufort Sea (SB) sub-populations.

Previous efforts to determine sub-population boundaries provided important initial estimates; however, the substantial physical changes taking place in the Arctic ecosystem mean that boundary estimates have the potential to become outdated quickly. Amstrup et al. (2005) used a multi-stage non-parametric clustering procedure to estimate boundaries based on utilization distributions estimated from telemetry observations. While the methods employed provide a probabilistic basis for inference, there are three primary opportunities for continued development.

First, Amstrup et al. (2005) did not make use of all observations. They discarded observations with measurement error classes deemed insufficiently precise, and thinned observations so that only one location is used per bear per 6-day duty cycle. Thinning the data represents an effort to diminish the dependence between consecutive observations so that locations can be treated as independent in the analysis. In our analysis, we directly incorporated the dependence between consecutive observations. We also modeled the wide variety of measurement error structures and severities present in the data allowing us to include all observations.

Second, the Amstrup et al. (2005) analysis did not provide probabilistic estimates of sub-population membership for each observed individual. Our analysis allowed us to compute the probability of sub-population membership for each polar bear, with estimates of uncertainty. Moreover, our methodology allowed us to estimate the sub-population membership probabilities for a newly observed individual, conditioned on the data used to fit the model.

Finally, Amstrup et al. (2005) did not directly incorporate the effect changing sea ice had on the observed locations of polar bears. Estimates of space use for each sub-population were summarized across all seasons, potentially confounding the roles of sea ice and sub-population boundaries. In contrast, we directly incorporated the effect of sea ice extent on movement through a RSF. Our novel linearization approximation allowed us to obtain inference from the resulting mechanistic model using an efficient sampling algorithm.

### 3.2 Movement process model

With the exception of mothers denning on land, most polar bears remain on the ice in the summer, following it north as it retreats away from continental land masses. Polar bears are specialist carnivores that use areas near the interface of sea ice and ocean to hunt seals during the sea ice melt season (Durner et al., 2009). However, polar bears do occasionally remain on land through the warmest part of the year and, in fact, Atwood et al. (2016) concluded that there is evidence the number of individuals exhibiting this behavior is increasing.

Two of the primary characteristics of polar bear movement are a tendency for most individuals to prefer portions of the landscape near the edge of the sea ice, and a tendency for individuals to occupy a general spatial region corresponding to the particular unobserved sub-population to which an individual polar bear is a member. To address the effect of the sea ice, we specify a model for movement that incorporates preference for areas near the sea ice boundary as it changes in time. We account for polar bears that remain on land during the summer by also including coastline as a feature on the landscape associated with increased rates of use, as individuals that spend the summer on land tend to remain on islands and/or near the coastline for most of the season (Rode et al., 2015). By accounting for habitat variability, we hope to avoid inadvertently clustering bears into sub-populations that confound true sub-population spatial regions with movement responding to changing sea ice.

We take a Bayesian hierarchical modeling perspective, allowing us to specify models for both measurement error and the true unobserved movement process in a single coherent framework. We use a discrete-time approach to specify a model for polar bear movement, in which the conditional probability density for the location of an individual is proportional to a product of components corresponding to resource availability and selection. The discrete-time approach is motivated by the scale at which we have measurements of Arctic sea ice, which features as a direct effect in our model for movement. We use estimates of the extent of sea ice provided by the National Snow and Ice Data Center (Fetterer et al., 2010) that are available at daily intervals; thus, we model movement as a discrete-time stochastic process on a daily scale.

Let denote the location of individual at time , where is a consecutive set of times from the reference set , and let

be a binary random variable that equals 1 if individual

belongs to the CS sub-population, and 0 if it belongs to the SB sub-population. Additionally, we denote by the set of points in the plane defined by the union of coast line, and the edge of the sea ice. We model the conditional distribution of as(3) |

Each component in (3), captures a different feature of the movement process. Namely, these are (i) the temporal dependence between locations on consecutive days, (ii) the association of each individual bear with a sub-population-level central place, and (iii) a RSF that appropriately weights locations near a coastline or the edge of the sea ice. The first two terms can be thought of as a two-component availability function that incorporates movement constraints and a sub-population activity center, similar in many respects to the modeling specification of Christ et al. (2008) and Johnson et al. (2008). The third term is a RSF that models the preference polar bears exhibit for habitat near either a coastline, or the sea ice boundary.

The model for movement specified in (3) also captures a documented secondary effect that changing sea ice has on polar bears (Durner et al., 2017). Especially in the summer months, polar bears that remain on the sea ice expend a substantial amount of energy merely keeping up with the retreating ice sheet. Increased rates of sea ice retreat associated with global climate change have been estimated to impose an energy cost on polar bears of 1-3 additional seals per year compared to historic norms (Durner et al., 2017). Sea ice does not retreat at a uniform rate at all points along the boundary. Rather, the change in shape is highly variable with some regions retreating more slowly than others (Steele et al., 2015). Our proposed RSF-based movement model has the effect of placing higher probability mass on paths that track with regions of the ice that have the slowest daily rates of change. Our model therefore captures an inherent incentive for polar bears to travel along those portions of the ice that minimize their expected energy expenditure associated with keeping pace with the changing ice.

### 3.3 Activity centers and sub-population membership

The second component in (3), for , and first component for is a bivariate Gaussian distribution centered on one of two central places, and , corresponding to the centers of the SB and CS sub-populations, respectively. The covariance matrices, and

, control the strength of the effect the sub-population center has on the movement of each individual. As the marginal variances of

and increase, individuals are allowed to range farther from their central place. We specify prior distributions for and based on their eigen decompositions. Details are provided in the Supplementary Materials.The binary random variables, , indicate the particular sub-population with which individual associates. We specify Bernoulli prior distributions for each with probability 0.5, corresponding to balanced a priori classification of each individual.

### 3.4 Time-varying RSF

For the application to polar bear movement, we modified the RSF introduced in (2) slightly to account for an important aspect of polar bear ecology. Polar bears’ preference for habitat near the coastlines and boundary between ocean and sea ice is highest during summer and autumn months and weakest in the winter (Durner et al., 2009). Thus, we define a RSF with the same shape characteristics as in (2), but with an added indicator function that accounts for the difference in seasons when the landscape feature is most relevant, and when it has negligible impact on movement. Let and denote the days of the year on which this “summer” season begins and ends, respectively. Then the modified RSF is defined as

(4) |

We specify Gaussian prior distributions for the end-points of the summer such that

corresponding to prior start and end date distributions with means of May 15th and November 15th, respectively, and a variance of two weeks.

One way to interpret this generalization of (2) is to define the RSF in the original way, but let vary in time. The modified RSF given by (4) arises for the case of a dynamic parameter that is constant and finite during the summer season, and infinite during the rest of the year. Thus, outside of the season defined by the landscape feature defined by has no effect on movement.

### 3.5 Measurement error

We analyzed telemetry observations of polar bears made by the USFWS (USFWS unpublished data) and the U.S. Geological Survey (USGS) (USGS unpublished data) using a variety of different telemetry device types contaminated with measurement error of varying severity. In each case, we model the observed locations as centered on the true, unobserved location of the individual with an additive measurement error process as

where the distribution for depends on the particular device used to make the observation , and we allow measurements to be made at any point on the continuous interval . We address the misalignment between this continuous-time scale and the discrete scale used to model the movement process in Section 3.6.1.

The simplest of the three measurement error structures is that associated with GPS devices, which generate small measurement errors that are well-modeled by a circular Gaussian distribution as

The second and third types of measurement errors are those associated with Argos transmitters (Service Argos, 2016). Argos devices use polar orbiting satellites to record the location of the device and have the advantage of providing values to researchers in real time. However, the satellite-based location estimates are often made with a substantial amount of error, and the structure of the errors can be highly non-circular (Costa et al., 2010). The severity of Argos measurement errors are related to the direction of satellite travel, leading to elliptical error distributions with an orientation that varies with each observation. Service Argos provide categorical error-class labels (3, 2, 1, 0, A, B, or Z) that correspond to measurements made with increasing severity levels. Recently, Argos has begun to include estimates of the maximum and minimum axis lengths, and angle of rotation for the error ellipse associated with each observation (e.g., McClintock et al., 2015).

For Argos observations that lack auxiliary ellipse information, we use the provided classes to model the measurement error process as

where indexes unique variance parameters for each Argos error class. We specified the means in the log-normal prior distributions for to match estimates obtained by Costa et al. (2010), which are provided in the Supplementary Materials. For an alternative treatment of Argos error structure in the absence of auxiliary ellipse information, see Brost et al. (2016) and Buderman et al. (2016).

For Argos observations with auxiliary ellipse information, we model the measurement error process as

where is a covariance matrix corresponding to a bivariate Gaussian distribution with elliptical contours that match the ratio of axis lengths and angle of rotation provided by Argos for each observation time.

### 3.6 Estimation

#### 3.6.1 Process imputation: Stage 1

Our proposed model for the movement and measurement error process are based on discrete- and continuous-time scales, respectively. To reconcile this difference, we employ a two-stage estimation procedure called “process imputation” (e.g., Hooten et al., 2010; Hanks et al., 2011; Scharf et al., 2016; Scharf et al., 2017). The procedure is based on the well-established method of multiple imputation for missing values (Rubin, 1996), and fits within a larger MCMC algorithm. In essence, it consists of first fitting a flexible, continuous-time model to the telemetry data and making a finite number of draws from the posterior distribution of the true continuous process at times . The realizations from the first stage are selected uniformly at random and treated as the true discrete-time paths in the second stage of the process imputation procedure. Sampling uniformly from the collection of first-stage realizations propagates the measurement error incorporated in the first stage through to the second stage.

Using the process imputation procedure allows us to account for misaligned time scales in the measurement and movement processes, account for measurement error, and ease the computational demand of model fitting in exchange for a small amount of approximation error in the overall estimation procedure. We refer the reader to Scharf et al. (2017) for a more detailed treatment of the procedure, and justification for its use.

We used the Ornstein-Uhlenbeck (OU) movement model of Johnson et al. (2008), implemented in the `R`

package `crawl`

(Johnson, 2016), as the first stage of a process imputation estimation procedure. Johnson et al. (2008) modeled the velocity of an individual using Brownian motion, integrating the velocity process over time to yield the movement process. The OU model has been shown to perform well as a first-stage model for process imputation (Scharf
et al., 2017), and accommodates a broad range of measurement error models. In stage one, we fit each individual separately, using the measurement error structures defined in Section 3.5 as appropriate for each individual’s device type, and drew realizations of the continuous path evaluated on the daily discrete-time scale for all days falling between the first and last observation times. Additional details associated with model fitting are provided in the Supplementary Materials.

#### 3.6.2 Process imputation: Stage 2

Conditioned on the true discrete-time locations, , we estimate the posterior density, , for model parameters with a MCMC algorithm, using standard Metropolis within Gibbs samplers for each parameter. We used visual examination of trace plots, as well as estimates of effective sample size to assess convergence. A single chain of length was generated, with the first half discarded as burnin. Additional model fitting details may be found in the Supplementary Materials.

#### 3.6.3 Pre-computation

In addition to providing a method for reconciling the time domains at which observations are made and the movement process is modeled, the process imputation procedure provides opportunities for improvements in computational efficiency. In the second stage of our fitting procedure, all parameter updates are conditioned on one of known draws from the distribution of true paths generated in stage 1. Hence, there are only possible values of , the location on closest to , which we can compute and store before running the MCMC algorithm. Computing requires evaluating the distance between and every point in , and must be done for all individuals and all time points. Pre-computing all possible values of reduces the number of times we must perform the complete search from (the number of iterations used in the MCMC algorithm) to (the number of imputed paths) resulting in a substantial decrease in the time required for estimation.

### 3.7 Results

Using the two-stage approach, we fit the model to all observations made by USFWS and USGS from 2012–2016 (186 unique individuals). Table 1

gives the posterior medians and equal-tailed 95% credible intervals for each variance parameter, as well as relevant prior distributions and hyperpriors. Figure

3 shows the posterior mean of each grouped by the agency that tagged the individuals. Posterior means of the class indicator variables can be interpreted as the posterior probabiltiy that individual is a member of the Chukchi sea sub-population.posterior summary | prior | ||

parameter | median | (2.5%, 97.5%) | density |

220 | (218, 223) | ||

7300 | (6890, 7920) | ||

see Figure 3 | |||

62.9 | (60.2, 64.0) | ||

334.1 | (332.1, 348.9) |

The posterior distribution of the parameters defining the “summer” season suggests that the feature , defined as the union of the interface between sea-ice and ocean and the continental coastlines, has greatest impact on the behavior of polar bears between March 4 and November 30.

To produce a meaningful spatial delineation of the two sub-populations from which our study animals were drawn, we use a derived quantity related to the infered locations of the sub-population activity centers. If we consider the observation of a single new location,

, and integrate across all arbitrary features, effectively removing the effect of the RSF, it can be shown that the posterior probability that

is given byFor the case of an uninformative prior, , the contour corresponding to is defined by the points in the plane where the two normal densities are equal. It can be shown that these contours are the roots of quadratic polynomials in two dimensions. By computing the contour at each iteration in the MCMC algorithm, we can obtain draws from the posterior distribution, a summary of which provides wildlife managers with a way to delineate the boundary between the two sub-populations.

Figure 4 shows a map of the region encompassing the Chukchi and southern Beaufort seas. In the background of Figure 4, weekly measurements of sea ice extent for March-September 2016 are shown as light blue polygons, with the darkest polygon corresponding to open ocean, the second darkest polygon to March 1st, and the lightest polygon to September 30th. Orange and purple lines show the paths drawn in the first stage of the process imputation framework. The colors correspond to the posterior means of , with dark orange hues corresponding to values close to 0, white corresponding to values close to 0.5, and dark purple corresponding to values close to 1. The solid black line corresponds to a central-measure summary of the posterior distribution of the derived spatial boundary. The dashed lines show equal-tailed 95% pointwise credible intervals computed orthogonal to the solid line (more details are given in the Supplementary Materials). The portion of the inferred boundary most relevant in this application is in the bottom right quarter of the map and suggests at most a small shift from the currently accepted sub-population delineation, denoted by the large polygons with thin black lines (also shown in Figure 2).

## 4 Simulation Study

The linearization approximation implemented in the estimation procedure (Section 2.2) introduces approximation error into the likelihood calculation unless the shape of the linear boundary is exactly a straight line. In many applications, the shape of a particular feature may be quite complex at the scale of the entire landscape, but well-approximated by a linear shape at the scale of movement increments. If the boundary of the feature is straight enough in the relevant vicinity of the individual reacting to it, then it is reasonable to assume that the approximation error will not adversely affect estimation.

To check the assumption of local linearity, we investigated differences between data simulated from the exact and approximated movement models. We used the same observed union of sea ice boundary and coastline from 2016 that we used in our analysis of the movement of polar bears (Section 3), and the posterior medians of model parameters and to simulate 16 paths over all 366 days from both the exact and approximate movement process models. For the purposes of the simulation study, we disregarded the effect of sub-population activity centers, effectively taking the limit of the movement model in Section 3 as the marginal variances in and tend to infinity. We then compared the distributions of three summary statistics: (1) the distance between each simulated location and the nearest point on , (2) the distance between consecutive locations, and (3) the turning angle between consecutive locations. The first statistic was chosen to reveal important discrepencies in the way the exact and linearized models account for the effect of the RSF. The second and third statistics are commonly investigated quantities in animal movement studies.

Figure 5

shows kernel density estimates of the distributions for each statistic. Each of the narrow, faint lines represents the distribution of a statistic for a single individual. The thick, dark lines represent distributions taken across all simulated individuals. The strong similarity between the solid and dashed lines suggest that our linearization method is providing an adequate approximation to the true likelihood.

## 5 Discussion

Our novel approximation method based on the linearization of a potentially time-varying spatial feature, , reduces the computational burden of fitting models in the common RSF framework, allowing researchers increased flexibility in the types of RSFs they can specify in mechanistically-driven models for movement. We demonstrated our approach in an application involving the movement of polar bears as they respond to seasonal shifts in sea ice over the course of 2012–2016. Recently, there has been a significant increase in research focused on the so-called “greenwave” hypothesis (e.g., Aikens et al., 2017), which posits that herbivorous animals align their movement with bands of high-quality forage that shift throughout spring as different elevations and latitudes experience phenological changes. Our modeling approach represents a way to validate the hypothesis if information about the shape of the greenwave is known, or potentially estimate the location of the posited band of high quality forage based on the observed movement patterns of herbivores.

The linearization approximation methodology can also be extended to higher-dimensional spaces and features. For instance, in marine environments, RSFs based on two-dimensional features, such as isotherms, may be locally approximated using rank-deficient Gaussian distributions corresponding to infinite planes. One-dimensional features in three-dimensional spaces, such as wind or ocean currents, can also be approximated with improper distributions, where the rank of the covariance matrix is deficient by a degree of 2. In practice, researchers will need to evaluate the appropriateness of a linearization approximation; however, in many cases, our methodology offers a way to include complex drivers of movement that might otherwise have been computationally inaccessible.

## 6 Suplementary Materials

## Appendix A: MCMC Implementation details

### Stage 1

We used the Ornstein-Uhlenbeck movement model of Johnson
et al. (2008) as the process model for the imputation distrbution in the process imputation framework (Scharf
et al., 2017). For each individual in the study, we define the first-stage model by pairing the process model in Johnson
et al. (2008) and the measurement error model described in Section 3.4 that matched the individual’s device type. We then used the `R`

package `crawl`

(Johnson, 2016) to fit the hierarchical model to the individual’s telemetry observations, and draw realizations from the posterior distribution of the continuous movement process. Each draw was evaluated on a daily time scale, corresponding to the rate at which we have observations of sea ice extent.

Priors for the sub-population covariances and

defined in Section 3.3 were defined through the eigenvalues and angle of rotation of the ellipses. We specified independent priors from a Gamma distribution with mean equal to 150 and variance equal to 10,000 for each eigenvalue, and a uniform distribution on

for the angle of rotation. This corresponds to our belief that the sub-population effect should have diminishing influence beyond distances of approximately 150km from the activity centers, and that these sub-population distributions may have arbitrary orientation.Mean parameters for the prior distributions on measurement error variance processes, , corresponding to cases in which an Argos device was used but no auxilliary ellipse information is available were specified as for Arogs error class respectively (Costa et al., 2010) (see Section 3.5).

### Stage 2

We drew realizations from the posterior distribution of the parameters using Metropolis within Gibbs updates in an MCMC procedure. We ran the MCMC algorithm for a total of iterations, discarding the first half as burnin.

## Appendix B: Summarizing the spatial boundary

Summarizing the distribution of the derived spatial boundary requires defining a meaningful center and spread of a 6-dimensional random variable. Each spatial boundary is a conic section corresponding to the points satisfying for some uniquely defined coefficients . We define a central measure of the distribution of as the conic section defined by the coefficient-wise means , where each mean is computed for a sample from the posterior.

How to define a useful summary of the variability around this central curve is not obvious. We elect to define variability through a point-wise, equi-tailed 95% confidence band orthogonal to the central measure defined above. This summary may be regarded as somewhat conservative, because the probability of drawing a realization from the posterior distribution that lies entirely withing the confidence band is at least 0.95.

## Acknowledgments

The authors thank Franny Buderman, Perry Williams, and Christopher Peck for insight that improved the article. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

The views of the USFWS authors in this publication are solely those of the USFWS authors, and do not necessarily represent the views of the USFWS. This article has been peer-reviewed and approved by USGS under their Fundamental Science Practices policy (http://pubs.usgs.gov/circ/1367). This research was permitted under the Marine Mammal Protection Act and Endangered Species Act under U.S. Fish and Wildlife Service permits (MA046081, MA 690038) and followed protocols approved by Animal Care and Use Committees of the U.S. Fish and Wildlife Service and U.S. Geological Survey.

The authors acknowledge support for this research from NSF DMS 1614392 and USFWS G17AC00068.

## References

- Aikens et al. (2017) Aikens, E. O., Kauffman, M. J., Merkle, J. A., Dwinnell, S. P., Fralick, G. L., and Monteith, K. L. (2017). The greenscape shapes surfing of resource waves in a large migratory herbivore. Ecology Letters 20, 741–750.
- Amstrup et al. (2005) Amstrup, S. C., Durner, G. M., Stirling, I., and McDonald, T. L. (2005). Allocating harvests among polar bear stocks in the Beaufort Sea. Arctic 53, 247–259.
- Atwood et al. (2016) Atwood, T. C., Peacock, E., McKinney, M. A., Lillie, K., Wilson, R. R., Douglas, D. C., Miller, S., and Terletzky, P. (2016). Rapid environmental change drives increased land use by an Arctic marine predator. PLoS ONE 11, e0155932.
- Boyce and Mcdonald (1999) Boyce, M. S. and Mcdonald, L. L. (1999). Relating populations to habitats using resource selection functions. Trends in Ecology & Evolution 14, 268–272.
- Brost et al. (2015) Brost, B. M., Hooten, M. B., Hanks, E. M., and Small, R. J. (2015). Animal movement constraints improve resource selection inference in the presence of telemetry error. Ecology 96, 2590–2597.
- Brost et al. (2016) Brost, B. M., Hooten, M. B., and Small, R. J. (2016). Leveraging constraints and biotelemetry data to pinpoint repetitively used spatial features. Ecology 98, 12–20.
- Buderman et al. (2016) Buderman, F. E., Hooten, M. B., Ivan, J. S., and Shenk, T. M. (2016). A functional model for characterizing long-distance movement behaviour. Methods in Ecology and Evolution 7, 264–273.
- Christ et al. (2008) Christ, A., Hoef, J. V., and Zimmerman, D. L. (2008). An animal movement model incorporating home range and habitat selection. Environmental and Ecological Statistics 15, 27–38.
- Costa et al. (2010) Costa, D. P., Robinson, P. W., Arnould, J. P. Y., Harrison, A. L., Simmons, S. E., Hassrick, J. L., Hoskins, A. J., Kirkman, S. P., Oosthuizen, H., Villegas-Amtmann, S., and Crocker, D. E. (2010). Accuracy of ARGOS locations of pinnipeds at-sea estimated using fastloc GPS. PLoS ONE 5, e8677.
- Durner et al. (2017) Durner, G. M., Douglas, D. C., Albeke, S. E., Whiteman, J. P., Amstrup, S. C., Richardson, E., Wilson, R. R., and Ben-David, M. (2017). Increased Arctic sea ice drift alters adult female polar bear movements and energetics. Global Change Biology 23, 3460–3473.
- Durner et al. (2009) Durner, G. M., Douglas, D. C., Nielson, R. M., Amstrup, S. C., McDonald, T. L., Stirling, I., Mauritzen, M., Born, E. W., Wiig, Ø., DeWeaver, E., Serreze, M. C., Belikov, S. E., Holland, M. M., Maslanik, J., Aars, J., Bailey, D. A., and Derocher, A. E. (2009). Predicting 21st-century polar bear habitat distribution from global climate models. Ecological Monographs 79, 25–58.
- Fetterer et al. (2010) Fetterer, F., Savoie, M., Helfrich, S., and Clemente-Colón, P. (2010). Multisensor analyzed aea ice extent - Northern hemisphere (MASIE-NH), version 1, updated daily. National Ice Center and National Snow and Ice Data Center 2008-2016, Boulder, Colorado USA. accessed August 13, 2017.
- Hanks et al. (2011) Hanks, E. M., Hooten, M. B., Johnson, D. S., and Sterling, J. T. (2011). Velocity-based movement modeling for individual and population level inference. PLoS ONE 6, e22795.
- Hjermann (2000) Hjermann, D. O. (2000). Analyzing habitat selection in animals without well defined home range. Ecology 81, 1462–1468.
- Hooten et al. (2010) Hooten, M. B., Johnson, D. S., Hanks, E. M., and Lowry, J. H. (2010). Agent-based inference for animal movement and selection. Journal of Agricultural, Biological, and Environmental Statistics 15, 523–538.
- Hooten et al. (2017) Hooten, M. B., Johnson, D. S., McClintock, B. T., and Morales, J. M. (2017). Animal Movement: Statistical Models for Telemetry Data. Chapman & Hall/CRC, Boca Raton, Florida, USA.
- Johnson et al. (2006) Johnson, C. J., Nielsen, S. E., Merrill, E. H., McDonald, T. L., and Boyce, M. S. (2006). Resource selection functions based on use-availability data: theoretical motivation and evaluation methods. Journal of Wildlife Management 70, 347–357.
- Johnson (2016) Johnson, D. S. (2016). crawl: Fit continuous-time correlated random walk models to animal movement data.
- Johnson et al. (2008) Johnson, D. S., London, J. M., Lea, M.-A., and Durban, J. W. (2008). Continuous-time correlated random walk model for animal telemetry data. Ecology 89, 1208–1215.
- Johnson et al. (2008) Johnson, D. S., Thomas, D. L., Ver Hoef, J. M., and Christ, A. (2008). A general framework for the analysis of animal resource selection from telemetry data. Biometrics 64, 968–976.
- Lele and Keim (2006) Lele, S. R. and Keim, J. L. (2006). Weighted distributions and estimation of resource selection probability functions. Ecology 87, 3021–3028.
- Manly et al. (2002) Manly, B. F. J., Mcdonald, L. L., and Thomas, D. L. (2002). Resource Selection by Animals: Statistical Design and Analysis for Field Studies. Springer, second edition.
- McClintock et al. (2015) McClintock, B. T., London, J. M., Cameron, M. F., and Boveng, P. L. (2015). Modelling animal movement using the Argos satellite telemetry location error ellipse. Methods in Ecology and Evolution 6, 266–277.
- Niu et al. (2016) Niu, M., Blackwell, P. G., and Skarin, A. (2016). Modeling interdependent animal movement in continuous time. Biometrics 72, 315–324.
- Northrup et al. (2013) Northrup, J. M., Hooten, M. B., Anderson, C. R. J., and Wittemyer, G. (2013). Practical guidance on characterizing availabilty in resource selection functions under a use-availability design. Ecology 94, 1456–1463.
- Obbard et al. (2010) Obbard, M. E., Thiemann, G. W., Peacock, E., and Debruyn, T. D. (2010). Proceedings of the 15th working meeting of the IUCN/SCC polar bear specialist group, 29 June - 3 July 2009, Copenhagen, Denmark. In IUCN, volume vii, page 235.
- Rode et al. (2014) Rode, K. D., Regehr, E. V., Douglas, D. C., Durner, G. M., Derocher, A. E., Thiemann, G. W., and Budge, S. M. (2014). Variation in the response of an Arctic top predator experiencing habitat loss: Feeding and reproductive ecology of two polar bear populations. Global Change Biology 20, 76–88.
- Rode et al. (2015) Rode, K. D., Wilson, R. R., Regehr, E. V., Martin, M. S., Douglas, D. C., and Olson, J. W. (2015). Increased land use by Chukchi Sea polar bears in relation to changing sea ice conditions. PLoS ONE 10, e0142213.
- Rubin (1996) Rubin, D. B. (1996). Multiple imputation after 18+ years. Journal of the American Statistical Association 91, 473–489.
- Scharf et al. (2016) Scharf, H. R., Hooten, M. B., Fosdick, B. K., Johnson, D. S., London, J. M., and Durban, J. W. (2016). Dynamic social networks based on movement. Annals of Applied Statistics 10, 2182–2202.
- Scharf et al. (2017) Scharf, H. R., Hooten, M. B., and Johnson, D. S. (2017). Imputation approaches for animal movement modeling. Journal of Agricultural, Biological, and Environmental Statistics 22, 335–352.
- Scharf et al. (2018) Scharf, H. R., Hooten, M. B., Johnson, D. S., and Durban, J. W. (2018). Process convolution approaches for modelling interacting trajectories. Environmetrics .
- Service Argos (2016) Service Argos (2016). Argos user’s manual. URL http://www.argos-system.org/manual/.
- Steele et al. (2015) Steele, M., Dickinson, S., Zhang, J., and W. Lindsay, R. (2015). Seasonal ice loss in the Beaufort Sea: Toward synchrony and prediction. Journal of Geophysical Research: Oceans 120, 1118–1132.
- U.S. Fish and Wildlife (2016) U.S. Fish and Wildlife (2016). Polar bear (Ursus maritimus) conservation management plan, final. U.S. Fish and Wildlife, Region 7, Anchorage, Alaska pages 1–104.
- Ware et al. (2017) Ware, J. V., Rode, K. D., Bromaghin, J. F., Douglas, D. C., Wilson, R. R., Regehr, E. V., Amstrup, S. C., Durner, G. M., Pagano, A. M., Olson, J. W., Robbins, C. T., and Jansen, H. T. (2017). Habitat degradation affects the summer activity of polar bears. Oecologia 184, 87–99.
- Warton and Shepherd (2010) Warton, D. I. and Shepherd, L. C. (2010). Poisson point process models solve the ”pseudo-absence problem” for presence-only data in ecology. Annals of Applied Statistics 4, 1383–1402.

Comments

There are no comments yet.