A Bayesian Analysis of Abundance, Trend, and Population Viability for Harbor Seals in Iliamna Lake, Alaska

Harbor seals in Iliamna Lake, Alaska, are a small, isolated population, and one of only two freshwater populations of harbor seals in the world, yet little is known about their abundance or risk for extinction. Bayesian hierarchical models were used to estimate abundance and trend of this population. Observational models were developed from aerial survey and harvest data, and they included effects for time of year and time of day on survey counts. Underlying models of abundance and trend were based on a Leslie matrix model that used prior information on vital rates from the literature. We developed three scenarios for variability in the priors and used them as part of a sensitivity analysis. The models were fitted using Markov chain Monte Carlo methods. The population production rate implied by the vital rate estimates was about 5% per year, very similar to the average annual harvest rate. After a period of growth in the 1980s, the population appears to be relatively stable at around 400 individuals. A population viability analysis assessing the risk of quasi‐extinction, defined as any reduction to 50 animals or below in the next 100 years, ranged from 1% to 3%, depending on the prior scenario. Although this is moderately low risk, it does not include genetic or catastrophic environmental events, which may have occurred to the population in the past, so our results should be applied cautiously.


INTRODUCTION
For the management of a natural population of animals, the goals of analysis often include learning about (1) demographic parameters, (2) the population's relationship to the environment, (3) abundance of the population, (4) trend of the population, and (5) viability of the population. All of these goals may be important when making decisions about a population's risk of extinction, such as a listing decision Marine Mammal Laboratory, Alaska Fisheries Science Center, National Marine Fisheries Service, National Oceanic and Atmospheric Administration, Seattle, WA, USA.
* Address correspondence to Peter L. Boveng under the U.S. Endangered Species Act. Often, different data sets and analyses are used for the different goals. However, information on any component is typically expensive to collect and may not be readily available. In those cases, in order to make the best decision, it may be possible to develop prior distributions on various components from other literature and expert opinion. In addition, demography, environment, abundance, trend, and viability are all related, and a single model that borrows information from each component and treats uncertainty globally is preferred (Goodman, 2004). Although not the only approach, Bayesian models are attractive in this regard and have been used successfully in the past as integrated analyses (Goodman, 2002;Maunder, 2004;Wade, 2002). In this article, we develop a Bayesian model for incorporating relatively sparse information on counts, harvests, demography, and environmental factors to improve quantitative understanding of abundance, trend, and population viability for a small population of harbor seals, Phoca vitulina, inhabiting a large freshwater lake in Alaska. Iliamna Lake is located in southwest Alaska, and is drained by the Kvichak River to Bristol Bay. It is Alaska's largest freshwater lake, approximately 121 km long, up to 32 km wide, and as deep as 300 m (Anderson, 1969;Hauser, Allen, Rich, & Quinn, 2008). The seals in Iliamna Lake have been genetically confirmed as harbor seals, and there has been no scientific documentation of any rationale for considering them to be a subspecies other than the eastern North Pacific subspecies (P. v. richardii) (Burns et al., 2016). These are one of just five freshwater seal populations in the world, and one of just two harbor seal populations that persist in freshwater, the other consisting of the subspecies P. v. mellonae, in the Lacs des Loups Marins of the Ungava Peninsula, Québec (Smith, Hobson, Koopman, & Lavigne, 1996). The seals in Iliamna Lake appear to be discrete, though there could potentially be interchange with the saltwater seal population in Bristol Bay via the Kvichak River. Genetic analysis and local understanding of the lake seals' isolation suggest that they are a separate population (Burns et al., 2016;National Marine Fisheries Service, 2016).
Seals are observed in the lake year-round, and the peak of births occurs in mid-July (Burns et al., 2016). The harbor seals in Iliamna Lake make use of abundant salmon runs (Hauser et al., 2008), at least during the summer, as their main prey. Iliamna Lake typically freezes over for much or all of the winter and provides a unique challenge for year-round persistence. When ice covers the lake, fewer seals are observed and little is known about how seals adjust their behavior in response to the ice cover. During that time, harbor seals in the lake likely make use of near-shore leads, cracks, and spaces under the ice for breathing (Burns et al., 2016). In the winter, seals are occasionally observed resting on the ice, and harvest by Alaska Natives usually occurs on the ice in spring, before the ice cover breaks up.
The conservation and management of harbor seals in Iliamna Lake requires decisions that must be made with limited empirical data on demographic parameters, abundance, trend, and population viability. Additionally, the data that are available originated from a variety of research projects and methodologies that were not necessarily designed for our objectives. We used aerial surveys of harbor seals in Iliamna Lake conducted by various agencies and or-ganizations, and harvest information, also collected from a variety of sources. Our Bayesian analysis allows flexible inclusion of data from varied sources and allows us to borrow useful information from other harbor seal studies. This framework not only provides the basis for a population viability analysis but sheds light on this unusual ecological setting for harbor seals. We now summarize some of the literature that either informed our development of priors for the Bayesian model or expresses our expectations about the results.
We examined two main factors associated with variation in observed seal numbers: time of year and time of day for aerial surveys. These factors are generally related to when seals haul out of the water. From studies of seals in saltwater, there is evidence that seals prefer to haul out near the middle of the day, or shortly thereafter Mathews & Pendleton, 2006;Ver Hoef, London, & Boveng, 2010). However, many other factors are also important, so the time-of-day effect may be weak and complex with considerable variation (London, Ver Hoef, Jeffries, Lance, & Boveng, 2012;Ver Hoef & Frost, 2003). Nonpup seals undergo molt each year, generally during the month of August in Alaska. It is well established (see the literature listed above) that harbor seals in marine environments spend more time hauled out during molt, so we expected to see similar patterns for the harbor seals of Iliamna Lake, with a peak sometime in early to mid-August. We expected to see an earlier peak for pups, with peak pupping in mid-July. However, pups do not molt in their first summer and will start to spend more time foraging (less time hauled out) as they grow. In addition, the pups grow rapidly and are difficult to distinguish from adults by mid-August when they wean, so we expected their counts to decrease by the time nonpup counts were peaking.
No formal abundance estimates have been made for the harbor seals in Iliamna Lake. Anecdotal evidence (Burns, 1978) suggested that there may have been around 300 individuals in the 1960s, and a subsequent drop to about 50 after a series of hard winters. Aerial surveys commenced in 1984 and have been conducted more frequently in recent years. Because harbor seal counts typically vary greatly from day to day, Mathisen and Kline (1992) suggested summing the highest counts from each haul-out site from multiple surveys in a short time interval, yielding an estimate of 137 seals in the lake in 1991. Later surveys indicated greater abundance to near 400 seals, with a maximum single-day count of 357 in 2008 (ABR, Environmental Research and Service, 2011). Summing the highest site-wise counts as recommended by Mathisen and Kline (1992) produced estimates of 367, 368, and 389 seals in 2005, respectively. All estimates lacked an explicit correction for seals that were missed because they were in the water. In saltwater, Simpkins et al. (2003) estimated that at least 20-35% of seals would be in the water, even when correcting to hypothetical ideal conditions (of date, time of day, tide, and weather) for hauling out. Using just date, time of day, and tide information from satellite transmitters affixed to seals in nearby Cook Inlet, Alaska, Ver Hoef et al. (2010) found that about 40% of seals would be hauled out under optimum conditions. Those data lacked information from the molt period because satellite transmitters were attached to hair, which falls off during molt. In summary, previous estimates based on counts are likely to be low, by perhaps 20% or more. Thus, our working hypothesis was that seal abundance was slightly over 400 individuals in 2005-2009. With so little documentation of historical abundance, we could only surmise that the population trend may have encompassed an increase since the 1970s, from a low of around 50 seals to recent numbers of around 400. In addition to using observed changes in abundance to estimate trends, rates of survival and reproduction can be used to infer population trend. However, even less is known about these demographic parameters than abundance for harbor seals in Iliamna Lake. In fact, very little is known about demographic parameters of harbor seals along the Pacific Rim. An exception is a detailed study on survival for various sex and age classes at a nearby saltwater site, Tugidak Island, by Hastings, Small, and Pendleton (2012). Little is known about reproduction either, though Pitcher and Calkins (1979) provided data from females in the greater Gulf of Alaska region. Local residents along the lake perceived recent seal numbers to be relatively constant (Burns et al., 2016), so we predicted that new data and analysis would indicate a stable population.
Even a population that has been stable faces certain risks if it is small. Four main risks for species extinctions were outlined by Shaffer (1981): genetic effects in small populations, demographic stochasticity, environmental stochasticity, and catastrophic events. All of these apply to populations of any size, but genetic effects and demographic stochasticity are especially important in small populations (Goodman, 1987). A general rule of thumb was established by Franklin (1980) that a population size greater than 50 was necessary to avoid greatly increased chances of extinction due to genetic effects (but see Taylor & Rojas-Bracho, 1999 for caveats). After establishing such thresholds, one goal of population viability analysis is to assess the risk, or probability, of crossing such a threshold, often called quasi-extinction (Ginzburg, Slobodkin, Johnson, & Bindman, 1982), based on current understanding of a population and the ability to project it into the future. One-hundred years is a common time frame for assessing extinction risk of animals, which represents about 11 generations for harbor seals (Swart, Reijnders, & Van Delden, 1996). Based on the relatively small population size, we expected that the harbor seals in Iliamna Lake are at considerable risk for crossing into a realm of critically low population size.
The goal of this article was to develop a Bayesian hierarchical model based in part on a population process model that includes stochastic variation in vital rates and harvest of seals. Linked to the process model is an observational model for counts from aerial surveys that includes time-of-day and timeof-year effects, and the proportion of seals that are in the water and unavailable to be counted during surveys. From this model, we wished to obtain estimates of abundance and trend and project the population into the future to assess the probability of quasi-extinction for different time horizons up to 100 years. We also wished to perform sensitivity analyses to better understand our prior assumptions. The practical goal of obtaining these results was to provide information for sound management of the harbor seal population in Iliamna Lake.

