Linking northern fur seal dive behavior to environmental variables in the eastern Bering Sea

Northern fur seals (Callorhinus ursinus) breeding on the Pribilof Islands, Alaska have declined dramatically over the past 40 years. Effective conservation of northern fur seals depends on understanding the foraging behavior of adult females whose foraging success is linked to pup survival. We determined the foraging behavior for 11 tagged lactating female northern fur seals from the Pribilof Islands using a statespace modeling approach with an autoregressive movement model. To interpret at-sea behavior in the context of oceanic habitat, we spatially and temporally matched high-resolution reconstructed tracks to a set of environmental covariates that included: commercial groundfish catch, sea surface temperature, primary productivity, wind speed, depth and time of day. We used a Bayesian hierarchical framework to implement a multinomial regression model to link behavior to environmental covariates and account for the mismatch of scale between fur seal behavior and the environmental variables by incorporating an errorin-covariates approach into the hierarchical model. The Bayesian framework allowed us to build a single model to synthesize the information from all the northern fur seal foraging tracks and the available information about the underlying environmental conditions. Application of the approach indicated that the behavioral states for the northern fur seal were significantly related to the Alaska commercial groundfish catch, particularly walleye pollock (Gadus chalogramma).


INTRODUCTION
The population of northern fur seals (NFS) in the Pribilof Islands of Alaska has declined dramatically during the past 40 years, and continues to decline without any obvious reason yet identified (Towell et al. 2006, Lee et al. 2014, Towell et al. 2014).Conservation efforts require a foundation of scientific understanding about NFS ecology, a key element of which is foraging strategy.The foraging strategy of adult female NFS is of particular interest as success in finding enough prey to eat can be linked not only to adult survival but also to pup survival (Antonelis et al. 1997).For a lactating female, foraging success is a time-limited trial in which enough energy must be obtained from prey caught in a changing ocean environment to maintain a dependent and growing pup left behind at her natal rookery.Failure implies trouble meeting the energy needs of her pup and/or herself, and could affect their survival.The relationship between lactating NFS movement and habitat is shaped by foraging success and the physiological constraints of feeding a stationary, land-based pup.Understanding the relationship between movement, behavior, prey species, and the oceanic environment is a critical part of understanding population processes (Bowler andBenton 2005, Patterson et al. 2009).Indeed, the National Marine Fisheries Service's NFS conservation plan has identified this as a highest-level priority (NMFS 2007).
We examined a set of lactating northern fur seals that were tagged on the Pribilof Islands in the Bering Sea during the breeding season of 2005 and 2006 in order to link foraging behavior to environmental conditions and prey fields.Marine animal movement is a complex 3-D process that does not always simplify into lower dimensions, and there is mounting evidence to demonstrate the potential perils of inferring animal behavior based on horizontal trajectories alone (e.g., McClintock et al. 2013 and references therein).For instance, horizontal straightness indices (i.e., area restricted searches) poorly correlate with feeding behaviors (Austin et al. 2006, Weimerskirch et al. 2007).First passage time is a movement metric that measures animal passage through a horizontal region of fixed radius and has been linked to environmental variables, but it can be confounded by the slow speeds of resting behaviors, and the fast speeds of foraging behavior along tortuous paths; both behaviors lead to similar times to traverse a fixed radius region.Depth information has been widely used to determine dive type, but is open to subjectivity in determining what dive types or dive metrics to link to foraging behavior (e.g., Gentry et al. 1986, Goebel 2002).More recently, classification methods that use state-space mod-els have been used to infer behavioral states (and switching) using horizontal location information (e.g., Jonsen et al. 2007).Extensions of standard state-space methods have been developed to link behavior to environmental covariate information (Dragon et al. 2012, Breed et al. 2012), but it remains a difficult problem to capture the vertical dimension of animal movement and its link to the environment.In this study, we used both horizontal and vertical movement information: high-resolution tracks were constructed from tag orientation and speed data, and depth information was used to classify fur seal behavior using the approach proposed by Dowd and Joy (2011).
The dominant prey in fur seal diets is juvenile walleye pollock (Gadus chalogramma; Perez and Bigg 1986, Sinclair et al. 1994, Zeppelin and Ream 2006, Zeppelin and Orr 2010), a species of groundfish that is commercially caught as adults.However, direct information on prey fields is difficult to collect and is not always available.This is also true for other environmental variables that affect foraging such as bathymetry (Antonelis et al. 1997, Call et al. 2008) and the shelf break (Goebel et al. 1991, Robson et al. 2004, Sterling and Ream 2004), lunar cycle (Ream et al. 2005), thermocline depth, and surface fronts (Nordstrom et al. 2013a, Sterling et al. 2014).The extent to which these variables influence prey fields is poorly understood (Ream et al. 2005, Kuhn et al. 2010).
Quantitatively linking behavior to changes in the environment at the landscape level is important for understanding foraging strategies and fur seal ecology.Bayesian hierarchical models have been suggested as an appropriate statistical framework for doing this (Schick et al. 2008).In this study, we developed and applied such a Bayesian hierarchical model.The primary advantage of the framework was that the uncertainty in both seal behavior and the environmental covariates could be fully accounted for by building a hierarchy of conditional models to describe the complexity of our data and the processes that generated them (Cressie et al. 2009).Our framework emphasizes the uncertainty in modeling behavior through the incorporation of errors-in-covariates, and population inference using individual information.While our application of this method was on the behavior of northern fur seals at a breeding site (rookery) on the Pribilof Islands, it is important to emphasize that our approach for analyzing animal movement in relation to the environment is adaptable for other populations and species.

MATERIALS AND METHODS
Our central goal was to understand how the ocean environment influences the behavior of a population of northern fur seals in the Bering Sea.We used biotelemetry data as well as external sources of information on the oceanic environment.The approach, as detailed below, included determining high-resolution positional information from the tag data, identifying environmental variables along the tracks, performing behavioral inference using a state-space model, and synthesizing all this information using a hierarchical Bayesian statistical model to quanti-fy the relationship between environment and behavior.

