Model-based clustering reveals patterns in central place use of a marine top predator

Satellite telemetry data are commonly used to quantify habitat selection, examine animal movements, and delineate home ranges. These data also contain valuable information concerning dens, nests, roosts, and other central places that are often associated with important life history events and may exhibit unique characteristics; however, using satellite telemetry data to study central places is complicated by common nuances like locational error and animal movement. We coupled a novel modeling framework that accounts for these nuances with an Argos satellite telemetry dataset to examine the spatiotemporal behavior associated with harbor seal haul-out sites on Kodiak Island, Alaska, USA. The methodology incorporates an observation model that accommodates multiple sources of uncertainty in telemetry data and a flexible Bayesian nonparametric model to uncover latent clustering in the telemetry locations. We also contribute extensions to examine the effect of covariates on site selection and to obtain populationlevel inference concerning central place use. Harbor seal haul-out sites generally occurred in inlets and bays, areas that are isolated from the open water of the Gulf of Alaska. Most individuals selected haul-out sites that were protected from wave exposure. The effects of bathymetry and shoreline complexity on haulout site selection were variable among individual seals, as were the effects of time of day, time since low tide, and day of year on temporal patterns of haul-out use. As repositories of satellite telemetry data on a wide variety of species accumulate, so do opportunities for using this information to learn about the locations of central places, as well as the temporal patterns in their use. The model-based approach we describe offers a practical and rigorous means for gaining insight concerning these sensitive locations, knowledge of which is important for the effective management and conservation of many species.

movements, and delineate home ranges. These data also contain valuable information concerning dens, nests, roosts, and other central places that are often associated with important life history events and may exhibit unique characteristics; however, using satellite telemetry data to study central places is complicated by common nuances like locational error and animal movement. We coupled a novel modeling framework that accounts for these nuances with an Argos satellite telemetry dataset to examine the spatiotemporal behavior associated with harbor seal haul-out sites on Kodiak Island, Alaska, USA. The methodology incorporates an observation model that accommodates multiple sources of uncertainty in telemetry data and a flexible Bayesian nonparametric model to uncover latent clustering in the telemetry locations. We also contribute extensions to examine the effect of covariates on site selection and to obtain populationlevel inference concerning central place use. Harbor seal haul-out sites generally occurred in inlets and bays, areas that are isolated from the open water of the Gulf of Alaska. Most individuals selected haul-out sites that were protected from wave exposure. The effects of bathymetry and shoreline complexity on haulout site selection were variable among individual seals, as were the effects of time of day, time since low tide, and day of year on temporal patterns of haul-out use. As repositories of satellite telemetry data on a wide variety of species accumulate, so do opportunities for using this information to learn about the locations of central places, as well as the temporal patterns in their use. The model-based approach we describe offers a practical and rigorous means for gaining insight concerning these sensitive locations, knowledge of which is important for the effective management and conservation of many species.