DATA
We conducted aerial surveys of harbor seals in Iliamna Lake between 2008 and 2013. Surveys were conducted in collaboration with local community participants and researchers at the University of Alaska (Burns et al., 2016). Surveys were flown using high-wing, twin-engine aircraft (Aero Commander 680, 690 or de Havilland Twin Otter). Survey altitude was typically 330 m, and aircraft speed was typically 120 kt. These abundance surveys were performed several times yearly for most years (Table I). Surveys were timed so that one survey was conducted while the lake was mostly frozen (late March/early April), one during pupping (mid-July), and often several during the August molt, when the greatest numbers of seals typically haul out on shore and are visible to aerial observers. Over the years, a comprehensive list of all known haul-out sites was developed. These sites were searched extensively on each flight. Surveys were flown, weather allowing, in the mid to late afternoon, when the number of seals hauled out was expected to be highest. Aircraft flight tracks were recorded by GPS, and all seals sighted were digitally photographed using a SLR camera (Nikon D700) with zoom lens (Nikkor 80-400 mm). Time, date, latitude, longitude, and altitude were automatically embedded into the image metadata. The observer manually recorded weather and visibility information, though the surveys were typically conducted during good conditions, and these variables were not used in our analysis. The total number of seals at each haul-out site was counted from the digital photographs. Pups were determined by their smaller size and close proximity (<1 body length; either nursing or laying right next) to a larger seal. Pups were no longer recorded beyond about mid-August when many had been weaned. By that date the pups were large enough that they could not reliably be distinguished from yearlings and two-year-olds, which we lumped into the nonpup age class for our analysis. In 2009, researchers from the Newhalen village provided 10 additional flights in collaboration with our survey project, and similar techniques were used.
We also used counts obtained prior to 2008 by other organizations. Between 2005 and 2008, ABR, Environmental Research and Service (2011) conducted a series of aerial surveys for harbor seals in Iliamna Lake. In addition, earlier counts from surveys conducted by the Alaska Department of Fish and Game (ADFG) (Small, 2001) and surveys from 1984 to 1991 by Mathisen and Kline (1992) were incorporated into the analysis to expand the historical reach. Geographic coordinates were provided (or, when not provided, determined based on descriptions or physical maps) for each survey site, and these sites were compared and merged with locations identified by us. In some cases, sites in very close geographic proximity were combined into a single site. In these surveys, as well as our own, seals observed in the shallow water near the haul-out site (typically, within 50 m) were included in the counts. It was rare to observe seals in the water away from haul-out sites, and the very few of those that were seen were not included in the counts. A consistent, collated data set of all five survey projects served as the basis for our analysis. The data we used are archived with the National Centers for Environmental Information (Withrow, London, Yano, & Boveng, 2014) and available for download. Table I summarizes the aerial survey effort across the various agencies and researchers.
Seals in Iliamna Lake have cultural and nutritional importance to residents of the region (Burns et al., 2016) and are harvested by local Alaska Native hunters. Hunting primarily occurs in the late winter and early spring when the lake is still ice covered. As the ice thins and cracks form, the seals haul out along the cracks and are accessible to hunters. Some hunting may go on during summer, but these animals are taken more opportunistically while the earlier hunt is a directed subsistence activity. The numbers of seals harvested from Iliamna Lake have been estimated by household surveys of hunters in the local communities, conducted by the ADFG Subsistence Division (Fall, Holen, Davis, Krieg, & Koster, 2006;Krieg, Holen, & Koster, 2009;Burns et al., 2016). We used harvest estimates for the various years between 1992 and 2011 (Table II) that were summarized by Burns et al. (2016). We also used estimates from Burns et al. (2011), which were obtained separately from, but by similar methods to, the ADFG studies. We did not use values from 1982-1991 that were reported by Burns et al. (2016) because we were less certain about whether some of those harvests were from the Kvichak River, rather than Iliamna Lake. Note that because not all households in each community were surveyed, harvest data were estimated such that the values are not necessarily represented as whole numbers of seals per year.