Tag data
During the 2005 and 2006 breeding seasons, 18 lactating northern fur seals were captured at Reef Rookery on St Paul Island (57.188 N, 170.278 W;left panel of Fig. 1;5 in 2005, 13 in 2006).Tags were attached mid-dorsally to each fur seal using methods described in Boyd and Croxall (1992).Three types of tags were used: (1) An archival dead-reckoner tag (Driesen & Kern GmbH, Bad Bramstedt, Germany), (2) An ARGOS (Advanced Research and Global Observation Satellite Platform Transmitter Terminal) satellite transmitter tag (Spot5, Wildlife Computers, Redmond, Washington, USA), and (3) A VHF radio transmitter tag (A2920 Glue On, Advanced Telemetry Systems, Isanti, Minnesota, USA).The deadreckoner tags were 10-channel loggers with a There is significant missing catch data off-shelf as there is no groundfish industry here, thus our analysis is biased towards on-shelf behaviors.Identified behaviors of northern fur seal along eleven at-sea foraging tracks of lactating female northern fur seals from St Paul Island, Alaska are also shown.
v www.esajournals.org32-MB archive that recorded time, depth, speed (using a swim paddle system), temperature, light, pitch, roll, compass heading (in threedimensions, 3-D), and body orientation (belly up or belly down).These were programmed to collect data at 2-second or 5-second intervals for the duration of the foraging trip (specifically, two of the 2005 tags were set to sample at 5-second intervals, and the remaining 2005 and all 2006 tags were programmed to sample at 2-second intervals).The satellite tag information (latitude and longitude) was used to calibrate deadreckoner route calculations.The radio tag allowed the animals to be relocated when they returned to the rookery.All fur seals were recaptured in order to remove the tracking devices, and retrieve the time series data logged during their at-sea foraging trips.

Track reconstruction
Information from the dead-reckoner tags and the ARGOS transmitters were combined to produce high-resolution foraging tracks.The dead-reckoner channels recorded compass bearing, inclination, body orientation, depth, and speed-through-the-water for each female.We used the dead-reckoning channel vectors of compass bearing, speed, and depth to reconstruct the 2-and 3-dimensional swim paths of the fur seals (e.g., Wilson et al. 1993, Ropert-Coudert et al. 2002, Shiomi et al. 2008).We processed the tag's speed channel to correct for cumulative errors due to bias in the speed paddle's position.Zero-speed periods were identified as regions with near-zero variance, and linear interpolation was used to re-calibrate the speed record.Each foraging track was then corrected for any speed and direction offsets and ocean drift by constraining the dead-reckoner track to lie between ARGOS satellite locations which have limitations related to accuracy of location and infrequency of satellite locations per day (Loughlin et al. 1999, Patterson et al. 2008).We used only high-quality locations (ARGOS Location Classes 1, 2, 3 i.e., location error , 1.0 km; Mate et al. 1998) to compute speeds to limit the possibility of large location errors.We did not formally accommodate this source of errors in our track reconstructions as 1 km errors were deemed small relative to the finest scale covariate data (1 minute of a degree or 1.9-km resolution).We then translated ARGOS locations into Universal Transverse Mercator (UTM) coordinates, and translated the dead-reckoner track into polar coordinates by rotating the angle of movement, and rescaling the radial coordinate to match the direction and great circle distance between consecutive AR-GOS locations.When applied to all fur seals, this yielded a 2-or 5-second resolution track reconstruction for each that in-filled between the ARGOS satellite locations.In this way, the track was linked in space and time to the environmental conditions encountered while foraging at-sea.All location and movement analyses were conducted in UTM units and back-transformed to geographic coordinate units for presentation purposes.

Fur seal behavior
Fur seal behavior along a foraging track was determined using the state-space modeling approach proposed by Dowd and Joy (2011).The analysis proceeds by differencing the tag's depth channel to create a measure of vertical velocity.A vertical movement model based on a second order auto-regressive model was used where w t represents the vertical velocity, and a 1 and a 2 are the autoregressive parameters.The error process was assumed to be a zero-mean white noise Gaussian process with variance parameter r 2 w .A unit time increment here corresponded to the 2-or 5-second resolution of the tag data.The data along each fur seal track is sectioned into 26-minute windows (with 13-minute overlap between adjacent windows) allowing us to identify approximately stationary ensembles of dives from which to determine behavior.For each time window, the parameters a 1 , a 2 and r 2 w are estimated, taking into account observation error.This was accomplished using an augmented state-space model and multiple iterated filtering (see Dowd and Joy 2011 for full details of the methodology).The along-track parameter estimates were then further smoothed using a locally optimized kernel smoother (Herrmann 2013, R library ''lokern'') to identify longer period behavior, and compensate for short-term random parameter fluctuations.
The estimated along-track movement parameters derived from the tag data were then classified into behavior types.We considered three discrete behavior classes: (1) non-diving; (2) active diving; and (3) exploratory diving.The non-diving behavior corresponded to near-zero values for the estimated process error variance r 2 w , and includes sleeping and resting states, as well as surface transiting; these behaviors are characterized by a lack of engagement in the immediate environment.The behavioral states corresponding to active and exploratory diving are diagnosed from the estimated values for the auto-regressive parameters a 1 and a 2 .Time series theory allows one to define the dynamic system behavior of an AR(2) model based on the values of a 1 and a 2 (Priestley 2004, Shumway andStoffer 2006).We defined the a 1 , a 2 values that corresponded to periodic behavior to represent active diving (i.e., regular and repeated dives).The other a 1 and a 2 values that corresponded to nonperiodic behavior were identified as exploratory diving (i.e., less regular and intermittent dives).
See Dowd and Joy (2011) for further details.
Fig. 2 provides a concrete illustration of this behavioral classification using the movement parameters.It shows the vertical velocity data, the estimated values of a 1 and a 2 , and the corresponding behavioral classifications for a single day (August 18, 2006) for one fur seal track showing a clear correspondence between the estimated parameter values and the behavioral type.Thus we have taken the output of a state space model that related movement characteristics in vertical speed to coherent bouts or segments of behavior.In the analysis that follows, we concatenate adjacent 26-minute windows of similar behavior into a single observation of the classified behavior that then acts as the response variable for that segment of track.For each fur seal, the series of behaviors observed along the length of track is then associated with a set of space and time positions in the Bering Sea that reflect the habitat where Fig. 2. Time series for a single day, August 18, 2006 for a single track.Shown are the estimated AR(2) model parameters a 1 and a 2 ; kernel smoothed (blue line) and original estimates (black dots).These are overlaid with the time series of vertical velocity (black lines).The yellow, pink and grey blocks correspond to active diving, exploratory diving, and non-diving, respectively (as diagnosed from values of a 1 and a 2 ).v www.esajournals.orgthe fur seal engaged in each behavior.This will be explicitly linked to the spatio-temporal environmental covariate information in the following section (At-sea environmental conditions), using the Bayesian hierarchical model (Hierarchical Bayesian model below).