INTRODUCTION
Satellite telemetry data provide an important source of animal distribution information and are commonly used to quantify habitat selection (Johnson et al. 2008b, Johnson et al. 2013, examine movements (Jonsen et al. 2003, Johnson et al. 2008a, McClintock et al. 2012, and delineate home ranges (Kie et al. 2010). These data also contain information concerning another key aspect of an animal's habitat use, namely the behavior associated with central places, or locations that are used repeatedly through time such as dens, nests, roosts, and rendezvous sites (Anderson and Lindzey 2003, Knopff et al. 2009, Brost et al. 2017. As repositories of satellite telemetry data on a wide variety of species accumulate, so do opportunities for using this information to learn about the locations of central places and the temporal patterns in their use. Yet, nuances like locational error (i.e., deviations between the true animal locations and the recorded telemetry locations) and animal movement require careful consideration in such applications . In this paper, we couple a novel modeling framework that addresses these challenges with Argos satellite telemetry data to investigate the spatiotemporal haul-out behavior of a marine top predator, the harbor seal (Phoca vitulina richardii).
Harbor seals are widely distributed in the temperate and arctic waters of the Northern Hemisphere (Scheffer and Slipp 1944). Harbor seals regularly leave the water and haul out on beaches, intertidal areas, and glacial ice in tidewater fjords to rest, molt, escape aquatic predators, give birth, and rear their pups (Ling 1984, da Silva and Terhune 1988, Thompson 1989, Watts 1992. The locations of haul-out sites may change seasonally to track available food sources (Montgomery et al. 2007, Cunningham et al. 2009), although harbor seals exhibit high levels of site fidelity over months to years and typically return to the same haul-out sites between at-sea foraging bouts (Harkonen and Heide-Jørgensen 1990, Thompson et al. 1997, Cunningham et al. 2009, Womble and Gende 2013.
Knowledge of harbor seal haul-out behavior is necessary to effectively manage and conserve this species. Human-caused disturbances can flush harbor seals from their haul-out sites, resulting in increased energy expenditure and potentially decreased fitness (Suryan and Harvey 1999, Jansen et al. 2010, Cordes et al. 2011, Blundell and Pendleton 2015. In particular, the impact of tourism vessels (i.e., cruise ships) on hauled-out harbor seals has received much recent attention (Jansen et al. 2010, Young et al. 2014, Blundell and Pendleton 2015, Mathews et al. 2016, and mitigation efforts require an understanding of harbor seal haul-out behavior. Furthermore, harbor seal population monitoring, which relies on counts of seals ashore to estimate abundance, benefits from additional information concerning when and where seals haul out of the water (e.g., to estimate and maximize detection probability during aerial surveys; Boveng et al. 2003, Small et al. 2003, Ver Hoef and Frost 2003. Although harbor seal haul-out behavior has been studied via direct observation (e.g., Cordes et al. 2011), remote locations and harsh conditions often inhibit extensive field investigations. In Alaska, for example, current knowledge concerning the location of haul-out sites is limited to aerial surveys conducted during August and September , Small et al. 2003, Ver Hoef and Frost 2003. Existing satellite telemetry datasets provide a practical means for learning about harbor seal haul-out behavior throughout the full annual cycle, without the expense of conducting additional fieldwork.
We examined an Argos satellite telemetry dataset, collected from harbor seals monitored near Kodiak Island, Alaska, USA, to understand the spatial distribution of haul-out sites used by individuals, as well as the temporal patterns in haul-out use. We use a general model-based approach that rigorously accommodates multiple sources of uncertainty in telemetry data while simultaneously uncovering latent clustering representative of central places ( Fig. 1; Brost et al. 2017). Ancillary behavioral data (e.g., accelerometer data) are further incorporated into the model to improve location estimation and quantify temporal patterns in central place use. Here, we extend the modeling framework to also examine the effect of covariates on haul-out site selection and to obtain population-level inference concerning haul-out use. Our study simultaneously offers insight into harbor seal haul-out behavior, provides important methodological advancements, and demonstrates how satellite telemetry data are useful for studying spatiotemporal patterns in central place use by marine and terrestrial animals.

Harbor seal telemetry data
Harbor seals were captured near haul-out sites on Kodiak Island, Alaska (Fig. 2), and equipped with satellite-linked depth recorders (SDRs; Wildlife Computers, Redmond, Washington, USA) during 1994 and 1995. Seals were sedated with a mixture of ketamine and diazepam administered intramuscularly, or with intravenous diazepam, and transmitters attached to the mid-dorsal surface using quick-setting epoxy as described by Lowry et al. (2001). The SDRs transmitted to Argos receivers onboard polar orbiting meteorological satellites, a system that uses the Doppler effect for geopositioning. The Argos least-squares positioning algorithm assigns each telemetry location to one of six quality classes based on the number of transmissions received during a satellite pass. In order of decreasing accuracy, the location quality classes are 3, 2, 1, 0, A, and B. The location quality classes have different error patterns and magnitudes, and some exhibit an x-shaped error distribution that has greatest error variance along the NW-SE and NE-SW axes (Costa et al. 2010, Douglas et al. 2012, McClintock et al. 2014, Brost et al. 2015. Argos location errors are often >10 km in magnitude, and in some cases >100 km (Costa et al. 2010, Brost et al. 2015, distances that exceed the typical extent of harbor seal movements , Cunningham et al. 2009). Locations obtained within 24 h of tagging were removed to mitigate tagging effects. All remaining locations were used in the analyses described below (e.g., no a priori filtering to exclude locations containing potentially large error).
The SDRs included a conductivity sensor that determined when the device was wet (low resistance) vs. dry (high resistance). These ancillary behavioral data were used as a proxy for haulout use. In other words, instances when the device was dry indicate the individual was out of the water at a haul-out site when the location was recorded, whereas locations collected while a.) b.) Fig. 1. Simulation demonstrating the signal of central place behavior in telemetry location data (gray closed circles). Because satellite telemetry devices record sequences of animal locations, repeated use of a site yields multiple telemetry locations collected at that site. Consequently, clusters of locations in mapped telemetry data are important indicators of a central place. (a) When telemetry location error is small (i.e., the observed telemetry locations are near the true animal locations), clusters of locations (large red open circles) are conspicuous and central places can be easily identified visually or by using a series of user-specified time and distance thresholds Lindzey 2003, Knopff et al. 2009). (b) When telemetry location error is large, a more sophisticated methodology is required because clusters are poorly defined and the number and location of central places are highly uncertain. Both figures are based on the same sequence of true locations, but vary in the magnitude of added locational error. the device was wet indicate the individual was at sea. The devices were programmed with a delay (10 consecutive readings at 45-s intervals) to prevent spurious wet/dry state transitions associated with splashing on a haul-out or short dry periods experienced by the sensor while a seal was surfaced in the water. Therefore, the wet/dry data reliably indicate a harbor seal's ❖ www.esajournals.org 4 June 2020 ❖ Volume 11(6) ❖ Article e03123 haul-out status. Although the transmitters were also equipped with pressure sensors for measuring dive depth, these data were not used to determine haul-out status. Tags were programmed to pause transmissions after an individual was hauled out for 6 h, but otherwise operated continuously (i.e., not duty-cycled).

Statistical notation
Let s ic t ð Þ s i;c;x t ð Þ; s i;c;y t ð Þ À Á 0 represent the pair of coordinates for an observed telemetry location collected at time t 2 T, where i indexes an individual harbor seal (i = 1, . . ., N) and c indexes Argos location quality class (c {3, 2, 1, 0, A, B}). Also let l i ðtÞ ðl i;x ðtÞ; l i;y ðtÞÞ 0 be the coordinates for the corresponding latent haul-out site. The matrices S i and M i comprise the locations of telemetry data and haul-out sites, respectively, for individual i (i.e., S i s ic t ð Þ; 8t f g and M i l i t ð Þ; 8t f g ). We denote the ancillary behavioral data as y i (t), where y i (t) = 0 indicates the telemetry device on individual i was wet at time t and thus the harbor seal was at sea, and y i (t) = 1 indicates the device was dry and the individual was at a haul-out site. We also denote the spatial domain within which haul-out sites can exist as S (e.g., the shoreline).

Haul-out site location estimation
We estimated haul-out site locations using a hierarchical model consisting of two general components, an observation model and a process model (Brost et al. 2017). The observation model, which accounts for potentially complex telemetry error structures and animal movement, was formulated as In Eq. 1, an observed telemetry location arises from a mixture of two multivariate normal distributions centered at l i (t) with variance-covariance a function of Argos telemetry measurement error (described by Σ ic orR ic ) and animal movement (described by U i (t)). The parameter p i (t) denotes the probability associated with the mixture components. The matrix Σ ic is parameterized to allow for various telemetry error structures: where r 2 ic quantifies measurement error in the longitude direction, a ic modifies r 2 ic to describe error in the latitude direction, and q ic describes the correlation between errors in the two directions (Brost et al. 2015. The matrixR ic is identical to Σ ic except for the off-diagonal elements which are Àq ic ffiffiffiffiffi a ic p . When q ic = 0, the model specification in Eq. 2 allows for circular (a ic = 1) and elliptical (a ic 6 ¼ 1) error patterns. Alternatively, the error covariance model allows for the x-shaped error pattern evident in Argos telemetry data when q ic 6 ¼ 0, where the mixture component with variance-covariance matrix Σ ic describes error along the SW-NE axis and the mixture component with variance-covariance matrixR ic describes error along the NW-SE axis. Note that telemetry error parameters are indexed by c, allowing the characteristics of locational error to vary for each Argos location quality class. We set p i (t) = 0.5 because the orbital plane of the Argos satellites changes continuously, and telemetry locations are equally likely to arise from either mixture component (i.e., locations have a 50% chance of arising from either axis of the x-shaped Argos error distribution; Brost et al. 2015).
When y i (t) = 0 and the harbor seal is at sea, uncertainty in the location of the haul-out site (l i (t)) is not only a function of telemetry error (Σ ic andR ic ), but also animal movement about the haul-out site. In this case, we let U i ðtÞ ¼ / 2 i I; where / 2 i describes the spread of the "home range" for individual i (l i (t) describes its center). Conversely, when y i (t) = 1, the true but unknown harbor seal location is the same as the location of the latent haul-out site. Uncertainty in l i (t) is thus exclusively due to telemetry error (Σ ic and R ic ), and we set Φ i (t) = 0, where 0 is a 2 9 2 matrix of zeros.
The process model used to estimate the true but unobserved locations of haul-out sites (l i (t)) consists of a clustering model known as a Dirichlet process (Brost et al. 2017). We represent the Dirichlet process as an infinite mixture ❖ www.esajournals.org where l ij are locations within the spatial domain S, d l ij is a point mass (or "atom") at l ij , p ij is the probability associated with the mixture component, and P 1 j¼1 p ij ¼ 1. The locations l ij and l i (t) distinguish possible haul-out sites from those actually used by a harbor seal. The l ij , for j ¼ 1; . . .; 1, are unique locations and represent the infinitely many possible sites (e.g., all locations along the shoreline) where a seal could haul out. Conversely, the l i (t) have a functional interpretation because they associate a l ij to each telemetry location s ic (t); they are locations that are used as a haul-out site. Put simply, Eq. 3 associates telemetry locations with haul-out sites probabilistically: each mixture component represents the location of a potential haul-out site, and p ij is the probability of the location being used by a harbor seal.
The number of mixture components necessary to generate the observed data can, at most, be the number of telemetry locations collected for an individual (e.g., a harbor seal never hauls out at the same location twice); however, only a handful of haul-out sites are actually used by any given seal. Accordingly, the mixture probabilities are formulated in a manner such that p ij decreases stochastically with increasing index j, favoring fewer haul-out sites with many locations per site (Sethuraman 1994, Ishwaran andJames 2001). In other words, as a discrete distribution that has probability mass concentrated over the first several mixture components, realizations of l i (t) from the Dirichlet process (Eq. 3) are generally not distinct, thereby creating clusters of telemetry locations that are defined by their association to a common haul-out site. Moreover, the rate of decrease in the p ij , and thus the number of haul-out sites used by a seal, is data-driven, allowing the complexity of the mixture to be tailored to each individual seal.
We used a custom Markov chain Monte Carlo (MCMC) algorithm written in R (R Core Team 2015) to estimate the observation and process model parameters in a unified framework. Markov chain Monte Carlo is an iterative approach to obtaining random draws, or samples, from the posterior distribution of the unknown parameters (e.g., r 2 ic , U i , and l i (t); Gelfand and Smith 1990). The full model statement, including prior distributions for all unknown parameters, is shown in Appendix S1. Inference was based on 100,000 MCMC samples after convergence (50,000 samples were discarded as burn-in).

Haul-out site selection
Replication in the values of l i (t), for t 2 T, not only partitions telemetry locations into clusters, but it also provides a measure of intensity of use. That is, sites associated with more telemetry locations receive more use. We relate intensity of use to the environmental characteristics of the haulout sites themselves with a simple and computationally tractable multiple (or "process") imputation procedure (Hooten et al. 2010, Hanks et al. 2011, Scharf et al. 2017. Conditioned on the locations of all haul-out sites (M i ), we obtain inference concerning haul-out site selection using the posterior distribution where the notation [Á] represents a probability distribution, b i is a vector of individual-level parameters quantifying the effect of covariates on haulout site selection, l b is a vector of population-level parameters that describe the mean effect across all harbor seals, and M i may be transformed by some deterministic function g(Á). Although application of Eq. 4 is straightforward, it does not account for uncertainty in our knowledge of l i (t). Therefore, we seek the posterior distribution of b i and l b given the observed telemetry locations recorded for individual i (i.e., S i ): Within a Bayesian hierarchical modeling framework, we perform the integration in Eq. 6 through composition sampling. That is, we make inference about the environmental drivers of haul-out site selection (i.e., b i and l b ) by conditioning on realizations from the posterior distribution of l i (t) (see below for specific procedural details; Hooten et al. 2010, Hanks et al. 2011). In the context of customary multiple imputation applications for missing data problems, the distribution [g(M i )|S i ] is akin to the imputation distribution (Scharf et al. 2017).
In practice, we constructed a model for w i = g (M i ), where the function g aggregates M i to obtain a counting process of the latent haul-out sites over raster cells in the spatial support S. The hierarchical mixture model accommodates zero inflation and variability among individual seals: where w ij is the number of times in which raster cell j is used as a haul-out site by individual i (i.e., w ij ¼ P t2T 1 fl i ðtÞ¼l ij g , where the indicator function 1 {arg} equals 1 if arg is satisfied and 0 otherwise), k ij is the intensity of the Poisson distribution, and z ij is a latent indicator variable that specifies the mixture component from which w ij arises. The mixture component consisting of a point mass at 0 (i.e., when z ij = 0) accounts for the preponderance of raster cells within S that are never used as a haul-out site (i.e., more instances of w ij = 0 than expected under the Poisson distribution alone; Welsh et al. 1996, Martin et al. 2005. We model the mean as a function of environmental characteristics using where x ij is a vector of covariates associated with raster cell j. The hierarchical specification is completed with individual-and population-level models: where Σ b and r 2 l b I are the respective covariance matrices. The multiple imputation procedure described above is implemented by sampling gðM ðkÞ i Þ $ ½gðM i ÞjS i on the kth iteration of an MCMC algorithm used to estimate the parameters in Eqs. 7-10, which are subsequently updated conditional on the value for gðM Appendix S2 for the full model statement, prior specifications, and details regarding model implementation. We investigated the effect of three covariates on the locations of haul-out sites, namely water depth, wave exposure, and shoreline complexity. Bathymetric data in the form of depth soundings from National Oceanic and Atmospheric Administration Electronic Navigation Charts (nauticalcharts.noaa.gov) were used to examine the relationship between estimated haul-out site locations and distance to water depth. We converted the depthsounding data (points) to a 100-m resolution raster and calculated the distance from each cell in S to the closest raster cell with water depth 20 m or greater (Montgomery et al. 2007). All distances were calculated as least-cost distance such that measurements were made exclusively through the water, thereby reflecting distances "as the seal swims" (distance measurements did not cross land; Dijkstra 1959). Biological wave exposure, obtained from ShoreZone aerial surveys (shorezone.org), was used to determine whether harbor seals selected certain exposures more than others. Biological wave exposure is assigned based on the presence and abundance of coastal biota that have known wave energy tolerances (Harper and Morris 2014). We combined the six biological wave exposure categories into two broader classes: protected (very protected, protected, and moderately protected categories) and exposed (very exposed, exposed, and moderately exposed categories). Biological wave exposure is considered a better index of exposure than physical wave exposure, which is based on fetch and coastal geomorphology (Harper and Morris 2014). We calculated shoreline complexity as the number of raster cells in S within 5 km of a focal cell. Thus, raster cells surrounded by more circuitous sections of shoreline have higher values for shoreline complexity. All inference was based on 100,000 MCMC iterations.

Temporal patterns in haul-out use
We obtained individual-and population-level inference concerning temporal patterns in haulout use with a hierarchical model that accommodates the simultaneous analysis of multiple individuals . Specifically, we modeled the behavioral data using a binary regression ❖ www.esajournals.org where y i (t) is the haul-out status of individual i at time t and w i (t) is the corresponding probability of being hauled out. We used the probit link to relate w i (t) to environmental conditions: where u i (t) are covariates measured at time t, c i are the corresponding individual-level coefficients, and Φ is the standard normal cumulative distribution function. In contrast to the more common logit link used in logistic regression, the probit link streamlines computation when fitting the model using MCMC (Albert and Chib 1993, Hooten et al. 2003, Dorazio and Rodr ıguez 2012, Johnson et al. 2012. The individual-level parameters were further modeled using where l c is a vector of population-level parameters that represent the average effect across all individuals, and Σ c and r 2 l c I are the respective variance-covariance matrices. The full model statement and details pertaining to model implementation are provided in Appendix S3. Previous studies indicated that temporal patterns in haul-out use are influenced by behaviors (e.g., breeding and foraging), physiological functions (e.g., thermoregulation and molting), and environmental conditions (e.g., tidal state) that operate at varying time scales , London et al. 2012. Accordingly, we examined the effect of the following covariates on harbor seal haul-out use: the number of hours from the nearest solar noon (i.e., 13:00 h near Kodiak Island), the number of hours from the nearest low tide, the number of days since 15 August, and the quadratic effect associated with days since 15 August. 15 August is approximately the beginning of the annual molting period in Alaska when harbor seals are most likely to be hauled out (Calambokidis et al. 1987, Daniel et al. 2003. Tide information was obtained from the National Oceanic and Atmospheric Administration water level monitoring stations nearest to the locations of monitored seals (tidesandcurrents.noaa.gov; Kodiak Island, Station ID: 9457292 and SW Terror Bay, Station ID: 9457493). Inference was based on 100,000 MCMC samples from the posterior distributions after convergence (50,000 samples were discarded as burn-in).

Harbor seal telemetry data
Twelve harbor seals were telemetered between October 1994 and June 1996, including six males and six females (Appendices S4 and S5). The age composition of harbor seals at the time of capture was seven adults, three subadults, and two pups. The average duration seals were monitored was 183 d (range: 76, 261 d), the average number of telemetry locations per individual was 760 (range: 301, 1460 locations), and the average time elapsed between telemetry locations was 5.8 h (range: 0, 499 h). On average, 37% of locations were recorded while an individual was hauled out of the water (range: 17, 61%).
Overall, 82% of the telemetry locations belonged to the Argos location quality classes that typically have the largest locational errors (classes 0, A, and B). Poor quality location classes were more likely to be associated with telemetry locations recorded while a harbor seal was at sea. For example, 65% of at-sea locations belonged to Argos location classes A and B, whereas only 35% of locations recorded while a harbor seal was hauled out of the water belonged to the same classes. Conversely, 31% of telemetry locations recorded while a harbor seal was hauled out belonged to the highest quality location classes (3, 2, and 1); only 12% of at-sea telemetry locations belonged to these higher quality classes.

Haul-out site location estimation
Inference concerning the intensity of haul-out site use (i.e., l i (t)) for all 12 harbor seals is summarized in Fig. 2. High posterior probability of haul-out use typically occurred in inlets and bays that were isolated from the open ocean (Fig. 2bd). One exception, however, was a subadult female that had high posterior probability of hauling out on an islet at the southeast corner of Kodiak Island, a location adjacent to the Gulf of Alaska (Fig. 2e). Inference concerning haul-out use for each individual harbor seal is presented in Appendix S4, and estimates for parameters in the observation model (i.e., parameters related to Argos telemetry error and animal movement) are presented in Appendix S5.

Haul-out site selection
Individual-level coefficients indicate a heterogeneous response to distance to 20 m bathymetric depth and shoreline complexity; responses were not consistent among individuals within sex or age classes (Fig. 3). Inference concerning l b , the population-level parameters, reflects this heterogeneity and suggests the individual harbor seals we examined lacked a common behavior relative to these two covariates (i.e., 95% credible intervals overlap 0; Fig. 3).
We were unable to evaluate the effect of wave exposure on haul-out site selection for most individuals due to complete separation in the counts w ij and the two exposure categories Anderson 1984, Hefley and. Complete separation occurred when all instances of w ij > 0 (or w ij = 0) for an individual were allocated to one of the exposure categories (protected or exposed), resulting in Markov chains that exhibited poor mixing and failed convergence. Instances of w ij > 0 were allocated exclusively to the "protected" category for six individuals, whereas all instances of w ij > 0 were allocated to the "exposed" category for one individual (shown in Fig. 2e). A model of the subset of five individuals for which complete separation did not occur indicated that haul-out site selection was negatively affected by "exposed" shorelines for three individuals; the remaining two individuals exhibited no effect for this covariate (see Appendix S6). Furthermore, the average effect of "exposed" shoreline was negative for these five harbor seals (95% CI for l b : À4.55, À0.20).

Temporal patterns in haul-out use
The effect of environmental conditions (time of day, time since low tide, and day of year) on temporal patterns of haul-out use was variable among the 12 harbor seals we examined (Fig. 4). A consistent pattern among individuals within sex or age classes was also not evident (Fig. 4). Distance to 20-m depth -4 -2 0 2 Shoreline complexity Fig. 3. Individual-and population-level inference concerning covariates examined in the haul-out site selection model. The top row (blue box) represents inference concerning the population-level parameter (l b ) that describes the average effect across the 12 harbor seals analyzed. The remaining rows show individual-level parameters (b i ), and individual seals are labeled according to their sex and age class. The points indicate the posterior mean, the thick lines represent the 50% credible interval, and the thin lines represent the 95% credible interval. was highly uncertain, and 95% credible intervals for most individuals included 0. Inference concerning the population-level parameters reflects the individual-level heterogeneity and shows that a common effect across individuals was generally lacking (i.e., 95% credible intervals for l c overlap 0; Fig. 4), although there appears to be a weak negative effect for the covariate hours since solar noon (95% CI: À0.37, 0.04).

DISCUSSION
We combined a previously collected Argos satellite telemetry dataset with a fully modelbased framework to examine the spatiotemporal behavior associated with harbor seal haul-out sites on Kodiak Island, Alaska. We adopted a recently developed approach to estimate the location of central places that rigorously accommodates large telemetry location error and animal movement (Brost et al. 2017), but also contribute extensions to examine the effect of environmental covariates and obtain populationlevel inference during the simultaneous analysis of multiple individuals. Our methods are general and can be applied to various types of telemetry data collected on terrestrial or marine species.
We used a multiple imputation procedure to relate haul-out site selection to shoreline characteristics (Hooten et al. 2010, Hanks et al. 2011), a computationally efficient approach that allows the site selection parameters (i.e., b i and l b ) to reflect uncertainty in the estimated location of the haul-out sites themselves (Eqs. 4-10). A model for haul-out site selection could, in principle, alternatively be specified in the base distribution of the Dirichlet process (i.e., the prior for l ij ; Hjort 2010). Such an approach allows habitat characteristics to directly inform the location of potential haul-outs using, for example, a point process that has a spatially heterogeneous intensity function (Aarts et al. 2012, Brost et al. 2015; however, this strategy has a couple of key shortcomings. As a model on l ij , the base distribution provides inference concerning whether a particular location is used as a haul-out and thus fails to distinguish between sites used a single time from those used repeatedly. Moreover, relatively few locations are used as a haul-out site by any given seal (i.e., <10), limiting the number of explanatory variables that can be incorporated into a model of the base distribution (Agresti 2002, pp. 212). Instead, we prefer to model a function of l i (t) (i.e., w ij = g(M i )), a strategy that examines how much a site is used and provides inference analogous to that obtained from common species distribution models (Aarts et al. 2012).
Our analysis couples positional data with behavioral (wet/dry) data to help inform a model of haul-out site location estimation. We could have also formulated a simpler model by exclusively using telemetry locations collected while harbor seals were hauled out of the water (i.e., y i (t) = 1) and uncertainty due to animal movement does not degrade inference; however, nearly two-thirds of the locations in our dataset were recorded while individuals were at sea (i.e., y i (t) = 0). These "at-sea" locations contain valuable information concerning the true location of the haul-out sites, making it important to incorporate both behavioral states using an observation model (Eq. 1) that explicitly accounts for animal movement. Other behavioral information, such as accelerometer data, could also be used to partition when individuals were using coastal vs. at-sea resources, or the model can be adapted to situations when no such ancillary data are available (e.g., Brost et al. 2017).

Harbor seal haul-out behavior
Almost all existing information concerning the location of harbor seal haul-out sites in Alaska has been acquired from aerial surveys that are typically conducted during the molting season in August and September , Small et al. 2003, Ver Hoef and Frost 2003, Womble et al. 2010. Satellite telemetry data, collected throughout the full annual cycle, provide an opportunity to learn about the location of haulout sites used at other times of year. Our results suggest that harbor seals favor haul-out sites in isolated bays and inlets (Fig. 2); however, one harbor seal in our study used a haul-out site on an islet adjacent to the open water of the Gulf of Alaska (Fig. 2e). Locations determined to have high posterior probability of haul-out use in our study match the locations of haul-out sites observed during aerial surveys conducted between 1993 and 2001 (see haul-out sites 7-10 in Figure 4 of Small et al. 2003 and Figure 1 of Boveng et al. 2003), serving as an informal evaluation of the methodology we used.
We are aware of only one previous study that examined the relationship between environmental variables and the selection of haul-out sites by harbor seals. Montgomery et al. (2007) used counts of harbor seals obtained from aerial surveys to model terrestrial habitat use in Cook Inlet, a 20,000-km 2 tidal estuary that opens into the Gulf of Alaska < 100 km north of Kodiak Island. They found that abundance of harbor seals was negatively related to distance to Cook Inlet communities, bathymetric depths of 20 m, and anadromous fish streams and that harbor seals tended to use haul-out sites with a rock substrate. Montgomery et al. (2007) also . Individual-and population-level inference concerning covariates examined in the temporal haul-out use model. The top row (blue box) represents inference concerning the population-level parameter (l c ) that describes the average effect across the 12 harbor seals analyzed. The remaining rows show individual-level parameters (c i ), and individual seals are labeled according to their sex and age class. The points indicate the posterior mean, the thick lines represent the 50% credible interval, and the thin lines represent the 95% credible interval. examined wave exposure but did not find evidence supporting an effect for this covariate. Our results do not indicate an effect of proximity to a bathymetric depth of 20 m and show a heterogeneous response to shoreline complexity (Fig. 3). Half of the harbor seals in our analysis exhibited complete separation such that estimated locations of haul-out sites only occurred on "protected" shoreline. This relationship could be coincidental or reflect selection for "protected" shorelines. An analysis of the subset of harbor seals for which complete separation did not occur revealed that three of five individuals selected against "exposed" shorelines (see Appendix S6). The areas in which seals were monitored for this study are isolated from the human communities on Kodiak Island and do not contain substantial variation in shoreline substrate (the majority of shoreline was sedimentary or of mixed types as determined by ShoreZone aerial surveys), and data concerning seasonal variation in prey availability were lacking. Therefore, we did not examine the effect of communities, shoreline substrate, or proximity to fish streams on haul-out use.
Several studies have used counts of harbor seal haul-out groups to investigate patterns in haulout use at multiple temporal scales. At daily time scales, the highest proportion of seals onshore are typically observed at times near low tides when favorable haul-outs are exposed (e.g., sites isolated from terrestrial predators; Schneider andPayne 1983, Pauli andTerhune 1987) and during midday when the air temperature is most conducive to thermoregulation (Stewart 1984, Calambokidis et al. 1987, Pauli and Terhune 1987see London et al. 2012 for an exception to these patterns). At annual time scales, temporal patterns in haul-out use are influenced by breeding and molting cycles that can be sex-and agespecific (e.g., adult females nurse pups onshore, and pups do not molt; Everitt and Braham 1980, Brown and Mate 1983, Calambokidis et al. 1987, Huber et al. 2001, Jemison and Kelly 2001, Daniel et al. 2003, as well as the distribution and availability of prey. Our results do not show a consistent effect of environmental conditions on haul-out behavior and thus do not corroborate the conclusions of these earlier studies (Fig. 4); however, haul-out behavior is known to vary regionally because seals likely adapt their behavior to local conditions . Patterns related to sex and age class resulting from differences in foraging strategy (e.g., Frost et al. 2001, Hastings et al. 2004, Carter et al. 2017, for example, were also not evident in our study.
Our study focused on the telemetry data available in an important region of the Alaskan coastline and included a balanced sex ratio of harbor seals. The model can be scaled up to include more individuals and larger spatial scales when additional data and computing resources are available; as with all studies, these steps increase the ability to identify demographic effects and reach more general conclusions. Although our analysis concerns a relatively small subpopulation of harbor seals, it also provides a solid and intuitive procedural framework that can be extended to larger populations of marine or terrestrial animals that use central places. For example, examining the central place use of many individuals of a variety of species would provide inference not obtainable through traditional space use and resource selection approaches.

CONCLUSION
Previously, central places were often studied using portable radiotelemetry (i.e., VHF) equipment to locate animals in the field; however, use of satellite telemetry has become standard practice in animal tracking studies, precluding "onthe-ground" searching as a means for learning about central places. Nevertheless, satellite telemetry data contain valuable information concerning these sensitive sites, which in turn can aid spatial planning for species of conservation concern (e.g., Cunningham et al. 2009, Womble andGende 2013). The model-based approach we describe offers a practical and rigorous means for gaining additional insight concerning the spatiotemporal behavior associated with central places. Future work could focus on incorporating temporal dynamics in the spatial model for central places (e.g., to examine seasonal patterns in haul-out site locations) or on the development of a hierarchical Dirichlet process for the concurrent analysis of multiple individuals.