ABUNDANCE, TREND, AND POPULATION VIABILITY
In this section, we describe the basic model for the data, and priors for the model parameters. We used Markov chain Monte Carlo (MCMC) methods for fitting the model (e.g., Gilks, Richardson, and  Spiegelhalter, 1996) which provide a sample from the posterior distribution of the parameters, such as survival and productivity. That sample also enables inference on derived quantities such as abundance, trend, and population viability.

The Model
For surveys in the earlier years when pups were not distinguished from nonpups, we denote the combined counts as Z t,i for the ith count in the tth year. Later, when pups were distinguished from nonpups, we denote pup counts as X t,i for the ith count in the tth year, and nonpup counts as Y t,i . Harvest data, as given in Table II, are denoted by H t,m for the tth year and the mth village. We used the following observational models for harbor seal counts in Iliamna Lake, where N (a, b) is a Gaussian (normal) distribution with mean a and variance b, and DU(·) is a discrete uniform where H t,m takes one of the values ν m = (v 1 , . . . , v C m ) with equal probability 1/C m for C m different classes. Here, C m is the number of different years with observations per village, and we simply allow for a discrete uniform distribution on the observed values for each village. We built our model hierarchically (Cressie, Calder, Clark, Ver Hoef, & Wikle, 2009). We assumed that there is some true underlying population value for pups, P t , and nonpups, A t , for each year t. However, some part of the population is in the water, away from the haul-out sites, and unobservable, so we denote the proportion that is observable as ψ. The observed values, P t and A t , are also affected by covariates; in this case day of year d t,i and hour of day h t,i . The expected values of the models in Equations (1) are set up as linear models, with regression parameters α d,1 , α d,2 , α h,1 , and α h,2 that control how day of year and time of day affect pup counts, and β d,1 , β d,2 , β h,1 , and β h,2 , which control how day of year and time of day affect nonpup counts. Notice that we could model some interaction between α d,• and β d,• because, as pups mature during the summer, they cannot be distinguished from nonpups. However, we have too little information and data for such estimates, and hence we leave α d,• and β d,• as free parameters to capture such interactions.
Note that the construction ψ P t and ψ A t leads to a nonidentifiability problem in the model. Any increase in P t can be counteracted by a decrease in ψ so that the product remains unchanged. We approached the problem as follows. First, we gave total abundance a lower bound by constrain- i to be greater than max{X i,t }/0.9 for every t, and like- i to be greater than max{Y i,t }/0.9 for every t, or their sum to be greater than max{Z i,t }/0.9, depending on how the data were collected in year t. Our justification for such a constraint follows from Simpkins et al. (2003). They found that, under optimal conditions, the expected proportion of seals hauled out of the water was around 0.815, but the maximum that was observed across time was between 0.84 and 0.92 (for two different study areas). Thus, we used a value of 0.9. Note that, for consistency with surveys conducted by others previously, we have included seals seen in shallow water near the haul-out sites as part of the fraction of the population that is "hauled out." This differs slightly from the convention used in most studies of marine harbor seals, in which only seals actually hauled out on shore are considered to be hauled out and observable.
The nonidentifiability problem is still unsolved if abundance is uniformly distributed from the lower bound to infinity while ψ is constrained between zero and one. To solve this, we used an informative beta distribution as a prior for ψ. With a proper anchoring of ψ, there is some information on the haul-out distribution from the surveys because the counts varied considerably. For example, if we could assume that the maximum number of counted seals reflected that 100% were hauled out, then a mean haul-out rate could be computed from all other surveys with less than the maximum count. The prior distribution helps anchor this distribution to realistic haul-out values, and we discuss our particular choice in the next section.
We also used a constant ψ rather than a temporally varying ψ t , for several reasons. First, the number of surveys per year varied considerably. Allowing ψ t to vary yearly would allow abundance to essentially "track" the yearly maximum, regardless of the number of surveys in that year. The haul-out time study of Simpkins et al. (2003) was based on many time periods, so using the maximum across all years was more appropriate. Second, the issue of nonidentifiability exists for the annual variability of ψ t as well, which would be confounded with the annual variability of P t and A t . While we do believe there is annual variability in ψ, because of the confounding we interpret ψ as an average across years and interannual variability is captured in the posterior variability of ψ.
Finally, note that we also tried formulating X t,i , Y t,i , and Z t,i in Equation (1) using negative binomial distributions, and then formulating Equation (2) as: However, we found that this multiplicative formulation suffered from overfitting and further confounding, where the regression coefficients α •,• and β •,• would attempt to model the zero counts with unrealistic fits that increased P t and A t and decreased ψ in order to obtain the fits. It seems that having ψ as a multiplicative effect, and α •,• and β •,• as additive effects in Equation (2), allows for better identification in the model. Equation (2) is imperfect in the sense that it is not constrained to be nonnegative, nor are the normal distributions used in Equation (1), even though our data are counts. We find this acceptable for the stability in model fitting, and it has no impact on our ability to project into the future based on the Leslie matrix model described next, and in Section 3.5. We embedded a very simple population model that creates a time series and a relation between the numbers of pups and nonpups. The model included density dependence (Caswell, 2001, p. 511) for fecundity, as this is likely where density dependence occurs for seals, where δ t−1 is an autoregressive parameter that, in a population model setting, is survival of nonpups from one year to the next; φ t−1 is fecundity, defined to be the proportion of nonpups that produce pups that survive to next year's survey period in mid-summer; c t−1 = exp(−ρ N t−1 /1000) is a densitydependent factor for fecundity, scaled by a parameter ρ, that shrinks φ t−1 toward zero as N t−1 = A t−1 + P t−1 increases; κ t−1 is the proportion of pups that survive from the survey period to become yearlings (and enter the nonpup age class); and H t−1 is the number of seals from the total population that were harvested in the previous year. In Equation (3), absent the harvest, we have used a very simple two-component Leslie matrix model with density dependence (Caswell, 2001, p. 511) matches the age class data that we have, namely, counts for pups and nonpups. The model is stochastic (Wade, 2002;Caswell, 2001, p. 452) because we allow δ t , φ t , and κ t to vary randomly with year t. Early Bayesian approaches to population modeling that included Leslie matrix models include Raftery, Givens, and Zeh (1995) and Poole and Raftery (2000). A review of embedding Leslie matrix models in a Bayesian inference framework, including stochastic demographic parameters, is provided by Gross, Craig, and Hutchison (2002) and Buckland, Newman, Fernández, Thomas, and Harwood (2007). Other approaches using stochastic demographic parameters include Cáceres and Cáceres-Saez (2013), who studied effective growth rates when demographic parameters are known, or can be assumed known, with random perturbations (whose distribution is part of the interest). They show how to obtain a stable population state with the effective growth rate. Our interest centered more on abundance and trend, and not an estimate or a long-term stable state, which is not very informative with just two age classes in Equation (4). Note that for harbor seal counts, population trends have often been modeled with loglinear models, such as Poisson or negative binomial regression (Mathews & Pendleton, 2006;Small, Pendleton, & Pitcher, 2003;Ver Hoef & Frost, 2003), where the trend was a linear or quadratic function of time in the mean of the model. In those regression models, the form of the trend component, and how it behaves near the end of the time sequence, largely determines Annual survival of pups to nonpup class Density dependent fecundity (4) ρ Scaling factor for density dependent fecundity (4) the trend into the future. There are advantages to demographic models, even very simple ones, especially for population viability analysis. For example, suppose that the trend in the data is slightly positive.
A regression approach will continue the trend, and a population would never go extinct. However, a demographic model has such a possibility, even for an overall positive trend, due to the stochastic nature of the parameters, which is a more realistic property of populations.