At-sea environmental conditions
We considered a variety of available environmental data from the Bering Sea that could plausibly be linked to northern fur seal foraging behavior.We included information on fisheries (as a proxy for prey fields), as well as physical and biological oceanographic characteristics that are detailed below.
Commercial groundfish catch and walleye pollock.-Amajor goal of the study was to ascertain the extent (if any) to which fur seal behavior was associated with fish abundance and, in particular, walleye pollock density.As a proxy for fish abundance, we used the US Department of Commerce domestic observer data of the Alaska groundfish fishery for 2004-2007(NMFS 2012)).We limited the fish catch data to be between July 9, the beginning date of the breeding season (Trites 1992), and November 11, the median dispersal dates for pups on St Paul Island in 2005 (Lea et al. 2009).We spatially linked the fur seal tracks to individual groundfish haul at the nearest minute of longitude and latitude ( 1.9 km).We selected two variables of interest: (1) the haul weight of walleye pollock (85% of the total catch weight was walleye pollock), and (2) the total catch weight (including both retained and discarded species).Where multiple hauls were linked to a single track segment, we took the median pollock haul weight and total catch weight to represent that segment of track (right panel of Fig. 1).If no catch was linked to a segment of track (i.e., off-shelf regions), this was considered missing covariate data.Missing catch data was not in-filled with zeros so that the catch data could be used as a proxy for abundance of prey, rather than directly representing fisheries catch (and confounding issues of 0 catch with 0 effort).
Sea surface temperature.-Seasurface temperature (SST) is a primary oceanographic variable that is easily measured (Nordstrom et al. 2013b), and may influence the behavior of northern fur seals (Nordstrom et al. 2013a).Here, temperature was recorded directly by the dead-reckoner tag at 2-or 5-second intervals.To isolate SST, we removed all the data that were observed outside of the top two meters of the ocean surface.The median sea surface temperature along each behavioral segment for each fur seal was extracted and assigned as an environmental covariate.
Primary productivity.-Oceanprimary productivity was obtained from the NOAA CoastWatch net primary productivity for the Gulf of Alaska and Bering Sea.This is based on satellitecollected chlorophyll-a concentration and photosynthetically available radiation (PAR) measurements, corrected for the amount of organic carbon used by planktonic organisms in respiration (Behrenfeld and Falkowski 1997).We used gridded spatial maps at 1/6th degree resolution (;18.5 km) processed as 8-day time-averages over the domain and period of interest.Primary productivity was assigned to the foraging track by linearly interpolating these data so as to match the time and location of the fur seal track.
Wind speed.
-To account for the potential effect of sea-state on foraging behavior, wind speed was determined along the fur seal tracks.We used the National Climatic Data Center wind product, which blends satellite wind speeds from multiple platforms such as scatterometers, and passive microwave radiometers (Zhang et al. 2006).For each 24-hour period, we extracted daily wind speeds using the highest resolution 1/ 4 degree (;28 km) gridded wind fields.These were linearly interpolated to match the location of each behavioral segment at the relevant time.
Ocean depth.-Thebathymetry of this region of the Bering Sea where the foraging tracks are located is characterized by shallow on-shelf waters ,200 m deep, and off-shelf waters with ocean depths of 3000 m or more.Note the maximum dive depth of a northern fur seal is ;200 m (Gentry et al. 1986).Hence, we simplified bathymetry into an on-shelf/off-shelf categorical variable (We also tried the definitions of Call et al. (2008) that defined a three-level categorical depth variable: inner-middle shelf (0-100 m), outer shelf (100-200 m), and off-shelf (.200 m), but the additional resolution was not informative in our study).Under this definition, ''off-shelf'' waters include the deeper regions at the shelf break associated with foraging behavior in northern fur seals (Goebel et al. 1991, Robson et al. 2004, Sterling and Ream 2004).We matched longitude and latitude of the track to the nearest depth measurement from the ETOPO1 bathymetry model (Amante and Eakins 2009), which gives bathymetry with resolution 1/60th of a degree, (;1.9 km).We then calculated the median depth for each of the behavior segments along the at-sea track.For example, if the seal was in an active diving mode between midnight and 4am, we took the median depth of the ocean during that period of the track, and classified it as on-shelf or off-shelf.
Time of day.-Time for each behavioral state was obtained from the dead-reckoner tag at the start of each behavioral segment.As northern fur seals have circadian patterns in behavior (Ream et al. 2005), we transformed the local time of day (t) into a circular variable of the form Time ¼ AcosðxtÞ þ BsinðxtÞ where x ¼ 2p/24 hr.This ensured that the interpretation of time at 24:00 would be equivalent to that at 00:00.This provides a sinusoidal covariate function representing time of day, wherein the coefficients A and B are the focus of inference.

Hierarchical Bayesian model
We followed a hierarchical Bayesian modeling approach to examine possible relationships between behavior and the environmental covariates that describe fur seal at-sea habitat.This statistical methodology addresses the linking of behavior to habitat by decomposing the overall problem into a set of simpler and more tractable problems.Specifically, the complex joint probability distributions is expressed as a series of conditional models that aim to explain the relationship between data and a causal process as follows (see Berliner 1996, Cressie et al. 2009

Observation model
Our observation model linked the set of observed outcomes, y ij , to their outcome probabilities, p ij , using a multinomial probability model.Specifically, the observation vector, y ij ¼ ðy ð1Þ ij ; y ð2Þ ij ; . . .; y ðKÞ ij Þ, characterized the behavioral state of the ith fur seal over the jth segment of track.Each y ðkÞ ij for k ¼ 1 to K was equal to 1 if the animal was observed in the kth behavioral state, and 0 otherwise.Correspondingly, the probability of observing the ith fur seal at the jth location engaged in the kth behavioral mode was written in vector form as

Local process model
The process model linked the probability of being in a particular behavioral state, or the latent probability vector p ij , to the environmental covariates.This was done using a multinomial regression model (McCullagh and Nelder 1989) according to the following procedure.Let X ij ¼ ðX ij1 ; X ij2 ; . . .; X ijp Þ T be the vector containing the true values for each of the p covariates that characterize the environment of the ith fur seal over the jth segment of track ( T designates the vector transpose).We modeled the logarithm of the ratio of the probability of each category relative to that of a baseline category; this is done since the total probability is one, and hence one of the behavioral categories is redundant.We selected the most commonly observed behavior (here, non-diving) as our baseline category (k ¼ 1), following Agresti (1990).The multinomial logit linear regression model is then log Pðy The process and observation models were combined by defining the joint likelihood of the three behavioral outcomes using the observation model and the multinomial regression, for all K track segments for the ith fur seal, can then be expressed as: where J i is the total number of behavioral segments for the ith fur seal.
One key aspect of this study was the need to account for errors in the covariate data (Stephens andDellaportas 1992, Carroll et al. 2006).The usual assumption in regression models is that all covariates have been measured without error, but violations of this assumption generate biased and inaccurate inferences (Gustafson 2003).In our case, these arose due to intrinsic measurement errors in the environmental variables, as well as the use of behavioral segments as the analysis units.These lead to errors of representativeness introduced through spatial and temporal averaging (i.e., aggregated data; see Gotway and Young 2002).
Here, we explicitly accounted for errors in the environmental covariates.To do so, we designate the observed value of the covariates along the behavioral segment as W ij , where W ij is related to its true value, X ij , according to where U ij describes the observation error associated with the ith fur seal on the jth track segment and we assumed that U ij ; N ð0; r 2 iU Þ so that the variance r 2 iU is specific to an individual fur seal.Note that this is a random effect model wherein the true values are determined through an errorprone but unbiased measure.
Errors in environmental covariates can be incorporated into the hierarchical model as a prior distribution for X i .We selected the multivariate normal distribution (MVN) as the form of that prior; i.e., PðX i Þ ; MVN ðW i ; r 2 iU I J i Þ where I J i denotes an identity matrix of dimension J i .The full Bayesian specification also requires us to specify a prior on its variance r 2 iU .Specifically, we select an independent conjugate inverse gamma distribution or r 2 iU ; IGða U ; b U Þ. The parameters we chose for this distribution were a U ¼ 3, and b U ¼ 1, and limited the sampling range for the prior to be within a factor of 5 standard deviations of the prior distribution mean.In practice, limiting the sampling range for the prior only affected a small number of cases for which there was limited information in the covariate data to estimate r 2 iU .A detailed description of all Bayesian priors used in our hierarchical model appear in the Appendix.
The following environmental covariates were modeled with error: total catch weight of groundfish; haul weight of walleye pollock; primary productivity; sea surface temperature; and wind.Total catch weight and haul weight of walleye pollock were positively correlated, r ¼ 0.85.Therefore, to limit the detrimental effects of multicolinearity, these two groundfish measures were never both incorporated into the same model.The following covariates were considered to be measured without error: bathymetry (onshelf/off-shelf ); and time of day.

Parameter model
The parameter model linked the lower-level models for each individual fur seal, and so may be viewed as a population model.The main idea is that by linking the regression models for the individual fur seals via a shared prior, we incorporate population variability (Cressie et al. 2009), and borrow model strength between animals (Ntzoufras 2009).This was done by selecting a joint prior for the lower-level regression parameters, b i .That is, the regression parameters specific to a single animal (e.g., b i ) are thought of as a random sample of coefficients from a distribution of possible values, and the posterior of each b i is a weighted mean of the animal's regression parameters and the overall population effect.
We assumed that each of the components of b i was normally distributed, and the joint prior distribution of b i was then a multivariate normal likelihood, i.e., Again, the fully Bayesian specification requires priors to be specified for the higher-level parameters B and R b of the MVN likelihood of Eq. 1.1.The following inverse Wishart (IW) and MVN distributions were selected as priors for these parameters, Note that the priors for both B and R b were conjugate to the higher-level MVN likelihood in Eq. 1.1.Full derivations of conditional distributions appear in the Appendix.

Bayesian computation and statistical assessment
We have presented a hierarchical multinomial Bayesian error-in-covariate model to describe the relationship between a fur seal's behavior and her environment.To solve for the target posterior distribution, we implemented a Markov Chain Monte Carlo (MCMC) Metropolis-within-Gibbs sampling algorithm that was coded using R statistical software (R Core Development Team 2012).The MCMC was run for 1,000,000 iterations (discarding 10,000 iterations for burn-in).The partial correlation coefficient was used to determine the decimation rate for the MCMC chain (i.e., we fit autoregressive models of successively higher-orders until the lag suggested non-significant partial correlations between chain components).As a consequence, chains were thinned to every 50th iteration.To check robustness, we initialized the chains from three different starting points using different random seeds.Convergence was assessed by inspection of trace plots, and using the Gelman-Rubin statistic ( R) as modified by Brooks and Gelman (1998; using R library ''coda'') as a quantitative measure of convergence.

Model comparison and model adequacy
Model comparison and selection was undertaken using information-theoretic approaches.In our study, we were interested in understanding fur seal behavior in general and not simply the behavior of the sampled fur seals.Since inference around population parameters was our focus, we compared the likelihood of each model's population parameters.Towards this end, we used the following metrics: (1) Akaike Information Criterion to compare model fit (AIC, Akaike 1973); (2) the posterior mean deviance as a Bayesian measure of fit (DðbÞ, Dempster 1974, Spiegelhalter et al. 2002); and (3) m DIC as a diagnostic of leverage (Spiegelhalter et al. 2002).Note that these quantities provide useful comparisons between candidate models, but do not provide insight into model adequacy.
We examined model adequacy through distributional summaries.We compared posterior predictive distributions of replicated data (produced under the model assumptions) with the observed data (Rubin 1984, Gelman et al. 2004).Assessment of the overall goodness of fit used posterior predictive P values ( p pp , Meng 1994, Carlin andLouis 2000).Due to uncertain large sample properties of posterior predictive P values (see discussion in Bayarri and Berger 2000), we assumed that a P value around 0.5 indicated that the distributions of the replicated and actual data were close, while a value close to zero or one indicated strong differences between them (Gelman et al. 2004).

Sample size effects
We ran simulations to examine the effect of sample size (i.e., the number of tagged fur seals) on model inference.This was done by drawing sets of simulated fur seals with various sample sizes from the model's (population) posterior.The analysis was then repeated for synthetic data sets ranging from 5 to 100 fur seals.We then compared these posteriors to the original posterior via the Kullback-Leibler (K-L) divergence statistic (Kullback and Leibler 1951).This allowed us to quantify the discrepancy in posteriors that was a consequence of sample size.The K-L divergences were calculated using the R library ''flexmix'' (Leisch 2004).

Foraging tracks
Of the 18 lactating female northern fur seals captured at Reef Rookery on St Paul Island during the pupping seasons of 2005 and 2006, 10 individuals yielded complete data sets from which the foraging paths and behavior information could be retrieved.One of these animals provided data from two complete foraging trips.The incomplete data records of the other animals were due to tag failures including: water entering the housing, the speed paddle breaking off, fish vertebrae lodging between the speed paddle and housing, and battery failure.Hence, our data set for the analysis comprised of 11 complete northern fur seal foraging tracks.Foraging trips ranged in length from 5.5 days to 11.2 days, with a median of 7 days spent at sea.An average of 40 ARGOS positional fixes and 404,400 archival data records were obtained per trip.
The right panel of Fig. 1 maps the track reconstructions of 11 fur seal at-sea trips (with colors representing behavioral classifications).The fur seals traveled an average maximum linear distance from St Paul Island of 279 km, with a maximum distance away of 391 km.In general, these fur seals dispersed from the rookery in all directions, except directly to the north.They covered a wide area of both on-shelf and off-shelf waters in the Bering Sea and showed no preference for foraging in any single area.Five animals stayed on the continental shelf, one animal moved along the shelf break foraging in the canyons, and five animals went across the shelf break into the deep water of the central Bering Sea.The one animal for which we had two successive trips swam to a similar off-shelf location in both her at-sea trips.

Behavioral inference
The track reconstructions of Fig. 1 show the results from the behavioral classification along the foraging tracks.The track segments corresponding to the three behavioral classes (active diving, exploratory diving, and non-diving) identified according to the state-space algorithm of Dowd and Joy (2011).Far from the rookery, a mix of behaviors with a preponderance of active diving can be seen for all tracks.There was relatively more exploratory diving nearer the rookery on the outgoing legs, while the nearrookery incoming legs are characterized by nondiving transiting behavior.Overall, the fur seals spent 35% of the time at sea engaged in either active (14%) or exploratory dives (21%).Fur seals spent 65% of time at the surface engaged in nondiving behaviors.The 11 foraging trips had a median of 41 segments of different behavior modes, with those animals taking longer foraging trips exhibiting a greater number of behavior changes.

Relating behavior to environment
We used the hierarchical Bayesian model (see Materials and methods: Hierarchical Bayesian model ) to investigate the link between fur seal behavior and the environment.A forward selection based model building strategy was undertaken using  1.
We first fit main effect models to each of our at-sea environmental covariates.We then fit a series of models that all included time of day (Time) and incorporated the other covariates one at a time (in an additive fashion).The time of day models fit with Pollock, Total Catch and Bathymetry (On/Off Shelf ) were found to be the best subset according to the model selection criteria.There was no relation to Primary productivity.We failed to get satisfactory convergence ( R ) 1) in the models with Sea surface temperature and Wind, and these results are omitted from Table 1.
We next fit a set of models that include Time and one other main effect, and a Time 3 main effect interaction.The two models with Time interacting with a measure of groundfish catch (both Pollock and Total Catch) had the best fit diagnostics.The model that included Time of day and Primary productivity failed the posterior predictive test, suggesting that the model did not fit the data well, despite having adequate higherorder parameter convergence in their respective MCMC's (i.e., R ' 1).
One additional step of model complexity was also considered.We fit a three-covariate model based on the best time-interaction model, and bathymetry (Time 3 Pollock þ On/Off Shelf; Table 1).This model was penalized aggressively by the AIC statistic, suggesting the additional complexity was not warranted despite a fair goodness of fit to the data.Table 1 reports a subset of the results from the model fitting.
Our final selected model from this procedure included the main effects of time of day and (log) haul weight of walleye pollock and included an interaction term of these covariates: Time 3 log(Pollock) (bold line in Table 1).It is notable that the model that substituted Total Catch for Pollock was a similarly good model with only marginally poorer fit (Time 3 log(Total catch) in Table 1).Coefficient magnitudes, and signs were similar for both models.Detailed summary statistics for posterior distributions of the regression coefficients can be found in Table 2).Fig. 3 graphically depicts magnitude and 95% credible intervals for each of the population parameters, and shows significance for the main effect of time of day and the interaction of time of day and (log) haul weight of walleye pollock.
To illustrate the effect of incorporating errors in covariates, we also fit the final selected model but assumed the covariates had no error (Fig. 3).By including errors in covariates, the regression parameters for those covariates with error were shrunk towards zero, and tended to show a loss of precision (i.e., there is increase in the uncertainty estimates).Also, this figure shows the parameter estimates of exactly measured covariates (i.e., those not modeled with error such as time of day) were significantly biased when the covariate errors in the remaining covariates were , or log odds of active diving vs. baseline non-diving).
à R is the Gelman-Rubin Bayesian measure of convergence; values near 1 implies good convergence of the MCMC chain.
v www.esajournals.orgignored.The implication is that more of the parameters would be significant if error in the covariates had been omitted, and those parameters would be biased.However, the model fit diagnostics for the case of no errors in covariates were found to be significantly worse (without error: AIC ¼ 297.0,DðbÞ ¼ 329.4,p pp ¼ 0.46 vs. with error: AIC ¼ 284.3, DðbÞ ¼ 323.2, p pp ¼ 0.72).
Fig. 4. Probabilities of behavior modes active diving, exploratory diving, and non-diving in response to time of day (shown for a 24-hour period starting at local noon), and increasing commercial catch size of walleye pollock.The left, middle and right panels of the graph show the predicted relative probabilities of behaviors in areas of small-sized walleye pollock hauls (0.5 tonnes), medium-sized hauls (2 tonnes), and large-sized hauls of pollock (10 tonnes).Note that the x-axis depicts time of day where the axis starts and finishes at noon to highlight the maximum amplitude of active foraging (in yellow) at night.This comparison leads to the conclusion that including covariate error in the model framework provides a significantly better model fit, as well as less biased parameter estimates.
The model results for the daily cycle of behavior in response to changes in time of day and commercial (log) haul size of walleye pollock are shown in Fig. 4.This figure highlights the increase in the relative probability of engaging in active diving behaviors in response to increasing size of pollock catch.This was accompanied by a near-zero probability of non-diving behaviors during those times of active diving in high catch areas.One caution in the interpretation of this figure is that our model did not predict the length of time that the northern fur seals engaged in particular behaviors, but rather predicted the probability of a behavior beginning at various times of the day.This figure suggests that the peak in the start of exploratory dive behaviors is around noon, and that fur seals are most likely to begin active diving behaviors at just past 8:00 in the evening, and most likely to end diving and begin non-diving behaviors around 5:00 in the morning.Comparing results for a haul size of 10 tonnes vs. 2 tonnes suggests there is a higher probability of starting active diving behaviors sooner in the evening, and hence a lower (in fact, near-zero) probability of non-dive behaviors at that time.
Central to the Bayesian paradigm is the notion that as the data quantity (and quality) increase; the posterior becomes less sensitive to prior assumptions (Cressie et al. 2009).Fig. 5 shows a result from the sample size assessment using the Kullback-Liebler divergence metric, K-L.This indicates that the gains in model performance stabilize as the number of tagged fur seals exceeds 20, after which the K-L divergence becomes a relatively flat function.The figure therefore suggests that significant improvement in the inferences could have been gained if about twice as many complete records were obtained for analysis.