Priors
A summary of model parameters is given in Table III. To complete the model, we specified prior distributions on the stochastic components, (5, 5) 30 where i.i.d. means independent and identically distributed. Note that in order to start the population model, we needed initial values for A 1 and P 1 , so each of these was given a flat uniform prior distribution.
Most Bayesian analyses require some discussion of the priors (Goodman, 2010). We used noninformative priors when we could, such as for η, σ , α ·,· , β ·,· , P 1 , and A 1 . However, our sample sizes were quite small, and the form and shape of the prior distribution had a large effect on the other parameters. To investigate those effects, and present a sensitivity analysis, we developed three scenarios for the demographic parameters fecundity (φ), nonpup survival (δ), and pup survival (κ). A detailed development of these priors is given in the Appendix, and we give an overview here.
In Table IV, the "base-priors" scenario, was developed from demographic parameters of harbor seal survival and fecundity given in Hastings et al. (2012)  Note: Sampling variation alone (base-priors), variation broadened from consideration of other expected sources of variation (wide-priors), and noninformative variation (uniform-priors). These are described in the text and fully developed in the Appendix. The parameters nonpup survival (δ), fecundity (φ), and pup survival (κ) are components of Equation (3). The columns labeled a and b are the beta distribution parameters, derived using moment matching from the mean and variance in the previous columns (see the Appendix for details). The uniform distribution in the uniform-priors extends from 0 to 0.5 for φ t and 0.5 to 1 for κ t and δ t . and Pitcher and Calkins (1979), respectively. Survival estimates for both males and females, at two different stages of development for pups, and ages classes 1 − 3 and > 3, were obtained from Hastings et al. (2012). Because estimates of pup survival were provided at two different stages, we were able to shift a detailed demographic model based on sex and age at the time of birth to one based on the time of surveys, and collapse it to just two classes, pups and nonpups.
In doing so, we preserved the stable age distribution and population growth rate of the nonshifted and uncollapsed matrix model, and we derived variance estimates for all parameters. Then, because φ, δ, and κ were constrained to greater than zero and less than one, we matched moments of the estimates and their variances with a beta distribution to provide a prior distribution that met the required constraints (see the Appendix for details).
While the base-priors were derived from the best information that we had on demographic estimates, the associated variances represented sampling error and did not necessarily reflect the variability associated with annual variation in the demographic parameters. Also, the base-priors came from a separate, marine population of seals. The data can inform the posterior distributions of φ t , δ t , and κ t , but when sample sizes are small there is some danger in making the priors too narrow because then those priors will largely determine the results. Hence, for the "widepriors" scenario, we tripled the standard errors of "base-priors." Also, under base-priors, the expectation for population growth rate was around 4% annually. To be neutral on population growth rate, we also adjusted the prior means of φ t , δ t , and κ t so that population growth rate was 0. The details are provided in the Appendix and are summarized in Table IV. For a third scenario, "uniform-priors," we used constrained noninformative priors. For φ, it would be very unlikely for fecundity to exceed 50% of all nonpups, considering some of the nonpups are male and some are too young to produce pups. It would also be very unlikely for nonpup survival to be less than 50% given the few studies done on harbor seals in Alaska (Hastings et al., 2012) and elsewhere (Cordes & Thompson, 2014;Härkönen & Heide-Jørgensen, 1990;Mackey, Durban, Middlemas, & Thompson, 2008), where they were all near 90%. We used a similar constraint for survival of pups from age at survey to nonpup status. Given these constraints, we assumed flat prior distributions (Table IV). Notice that under uniform-priors, the distribution on population growth rate is near -5% annually.
Finally, we considered the prior distribution for ψ, the expected proportion of seals that are observable. Research in ocean environments along coastal Alaska showed that haul-out rates averaged from 41% to 50% (Pitcher & McAllister, 1981), with ranges from 20% to 80% along the west coast of North America (Huber, Jeffries, Brown, Delong, & Vanblaricom, 2001;Simpkins et al., 2003;Sullivan, 1979;Ver Hoef et al., 2010). Using a prior of Beta(5,5) concentrated most of the prior mass between 0.2 and 0.8. Additionally, the prior was used multiple times, one for each year t; otherwise, there was a tendency for P t and A t to increase without bound, driving the posterior distribution of ψ toward zero.

Fitting the Model
The model was fitted with Markov chain Monte Carlo (MCMC) methods using the Metropolis algorithm (Metropolis, Rosenbluth, Rosenbluth, Teller, & Teller, 1953), which was custom coded in R (R Core Team, 2016). We used a burn-in of 50,000 iterations, and then kept each 10th sample from the next 100,000 iterations, yielding 10,000 samples. From the hierarchical model described in Sections 3.1 and 3.2, we obtained a sample from the posterior distribution, [η, σ, α, β, P 1 , A 1 , φ, δ, κ, ρ, ψ|x, y, z, d, h], where, from Equation (1) The samples from all posterior distributions were checked for convergence, and no problems were found (Jones, Haran, Caffo, & Neath, 2006).
There were missing values. In particular, we had harvest information for just six years. For the years that were missing harvest data, for each iteration of the MCMC algorithm, we imputed one of the observed values, and this was done by village. Recall that harvest was modeled as a discrete uniform distribution DU m (1/C m ). For example, in Igiugig where harvest data were 2, 5.4, and 4.2, all missing years were assigned one of these three values, at random for each year, with probability 1/3. Then, the total harvest for each year was computed as the sum across the observed and imputed values for all villages. Likewise, there were a few seal counts that were missing hour-of-day values, and so for each iteration of the MCMC algorithm we randomly im-puted one of the observed values. These imputations preserved the variation in the observed data for the missing values without explicitly modeling these variables, which were of no direct interest.