DISCUSSION
This study examined the at-sea foraging behavior of a population of lactating northern fur seals on St Paul Island, Alaska.Information on fur seal behavior was obtained from electronic tags for 11 complete foraging trips.Detailed 3-D tracks were reconstructed using high-resolution depth, speed and compass directional information from attached tags, calibrated to positional ARGOS fixes.Fur seal behavior was inferred with a state-space modeling approach using vertical velocity data and movement model; this was the first comprehensive application of the approach proposed by Dowd and Joy (2011) and allowed for successful identification of coherent behavioral segments (non-diving, active and exploratory diving) along the foraging tracks.A comprehensive hierarchical Bayesian framework was then developed to relate dive behavior to the marine environment in which they forage.
Our data were collected during the summer pupping seasons when female northern fur seals are tied to the rookery.Such restricted foraging constrains their ability to explore the foraging habitat in terms of both time spent away, and distance travelled.That is, they must select foraging habitat to obtain energy to feed themselves and their pup, while still returning to the rookery in time to feed a fasting pup.Female northern fur seals in our study embarked on relatively long pelagic journeys averaging 279 km in linear distance from the rookery, and lasting from 5.5 to 11.2 days.This distance is consistent with, but slightly longer than, the  Goebel et al. (1991) (200 km).The length of time of the foraging trip is consistent with reported lengths from other studies (e.g., Gentry andHolt 1986, Nordstrom et al. 2013a; from 1 to 14 days).The tagged fur seals did not appear collectively to prefer any one region around St Paul Island, other than avoiding the continental shelf immediately to the north.Individuals may, however, have preferences: the one female for which we had two consecutive tracks went to a similarly located feeding area off-shelf in both foraging trips.Call et al. (2008) showed that 27 of 36 female fur seals with repeated trips followed the same general direction they used in their previous trips.Other studies have suggested that northern fur seals from Reef Rookery use all hydrographic domains around St Paul Island (e.g., Robson et al. 2004), including both on-shelf and off-shelf habitats (Loughlin et al. 1987, Goebel et al. 1991, Sterling and Ream 2004, Call et al. 2008).Scat samples from the same rookery contained both on-shelf (e.g., walleye pollock, Pacific herring) and offshelf species (e.g., squid, myctophids; Zeppelin and Orr 2010).Our behavior maps showed constant movement and consistently changing diving behaviors throughout the at-sea foraging trips.Active and exploratory dive behaviors were found over much of the foraging paths, and non-diving transiting was most prevalent near the rookery, particularly on the return leg.
Our study suggests that the foraging behavior of lactating northern fur seals may be influenced by the abundance of walleye pollock.Walleye pollock of year-class 0 and 1 are the most common prey in scat studies of the northern fur seal diet (Perez and Bigg 1986, Sinclair et al. 1994, Sinclair et al. 1996, Zeppelin and Ream 2006), and walleye pollock comprised 89.3% of fur seal diet as measured in scats collected within a month of our 2006 tagged females foraging trips from Reef Rookery (Zeppelin and Orr 2010).We found that, in locations with more abundant walleye pollock, northern fur seals were more likely to be actively diving, especially at night.Our proxy for walleye pollock abundance was the US Department of Commerce domestic observer data of the Alaska groundfish industry aggregated over both time and space for the region of interest.It is unlikely that this reflects the real-time abundance of prey that the fur seals actually encountered along their foraging route, and is confounded with fishing effort.However, it may still reflect persistent prey distributions that are attractive to foraging northern fur seals.Nordstrom et al. (2013a) also linked foraging of northern fur seals to oceanic surface fronts (eddies and filaments) and Kuhn (2011), to thermocline depth.These features may also be proxies for the spatial distribution of prey species.
While we used commercial groundfish haul data, we had no direct measure of bycatch of juvenile pollock, nor did we have any small scale evidence that commercial boats were fishing at the same time and location as our tagged fur seals.Hence, we cannot draw conclusions about the potential for competition occurring between the commercial fishery and northern fur seals.However, we can say that being in the same area where fishing occurred would not necessarily imply competitive interactions unless it was known that the fur seals were feeding on the same fish targeted by the fishery, and whether the abundance of that target prey was limited.Evidence from scat analysis suggests that juvenile pollock of year-class 0 and 1 (2-20 cm) are preferred by foraging fur seals (Zeppelin and Ream 2006), which are smaller than the .40cm adults taken by the commercial fishery (Ianelli et al. 2007).As adult walleye pollock are known to prey on age-0 juveniles of their own species (Bailey 1989), it is reasonable to suggest that the link between the distribution of commercially harvested adult pollock and northern fur seals may be confounded due to their sharing of a common prey resource.Alternatively, Gudmundson et al. (2006) investigated diet overlap at a breeding rookery with commercially targeted year-classes of walleye pollock using prey remains in spews, and found significant overlaps with the fishery.Although the percentage of the diet that is regurgitated versus the percentage that is passed through the digestive tract is unknown, the spews indicate that northern fur seals do eat adult size-classes of walleye pollock.
Our results highlight the strong diel pattern in the foraging behavior of northern fur seals, and its link to prey abundance.We found that female fur seals showed strong preferences for active diving at night, while preferring non-dive be-haviors such as resting or transiting in the mornings.Afternoons were typically associated with exploratory dive behaviors.Kuhn et al. (2010) proposed that shallow foraging bouts at night may be related to nighttime migrations of juvenile pollock as they follow vertically migrating zooplankton to the surface.Deeper foraging behavior has been associated with fur seals targeting concentrated groups of juvenile pollock at the thermocline (Kuhn 2011, Nordstrom et al. 2013a).We speculate that if northern fur seals were aware of suitable feeding areas (from prior experience or environmental cues), they could reduce search time and energetic costs by employing strategies (e.g., active or exploratory diving) at particular times along foraging paths that were successful on a previous feeding trip.
We used an auto-regressive movement model to link dive information to foraging behavior.This provided a mechanistic way of categorizing dive types, using interpretable movement signatures as diagnosed from the parameter values of our AR(2) vertical movement model.This contrasts with other classification schemes that distinguished dive types based on depth (e.g., Gentry et al. 1986, Goebel 2002), or used the torturous paths associated with prey patchiness (Benoit-Bird et al. 2013a, Benoit-Bird et al. 2013b).Parameters corresponding to periodic solutions of the AR(2) process were termed active diving.These included shallow, repeated dives with short surface-time intervals, as well as deeper, repeated U-shaped dives with comparatively longer surface recovery times.Such dive types have been attributed to foraging behaviors in an array of top marine predators such as gray seals (Halichoerus grypus; Austin et al. 2006), southern elephant seals (Gallon et al. 2013), Australian fur seals (Arctocephalus pussillus; Arnould and Hindell 2001), harbor seals (Phoca vitulina; Baechler et al. 2002), and others.Parameters corresponding to aperiodic solutions of the AR(2) process, or a correlated vertical random walk, were termed exploratory diving.That is, the typical surface time between dives is longer and more irregular relative to the dive time, and the vertical profile of the dive may be more V-shaped.It is less clear what the underlying motivation for intermittent V-shaped dives is, but others have attributed them to foraging on larger pelagic fish or squid in northern gannets (Morus bassanus; Garthe et al. 2000), or non-foraging activities, including predator avoidance, and explorations in crabeater seals (Lobodon carcinophagus; Bengston and Stewart 1992).The non-diving and non-foraging behaviors such as resting, sleeping, grooming and surface transiting were described, in the AR(2) solution space, as corresponding to periods for which vertical movement followed a white-noise process.
A general hierarchical data analysis and modeling framework is proposed here for interpreting marine animal tag data.It starts with track identification and behavioral inference using a state-space model, and then uses these results within a Bayesian hierarchal model to link fur seal behavior to the ocean environment.The analysis is spatially implicit (not explicit).That is, we did not attempt to answer the question of where in the ocean did the fur seals chose to forage or not forage.With only 11 tracks over an area of .125,000km 2 , this would have been too ambitious.Instead, we addressed the simpler question: given that the fur seal swam through a particular location, could we predict her most likely behavior based on a set of local environmental variables?
We note also that it is possible to consider the entire approach within the Bayesian hierarchical framework (e.g., Bestley et al. 2013).However, for our case it made sense to separate the components, doing the behavioral inference first and then the Bayesian hierarchical model.The main reason was the scale mismatch between the high-resolution movement information (seconds), and that of the behavioral inference and linkage to covariates (hoursþ).This implies that incorporating the state-space-based behavioral inference as part of the hierarchy of the Bayesian model, while conceptually appealing, would have been computationally infeasible due to the MCMC approaches employed.However, by taking the output from the state-space model as input to the hierarchical Bayesian model, there is some risk of misclassification of behavior.
Although it might be possible to write a model that has misclassification error in the response as well as measurement errors in the covariates, these models would have non-separable model errors.The effect of misclassification of a binary response in a logit regression model is to introduce a misclassification parameter that must be additionally estimated (Neuhaus 1999) as well as to reduce significance of the naive estimates of regression parameters (Ku ¨chenhoff et al. 2006).Therefore in this model, the potential error in the classification of behavior would have the same effect as modeling the covariate error by shrinking regression parameters of model parameters towards zero.
Our Bayesian hierarchical model has some important features that should be incorporated into any study attempting to link animal behavior to environmental conditions.First, it focuses on population level inference using a collection of individuals.This was done by imposing a level in the hierarchy that accounted for individual variation.Second, errors in the environmental covariates were modeled explicitly.This is extremely important since not only are the variables themselves subject to measurement errors but they often have a complex error structure since they are frequently themselves data products derived from multiple sources.As well, there is also a great deal of uncertainty associated with spatially and temporally projecting them to match the fur seal tracks.In our study, this is the case for the groundfish catch data.Finally, we remark that since all our MCMC algorithms were coded directly, we had full control of the statistical model building process.This allowed us to tailor the Bayesian computational approaches to the specific problem at hand, and so account for the fact that these sampling-based algorithms are sensitive to assumptions made, and hence must be optimized for reliable posterior inference.
In summary, our data and analysis framework allowed us to link northern fur seal behavior to their ocean environment.The analysis could be further expanded.The sample size analysis suggests that having another 10 tracks would have generated enhanced precision.There are other potentially useful channels on the deployed tags that we have not made use of, particularly the high-resolution 3-D body orientation information that could provide further insight to behavior.We anticipate that other oceanographic information on, for example, frontal eddies and filament locations, could provide additional valuable insight into the association between fur seals and the Bering Sea environment.The general issues of scale, and temporal and spatial dependence, remain outstanding issues in all such problems linking fine-scale movement and behavior to large-scale environmental variables.It can, however, be treated with an errors-incovariate approach within a Bayesian hierarchical framework like the one used here.Finally, we feel that examining movement and foraging behavior of the maternal fur seal in relation to her pup weight gains, rookery residency times, and other life-history and bioenergetic characteristics would further contribute to the understanding of fur seal ecology and pup survival.

Fig. 1 .
Fig. 1.The left panel shows Reef Rookery, St Paul Island, Pribilof Islands, Alaska.Reef Rookery is located on a peninsula that extends off the southern coast of St Paul Island.The intensity of blue represents bathymetry of the surrounding ocean.Adjacent waters are considered ''on-shelf'' and are typically ,200 m.The right panel shows haul weight of walleye pollock calculated from the US Department of Commerce domestic observer data of the Alaska groundfish fishery for 2004-2007 date-restricted to be during northern fur seal summer pupping dates: July 9-November 11, 2004-2007.There is significant missing catch data off-shelf as there is no groundfish industry here, thus our analysis is biased towards on-shelf behaviors.Identified behaviors of northern fur seal along eleven at-sea foraging tracks of lactating female northern fur seals from St Paul Island, Alaska are also shown.
and a 0 ¼ 0.01 (to ensure the covariance matrix is uninformative).Here, I m denotes an identity matrix of dimension m (i.e., R b and I m are m 3 m matrices).

Fig. 3 .
Fig. 3.The 95% credible intervals for 12 regression parameters B from the population level model that includes time of day and (logged) haul size of walleye pollock (Time Pollock) modeled with and without error.Each pair of B coefficients show the credible intervals of the population-level parameters linking northern fur seal behavior to at-sea habitat.B 1:2 denotes the regression parameters corresponding to the logit response log p ð2Þ ij =p ð1Þ ij , or the log odds of exploratory diving vs. baseline non-diving.B 1:3 denotes the regression parameters corresponding to the logit response log p ð3Þ ij =p ð1Þ ij , or log odds of active diving vs. (baseline) nondiving.

Fig. 5 .
Fig. 5. Effect of increasing sample size on Kullback-Leibler (K-L) divergence for one selected populationlevel parameter (cos(Time) 3 log(Pollock): B 13 ).K-L divergence is 0 if and only if the posterior distribution for a simulated dataset is identical to the model from which it was generated.The dashed vertical line at 11 simulated fur seals corresponds to the sample size for this study. ):

Table 1 .
Diagnostics and information criteria for the various Bayesian hierarchical models fit.The final, selected model is noted in boldface.DðbÞ measures model adequacy (smaller implies better).§ m represents the number of population parameters in the upper level model.} m DIC is an estimate of leverage.# p pp is a measure of model predictive power; a value greater than 0.05 implies there is no evidence model is predicting poorly.Detailed results are given in Table

Table 2 .
Posterior summaries for higher level model coefficients, B, in the final, selected model (Time 3 log(Pollock); Table1).Parameters with significant Bayesian P values (* P , 0.05) are noted in boldface.