Trend and Abundance
Let P [k] 1 and A [k] 1 be the kth MCMC sample of the posterior distribution of the true, initial abundance of pups and nonpups, P 1 and A 1 . Then abundance in the following years is determined by the MCMC sample of the demographic parameters and the observed and imputed harvest values, t depends on the kth MCMC sample, ρ [k] , along with A t and P t . The kth MCMC estimate of total abundance then is N t for any t. We used two methods to estimate trend. To obtain a global estimate of trend for the model, we pooled all of the leading eigenvalues of M [k] t across both t and k. As a second method, we used simple linear regression on the posterior distribution of the annual abundance estimates (Johnson & Fritz, 2014); that is, we simply regressed abundance on year. We did this for the last 5, 10, and 15 years. An estimate of trend for the most recent five years, taken from the kth MCMC sample, was computed as the slope of the regression of N [k] t ; t = 26, . . . , 30 on the integers 1, . . . , 5, which we denote as τ [k] 5 . Likewise, an estimate of trend for the most recent 10 years, taken from the kth MCMC sample, was computed as the slope of the regression of N [k] t ; t = 21, . . . , 30 on the integers 1, . . . , 10, which we denote as τ [k] 10 ; and an estimate of trend for the most recent 15 years, taken from the kth MCMC sample, was computed as the slope of the regression of N [k] t ; t = 16, . . . , 30 on 1, . . . , 15, which we denote as τ [k] 15 .
For each t = 30, . . . , 129, we resampled 1 year, t , independently from t = 1, . . . , 29, with replacement, and then set δ for t = 30, . . . , 129 and k = 1, . . . , 10, 000. We also extended harvest into the future. Initially, we drew H [k] t,m randomly from DU m (1/C m ) for the mth village, and then summed them,H [k] t = 6 m=1H t,m for t = 30, . . . , 129 and k = 1, . . . , 10, 000. Based on analyses given in the subsequent sections, it appears that the harbor seal population averaged around 350-400 seals, so we allowed for simulated harvest rates (rather than simulated amounts) by scaling harvest to current abundance: H is the current estimate of abundance, as described below. Thus, as the population increases, more seals would be harvested, and also as the population decreases, fewer seals would be harvested. This scaling seems natural and keeps the population from becoming negative due to harvest.
The simulated population abundance, projected into the future for each k, t, is: t , for t = 30, . . . , 129 and k = 1, . . . , 10, 000. For any t, k combination, the simulated total population is N t , the population will never reach 0. In this case, a threshold needs to be defined that reflects a notion of a critical size, below which the population is very unlikely to avoid extinction. We chose to illustrate the effects of using thresholds of 10, 50, and 100 individuals. Population viability analysis, then, is concerned with the cases when N [k] t drops below a critical threshold, often called quasiextinction. Statistical methods similar to those presented here can be found in Holmes, Sabo, Viscido, and Fagan (2007). In the simulated predictions of population values, the first year that N [k] t falls below such a threshold was recorded, for various thresholds. For each k, an indicator that the population falls below ξ by year T is given by: where I(·) is the indicator function, which is 1 when its argument is true, and 0 otherwise. Across all simulated populations, the proportion of times that a threshold is crossed, accumulated by year, provides an estimate of the probability of crossing the threshold in the years following 2013. That is, an estimate of crossing some threshold ξ by time T is: q(T, c) = 1 10, 000 10,000 k=1 f (T, k, ξ).
Note that other uses of population viability analysis (PVA) include estimating time to extinction (Taylor, 1995), which we do not consider here.

RESULTS
We first show posterior distributions for the demographic parameters and covariate effects to confirm that the values are reasonable with respect to our understanding of harbor seal biology. We then turn our attention to derived quantities from the model, such as abundance, trend, and population viability analysis.

Distributions of Demographic Parameters
The main drivers of trend and abundance are the demographic parameters, δ t , c t , φ t , and κ t (see Table III). Smoothed posterior densities for each of these quantities are given in Fig. 1. For each parameter, all values across the t years were pooled, so the subscript is dropped when discussing these results. The first column in Fig. 1 contains the results under the base-priors scenario (Table IV). Note that there is very little change from the prior distributions to the posteriors, indicating a lack of information in the data compared to the priors. The average value of the posterior for nonpup annual survival rate (δ; first row, first column) was 0.893, for fecundity (φ; second row, first column of Fig. 1) was 0.232, and for pup survival to nonpup status, κ, it was 0.778. The posterior distribution of the density-dependent scaling factor, c (third row, first column of Fig. 1), depends on the posterior distribution of ρ and various values of N t . We show the posterior distribution of c for N t = 300, 600, 1200, which illustrates how fecundity cφ declines with increasing values of N t due to density dependence.
The prior distributions were somewhat broader for the wide-priors scenario, and the resulting posteriors diverged more from the priors, shown in the second column of Fig. 1. Here, the average of the posterior values for δ was 0.861, for φ was 0.268, and for κ was 0.782. The posterior for c for the wide-priors The posterior densities for density dependence c are also shown for various population sizes N. The first column of plots uses parameters from the base-priors scenario, and the second and third columns of plots use parameters from wide-priors and uniform-priors, respectively. scenario was very similar to that for the base-priors scenario.
The prior distributions were very wide for the uniform-priors scenario, and the resulting posteriors diverged substantially from the priors, as seen from the third column of Fig. 1. The average of the posterior values for δ was 0.851, for φ was 0.313, and for κ was 0.782. Once again, the posterior distribution of c was similar to both base-priors and wide-priors scenarios. In all three cases, most of the posterior distribution for c is very near 1.0 for population sizes less than 1,000.

Effects of Date and Time of Day
As expected, there was a strong effect of day of year on pup counts (Fig. 2a). The mean day of year for pup counts was 206 (26 July), and the estimated date of peak pup numbers was 20 July. Counts as early as day 110 (21 April), and as late as day 301 (29 October), are expected to have around 30 fewer pups than at the peak. The uncertainty in these temporal patterns is reflected in the variation in the fitted (gray) curves of Fig. 2a. Note that the simple quadratic model may not be fitting very well near the ends of the data because we would not expect to see any pups at those dates, nor at any dates more than about three weeks from the peak. In early April, no pups are born yet, and by late October, they are large enough to be difficult to distinguish from older seals, so these changes reflect both the annual pupping cycle and the decline in observability of older pups. There was a minor peak in the effect of time of day on pup counts at about 15:00, which was one hour later than the mean time of the surveys (14:00; Fig. 2b). However, substantial variation among MCMC samples indicated there is little evidence that more pups are expected to be counted at any particular time of day.
The peak in expected counts of nonpups occurred later than for pups, on 10 August (Fig. 2c). In comparison to the pups, counts of nonpups decreased more rapidly away from peak counts. Note that these are expected changes in absolute number, and there are many more nonpups than pups in the population, which can allow for greater increases and declines. The peak count of nonpups was at 13:00 (Fig. 2d), 1 hour earlier than the mean time of surveys. However, there was considerable variation in the curves and, similar to pup counts, there was only weak evidence for a peak time of day in nonpup counts.

Abundance
For k = 1, . . . , 10, 000 MCMC samples, the abundance estimates N [k] t for years t = 1, . . . , 30 are shown in Fig. 3a, with priors on the demographic parameters from the base-priors scenario (Table IV). For each t, the k different averaged values are denoted asN t . Notice the effect of the constraint that N [k] t must be greater than the maximum observed count in any year, inflated by a factor of 1/0.9. The highest counts in 1998 and 2008 were largely responsible for determining how many seals were missed due to being in the water, from average counts (note that the year ticks are set to 1 January, but surveys generally occurred in summer, and hence the offset in time). The thin solid line that is fitted to the observed data is essentially ψN t , with some correction for covariate effects, reflecting the expected number of seals that are "observable" (i.e., at haul-out sites). The posterior distribution of ψ is shown in Fig. 3b. Notice that although most of the mass of the prior was between 0.2 and 0.8, the posterior distribution of ψ had most of its mass between 0.35 and 0.5, with a mean of 0.417, indicating that on average slightly more than half of the seals were missed. Fig. 3c shows the same information as in Fig. 3a, but using the wide-priors on the demographic parameters. Even with the same prior on ψ, the greater uncertainty in the prior demographic parameters led to more uncertainty in the abundance, as might be expected (Fig. 3c). Also note that this greater uncertainty allowed the abundance estimates to "wrap" around, to a greater extent, the highest observed values in 1998 and 2008. For the wide-priors scenario, the posterior distribution of ψ is shown in Fig. 3d, and it appears to be quite similar to that from the base-priors scenario (Fig. 3b), with a mean of 0.423. Fig. 3e shows the same information as given in Figs. 3a and c, but using the uniform-priors scenario for the demographic parameters. This scenario had the greatest uncertainty in the prior demographic parameters and led to the greatest uncertainty in the abundance, especially in those years with few or no surveys. For the uniform-priors scenario, the posterior distribution of ψ (Fig. 3f) again appears to be quite similar to those from the other scenarios, with a mean of 0.434.

Trend
For each MCMC sample k, the values δ can be used to form the matrix M t (Equation (4)) for each time period t, yielding 290,000 such matrices. The posterior distributions of the leading eigenvalues from M t are given in the first column of Fig. 4 (matching the situations in Fig. 3). The leading eigenvalue is the annual factor of increase at the stable age distribution for each Leslie matrix. Despite the fact that the wide-priors scenario had parameters of M t leading to zero growth, and the uniform-priors scenario had parameters of M t leading to negative growth, the posteriors in all three scenarios had means showing positive growth: basepriors scenario = 1.062, wide-priors scenario = 1.055, and uniform-priors scenario = 1.074. The base-priors scenario had the least variability, and uniform-priors scenario had the most.
Linear regression estimates of trend were described in Section 3.4. In particular, the posterior distribution of the most recent 15 years, using the sample of τ [k] 15 , k = 1, . . . , 10, 000, is given in the second column of Fig. 4. For the base-priors scenario, there was a posterior probability of 0.5075 of a decreasing trend, indicating little evidence that the population was decreasing over this time period. For the widepriors scenario, there was a posterior probability of 0.8266, and for uniform-priors, there was a posterior probability of 0.7867, of decreasing trends, neither of which provide strong evidence for a decline. Nevertheless, these values highlight the difference between the annual factors of population increase computed from the Leslie matrix eigenvalues, which were positive and uncorrected for harvest, and linear trends on the last 15 years of abundance, which were negative and included harvest.
The posterior distributions of trend for the most recent 10 years, using the sample of τ   Fig. 4. Trend summaries from the posterior distribution, computed from the 10,000 MCMC samples. The first column is the posterior density of the first eigenvalue of the Leslie matrix, computed on the samples of δ, φ, and κ and used in Equation (4). The first row is from the base-priors scenario, and rows two and three are from the wide-priors and uniform-priors, respectively. The second through fourth columns are least-squares linear trends (individuals/year) through the abundance estimates (all of the gray lines in Fig. 3), for the last 15, 10, and 5 years of abundance estimates, respectively. All posterior densities were smoothed with a kernel density estimator. 1, . . . , 10, 000, are given in the last column of Fig. 4. For the base-priors scenario, there was a 0.4961 posterior probability for wide-priors a 0.584 posterior probability, and for uniform-priors a 0.1383 posterior probability of a decreasing trend; the latter indicates modest evidence of an increasing trend in the last five years.

Population Viability Analysis
When projecting into the future, as described in Section 3.5, the survival and fecundity rates estimated under the base-priors scenario produced little variability (Fig. 5a). Note that for the first 30 years, the lines N [k] t ; t = 1, . . . , 30 for k = 1, . . . , 10, 000, are the same as in Fig. 3a, but plotted on the log scale. The population predictions become more variable when projecting into the future, N [k, j] t ; t = 31, . . . , 130; k = 1, . . . , 10, 000. Under this scenario, there was almost no chance that the population would drop below any of the thresholds,q(T, c) (Equation (5)), of 10, 50, or 100 individuals (Fig. 5b). Parameters from the wide-priors scenario produced more variabilty (Fig. 5c). Under this scenario, some of the future population trajectories crossed the various thresholds. The cumulative probability that the population would drop below 100 by year 100 is about 9%, and there is a 1-2% chance of dropping below 50 (Fig. 5d). Under the more variable uniformpriors scenario (Fig. 5e), the cumulative probability that the population would drop below 100 is about 5% within 30 years, and around 15% by year 100 (Fig. 5f). There is also about a 4% chance that the population could drop below 50 by year 100, but still very little chance that the population would drop below 10. Hastings et al. (2012) compared their survival and productivity estimates to other literature on harbor seals in Alaska and elsewhere. Using their estimates in our base-priors scenario, along with the precision of those estimates, the data did little to move the posterior distributions from the prior distributions (Fig. 1, top row). However, with our widepriors and uniform-priors, there was some evidence that both adult survival, δ, and fecundity, φ, were higher than the prior values (Fig. 1, bottom two  rows). This may not be surprising based on traditional knowledge of hunters around Iliamna Lake, who assert that seals in the lake are larger and fatter than their saltwater counterparts (Burns et al., 2016). Indeed, in all three scenarios, the posterior distribution of the annual factor of increase (i.e., population growth rate) implied by the Leslie matrix parameters was greater than the prior; (basepriors = 1.062, wide-priors = 1.055, uniform-priors = 1.074). These growth rates make sense in light of the harvest information (note that the harvest data in the model helped determine them). On average, Table II indicates about 20 seals are harvested per year. Based on a population of approximately 400, that indicates about 5% harvest rate, which essentially balances the production rate, consistent with the apparent near-equilibrium in numbers over the past several decades.

Effects of Date and Time of Day
The influence of time of day and day of year on harbor seal counts in Iliamna Lake mostly matched our expectations from patterns described for marine harbor seals. Despite considerable uncertainty, there was a weak peak just after midday for both pups and nonpups around 1400 local time (Figs. 2b and 2d), which is very near solar midday for this location and exactly the peak time expected from marine harbor seal studies (e.g., Simpkins et al., 2003). There was a much stronger relationship to day of year. The August peak in Fig. 2c generally matched our expectations, based on the greater amount of time that harbor seals tend to spend ashore when they are most actively molting their pelage. Pupping generally peaks around mid-June in Alaska (Mathews & Pendleton, 2006), but aerial counts are confounded by changes in the absolute numbers of pups, the time that they spend in the water, and the difficulty in distinguishing them from nonpups as they grow during the summer, which generally explains the broad peak in Fig. 2a. A third covariate often modeled for seal counts in marine populations is the effect of tide, which is not present in Iliamna Lake.

Abundance and Trends
The results in Fig. 3 suggest a minimum of around 300 seals occurred in the mid-late 1980s, but there was a large amount of uncertainty, with only a few, relatively low counts obtained during that period (e.g., 77 seals in 1984). It is likely that the population grew through the 1990s and has stabilized at around 400 individuals. Using the wide-priors scenario as shown in Figs. 3c and 3d, the most recent abundance estimate, in 2013, was 398 with the 95% credible interval ranging from 326 to 485. The peak abundance estimate in the recent years of increased monitoring occurred in 2008 and was estimated at 460, with the 95% credible interval ranging from 412 to 536. This corresponds reasonably well with the 389 estimate by ABR, Environmental Research and Service (2011), which would be expected to be lower because it did not include a correction (ψ) for the proportion of seals in the water.

Population Viability
The probability of quasi-extinction in the next 100 years, using 50 as a threshold, was estimated to be around 1% under wide-priors and 4% for uniform-priors. We generally disregard the basepriors scenario, as the variances in those priors were developed based on sampling error alone, and not on annual stochastic variation in demographic parameters. In the Section 1, we reported that the concept of extinction risk due to demographic stochasticity is often separated from environmental stochasticity. However, our analysis includes both because the variation in estimated demographic parameters includes normal environmental variation. That leaves genetic and catastrophic environmental effects as unmodeled. A quasi-extinction threshold of 50 may be large enough to avoid genetic effects somewhat, but we acknowledge that populations of long-lived animals that remain at such low levels for very long will likely incur genetic challenges (Frankham, 1996). Finally, it is important to consider catastrophic events. For example, catastrophic population declines have occurred in otariids due to El Niño or unknown causes (Gerber & Hilborn, 2001). Our results are based on a span of 30 years, which is a relatively long time for monitoring a marine mammal population. However, a population's risk may depend on rare events that occur less frequently than any observed during the sampling time frame. Based on accounts by local residents, Burns (1978) and Mathisen and Kline (1992) reported that the numbers of seals in Iliamna Lake dropped after extremely severe winters in 1970-1971 and 1971-1972 to approximately 40-50 individuals. Although the evidence is anecdotal, harbor seals in Iliamna Lake may be subject to such catastrophic events, which increases their risk beyond those that are modeled in this article. It is difficult to quantify such catastrophic risk (but see Ward, Hilborn, Towell, & Gerber, 2007, for a quantitative approach when there are better data), so any management decisions should consider the risk presented in this article as the minimum risk, and our model could be used to evaluate the effects of hypothetical scenarios for the likelihood of various kinds of catastrophes.
For our data, changes in assumptions (priors) sometimes led to relatively large changes in inference. For example, the population parameters contained in the scenarios as shown in Fig. 1 led to the different trends (Fig. 4) and population trajectories (Fig. 5), with differing chances of crossing a quasi-extinction threshold. As such, these scenarios should be considered a compromise between projections and predictions, in the sense of Keyfitz (1972) and Caswell (2001, p. 29). That is, a prediction is an inference on what will happen, and a projection on what could happen given a certain set of assumptions. Here, the data modify the starting assumptions somewhat, but not enough to overwhelm those assumptions. We have presented several reasonable scenarios through the specification of prior distributions. The results largely match the hypotheses and expectations that we outlined in Section 1, but the analysis quantified their values and illustrated the uncertainty stemming from the relatively sparse data available for this poorly documented seal population. This should form a basis for current management decisions, designing future sampling efforts, and a comparison point for future analyses.

ACKNOWLEDGMENTS
This article is dedicated to the memory of Professor Daniel Goodman of Montana State University, who advocated strongly for the use of demography, Bayesian methods, and quantitative risk assessment in conservation decisions. Aerial surveys were authorized under a Marine Mammal Protection Act General Authorization (LOC No. 14590). We thank Paul Conn, Donna Hauser, Kelly Hastings, and Devin Johnson for reviews of early drafts, and two anonymous reviewers for helpful comments on the article. The findings and conclusions in the article are those of the authors and do not necessarily represent the views of the reviewers nor the National Marine Fisheries Service, NOAA. Any use of trade, product, or firm names does not imply an endorsement by the U.S. government.

DATA AND CODE ACCESSIBILITY
An R (R Core Team, 2017) package called IliamnaSealsTrendPVA was created that contains all data, code, and analyses. This document was created using knitr (Xie, 2014(Xie, , 2015(Xie, , 2016, and the manuscript combining L A T E X and R code is also included in the package. The package can be downloaded at https://github.com/jayverhoef/ IliamnaSealsTrendPVA.git , with instructions for installing the package.

APPENDIX: PRIORS FOR DEMOGRAPHIC PARAMETERS
To develop informative priors, we used the detailed information of sex-and age-specific survival obtained from Hastings et al. (2012), and the information on female productivity contained in Pitcher and Calkins (1979), which are given in Table AI. The average female fertility rate for all females of age ≥ 3 was obtained from Pitcher and Calkins (1979) in a study of harbor seals in the Gulf of Alaska, the nearest available source for such data. They collected 165 females aged 3+ at various times in the annual cycle between ovulation (early July) and the next pupping-weaning-ovulation season. Among those 165 individuals, 143 had ovulated. Pitcher and Calkins (1979) found that implantation of the embryo took place around early October. There were 132 females in their sample that were collected between the time of implantation and the birth season; 14 of those either failed to become fertilized or lost their embryo or fetus. We estimated the raw birth rate, then, as the ovulation rate (143/ 165 = 0.867) minus the pregnancy failure rate (14/ 132 = 0.106), which equals 0.761. Because the females collected for evidence of ovulation and failure of pregnancy were collected at various (unspecified) times in the annual cycle, this raw birth rate accounts incompletely for females that would not have been available to be sampled because of mortality during the interval (i.e., the data are censored). To adjust for this, we applied one-half of the annual mortality rate for the 3+ age class, (1 − 0.929)/2 = 0.036, to the raw birth rate to yield a fecundity rate appropriate for use in a Leslie-type population projection matrix, 0.761 * (1 − 0.036) = 0.733. The variance was developed assuming independent binomial proportions for raw birth rate and pregnancy failure, and used Goodman's result on direct products (Goodman, 1960), assuming again that mortality rate was independent of birth rates and pregnancy failures.
We began with a full sex-and age-specific Leslie matrix model for harbor seals, then adjusted it for the timing of our surveys, and finally collapsed it to the two classes for which we had information: pups and nonpups. The full Leslie matrix model, developed from data in Hastings et al. (2012) and Pitcher and Calkins (1979), is n i+1 = Mn i , where M is given by: and n i is: ⎛ where n is the number of animals in various sex/age classes, "M" indicates male, "F" indicates female, "0" indicates age at birth, "1" indicates first birthday, "2" indicates second birthday, and "3+" indicates birthdays 3 and higher. All of the survival values in this model were taken from Hastings et al. (2012,  Table 2, part a). We chose these values because they were from a detailed study of well-marked seals, temporally and geographically close to our study in Iliamna Lake. Note that, from genetic considerations, we assumed that the productivity rate was split equally between production of male and female pups, resulting in an expected 0.367 probability of a male pup, and the same probability of a female pup, per female aged 3+. Our surveys generally occurred around early August. If not, counts were modeled as being adjusted by a covariate for date as if they were counted in early August. Hence, there is a misalignment with the demographic parameters in Equation (A1) and our observations because seal pups are born prior to most of the aerial surveys. We corrected for that as follows: Hastings et al. (2012) gave two survival rates for pups; from the time they are very small and from the time that they are large. We take these to be equivalent to the time of birth and the time of our surveys. Harbor seal pups grow rapidly, and sometime in August it becomes unreliable to distinguish pups from nonpups in aerial surveys, so the survival of large pups aligns very well with our survey period. For male pups, survival from birth to yearling is reported as 0.405, but survival from a large pup to yearling is 0.717, and hence we take the estimate of survival from small to large pup as 0.405/0.717 = 0.565. Likewise, the female survival from small to large pup is 0.549/0.82 = 0.67. We adjusted the productivity values then not for birth, but for the production of pups that survive to the time of surveys. Hence, females aged 3+ have a 0.367*0.565 = 0.207 probability of producing a male pup that survives to the survey period, and females aged 3+ have a 0.367*0.67 = 0.245 probability of producing a female pup that survives to the survey period. This requires modifying Equation (A1) so that M is: There are several characteristics of this population model. The first (dominant) eigenvalue of the matrix, for both matrices A1 and A3, is 1.0566, indicating an annual growth rate of 5.66%. The first eigenvector associated with the dominant eigenvalue for matrix A3, when scaled so that all elements sum to one, is, which is the stable age distribution, shown with i = ∞ for the vector, n, of the number of animals in each sex/age class. Note that we have changed the notation on n to reflect the shift away from birthdays to the time of surveys, indicated by the subscript s on the pup class, and the survey time between birthdays using a dash for other age classes.
Our objective was to collapse the Leslie matrix in A3 into two classes. That is, we seek a simplified model: where φ is the probability, among all nonpups alive at one annual survey, of producing a pup that survives to the time of next year's surveys, κ is the survival of those pups from the time of aerial surveys to the nonpup class, and δ is the survival of the nonpups to stay within that class. Assuming the population is at a stable age distribution, from Equations (A3) and (A4), we used the following weighted averages:δ where S δ = v 2 M 3,2 + v 3 M 4,3 + v 4 M 4,4 + v 6 M 7,6 + v 7 M 8,7 + v 8 M 8,8 and D δ = v 2 + v 3 + v 4 + v 6 + v 7 + v 8 . This yields the matrix: Notice that, because we used weighted averages of the stable distribution of matrix A3, the eigenvalue of matrix A6 is 1.0566, which is exactly the same as that for matrix A3, and the stable proportion of pups from matrix A6 is 0.177, which is exactly the sum of the male and female stable proportions from A4. We also developed the variances for the parameters estimated in A6. To do so, we used the delta method (Dorfman, 1938;Ver Hoef, 2012) twice. First, consider the following Taylor series approximation for a vector of functions of random variables: where , and x = (X 1 , . . . , X m ) for random variables X i , and μ = E(x) where the ith column of μ is μ i . Then if is the covariance matrix of x, then: cov(f(x)) ≈ D D .
Notice that there are eight unique entries in A3, but some of them are functions of the quantities in Table AI. We assume that all variables in Table AI are independent, so that the covariance matrix is diagonal. However, using Equation (A8) for the functions of the estimates in Table AI contained in A3 results in the covariance matrix , which is:  Next, notice that δ, φ, and κ are functions of the entries in M, but that function is complicated and nonlinear because some entries occur multiple times in M, and because of the scaled eigenvector v. In fact, we can write them showing the dependence on m as δ(m), φ(m), and κ(m), where the exact dependencies are given in Equation (A5). We need the derivatives of δ(m), φ(m), and κ(m) with respect to every element of m. If we write f i (m j ), where m j is the jth element of m and f i (m) is the ith function, which is one of δ(m), φ(m), or κ(m), then the i, jth component of the matrix D will be a numeric approximation to the derivative by letting D i, j = ( f i (m j + h) − f i (m j ))/ h for some small value of h. We used 1 × 10 −5 for both positive and negative values of h for all i and j, and tried 1 × 10 −4 and 1 × 10 −6 but found there was very little change in D for smaller values of h. Then, using Equation (A8) again, with as given in A9, we obtained the following covariance matrix among the estimates of δ, φ, and κ, Finally, the prior distributions for these probabilities should be constrained between 0 and 1. We did this by matching moments with a beta distribution. A beta distribution, given as var(X) = ab (a + b + 1)(a + b) 2 .
For the estimated value of δ in Equation (A5), and its variance in A10, we obtained a = 272.63 and b = 33.52. For developing priors, we did not use the estimated covariances amongδ,φ, andκ. Using this basic framework, we developed three prior scenarios. "Base-priors" used the parameters as estimated above, and moment matching to a beta distribution. However, notice that the variances in basepriors are estimates for parameters that are assumed fixed, and from a nearby saltwater location. Our model allows for annual variation in the demographic parameters, and there is little reason to assume that the variance of estimation is related to annual variation. Therefore, it would be prudent to allow the variation to be wider. We decided to triple the standard errors (nine times the variance of each parameter). In addition, the base-priors biased the population toward positive growth because the dominant eigenvalue is 1.0566. In "wide-priors," we wanted to be neutral on population growth. Therefore, we wished to scale matrix A6 so that the eigenvalue was 1. That is, we wanted to find some α such that: Notice that when α = 1, λ is the eigenvalue of the matrix, so α = 1/λ yields an eigenvalue of 1. For the values in A6, it is easily determined that α = 0.946 creates the desired matrix. These values, and the variance for wide-priors, are given in Table IV , along with the derived a and b parameters for the beta distribution. As a third scenario, "uniform-priors," we used flat uniform priors over the interval 0.5-1 for δ and κ, and over the interval 0-0.5 for φ. Note that using the mean for each of these parameters yields a dominant eigenvalue of 0.948, biasing the population toward negative growth, but with the widest prior variances. Our three scenarios are summarized in Table IV.