Comparing estimates of population change from occupancy and mark-recapture models for a territorial species

Monitoring studies often use marked animals to estimate population abundance at small spatial scales. However, at smaller scales, occupancy sampling, which uses detection/nondetection data, may be useful where sites are approximately territories, and occupancy dynamics should be strongly correlated with population dynamics. Occupancy monitoring has advantages in that it is less expensive and invasive, and marked animals are not needed. Here, we used empirical data to determine whether and when change in occupancy is a good proxy for population change for a territorial species. As part of this overall goal, we also compared maximum-likelihood estimates using a model-averaging approach with a Bayesian MCMC approach. We used field data collected from 1993 to 2013 on three study areas for California spotted owls (Strix occidentalis occidentalis), a territorial species. Although correlations for trajectories of realized population change (Δt) between territory occupancy and Pradel models were moderate to high for Bayesian MCMC-based estimates and high for model-averaged estimates, magnitudes of the trajectories were different with the Pradel model reporting greater magnitudes of change. For the two areas showing a decline, Δt for the Pradel model was approximately 20–30% lower than for the occupancy model, and 25% higher in the area showing an increase. These differences can arise because the occupancy model is less sensitive, in that if two owls share a territory, the loss of one may be reflected in survival and, consequently in Δt by the Pradel model, but because the territory remains occupied it is not reflected by the occupancy model. Bayesian MCMC-based and model-averaged estimates of Δt were in close agreement in pattern (correlation ≥0.74) and magnitude (relative differences of last Δt were ≤5%) for both occupancy and mark– resight models. Results from the Pradel model may lead to conservation actions necessary to avoid high extinction or extirpation risk for small populations, while results from the territory occupancy model may result in status quo management. We found both Bayesian MCMC-based and model-averaged estimates of Δt robust approaches to evaluate population trends. However, we recommend the Bayesian MCMC approach for estimating risk (e.g., probability of declines) for retrospective analyses.


INTRODUCTION
Ecological and conservation monitoring methods tend to diverge depending on spatial scale. Vital rate and abundance monitoring are often reserved for populations or smaller areas because it is logistically difficult and prohibitively expensive to estimate changes in these parameters across large areas over time. When monitoring programs are large scale, such as at bioregional scales, detailed demographic information at a few locations may be exchanged for less detailed information about a larger area. At broad spatial scales, species distributions and their dynamics are often the metrics of interest, and these metrics can be described by the proportion of sites occupied (ψ) over time. Occupancy monitoring, which uses detection/nondetection data, can facilitate time-sensitive assessment of population status across large landscapes (Zuckerberg et al. 2009, Roney et al. 2015. In a monitoring context, occupancy can be used as a surrogate for population abundance or mark-resight estimates of rate of population change (λ), both of which are more expensive and require more effort to reflect the state of a population (MacKenzie et al. 2006).
Population performance can be robustly depicted by realized population change (Δ t ), which is an encompassing metric of population trend over a period of time (Franklin et al. 2004, Conner et al. 2013. In essence, it is the ratio of population size at a subsequent time period (t) relative to the initial population size. For longterm monitoring data, Δ t is an insightful metric for portraying the cumulative dynamics of a population across medium and long timescales, such as 5-20 yr (Gerrodette and Rojas-Bracho 2011). For example, if Δ t = 0.78 for an 18-yr period, then the population declined by 22% over that period. This metric can be estimated from population abundance estimates (N) as Ñ t ∕Ñ 1 , or from occupancy estimates as ψ t ∕ψ 1 , or from the product of ̂λ t based on vital rates of marked animals .
While marked animals are needed to estimate population Δ t based on abundance or λ t from vital rates, they are not required to estimate Δ t for the occupancy model because sites are the sampling units. For territorial species, a territory may represent a single animal, a pair, or a group (e.g., pack) of animals. If occupancy sampling is used with territories equaling single sites in an occupancy model, then occupancy changes should be strongly correlated with population changes, such that territory occupancy Δ t should be strongly correlated with population Δ t (Tempel and Gutiérrez 2013). Many wildlife and conservation agencies are faced with present budget reductions as well as likely future reductions, which may result in decreased monitoring. Occupancy monitoring offers the potential of providing information about population dynamics at reduced cost compared to more intensive mark-recapture approaches (MacKenzie et al. 2006(MacKenzie et al. , 2010(MacKenzie et al. , 2011. When estimating Δ t , Bayesian methods offer an advantage in that the posterior distribution of Δ t provides a robust method for detecting and describing declines (Wade 2000). Indeed, probabilities of declines can be addressed directly, providing information beyond the acceptance or rejection of a null hypothesis (Eguchi andGerrodette 2009, Gerrodette 2011). For example, for a population of California spotted owls (CSO; Strix occidentalis occidentalis), monitored from 1990 to 2011, Δ t = 0.78 with a CI of 0.54-1.08 and the statistical null hypothesis of no decline was not rejected (Conner et al. 2013). However, based on the posterior distribution of Δ t , the probability that the population declined was high (0.94) and the probability that the population declined ≥10% was substantial (0.80). The problem of having a large probability of decline, even when the null hypothesis of no decline is not rejected, is especially important for small populations. Smaller populations can drop, undetected by a null hypothesis approach, to low levels where the probability of extinction increases to unacceptably high levels due to demographic variability and stochastic environmental events (White 2000b, Lande 2001, Morris and Doak 2002. Thus, a Bayesian approach, if it performs well, is desirable for estimating probabilities of declines at specific levels of concern (e.g., ≥10%, ≥30%, and ≥50%).
Although correlations between occupancy and relative abundance have been found (e.g., Estrada and Arroyo 2012, Oliver et al. 2012, Roney et al. 2015, we are only aware of one comparison between trends in annual occupancy and population trends for the same study population (Tempel and Gutiérrez 2013). This v www.esajournals.org CONNER ET AL. comparison was specific; it examined the relationship between model-averaged values of Δ t from occupancy and the Pradel model for a declining population, where trajectories of occupancy Δ t and population Δ t were found to closely match (Tempel and Gutiérrez 2013). Here, we extend this work to three populations with three different trajectories: declining, relatively stationary, and increasing (Conner et al. 2013). In addition, because of the advantage of a Bayesian MCMC approach for Δ t , we wished to compare it to the model-averaging approach taken by Tempel and Gutiérrez (2013).
Our overarching goal was to more fully examine whether and when occupancy works as a surrogate for mark-resight population monitoring for territorial species using Δ t . Our second goal was to evaluate different methods for estimating Δ t for the occupancy and Pradel models. Specifically, we wished to compare the value of estimates based on Bayesian MCMC and modelaveraging approaches for evaluating population trajectories of species with long-term monitoring data. Here, we evaluate the relationship between the occupancy and mark-resight estimates using field data collected from 1993 to 2013 on three study areas for CSO, a territorial species with high site fidelity (Blakesley et al. 2006, Seamans and Gutiérrez 2006, Seamans and Gutierrez 2007 as a case study. The CSO remains a focal species of conservation concern, and knowledge of population trends is an important component of assessing status and informing forest management planning efforts.

Study area
We used data for spotted owls on three study areas in the southern Cascade and Sierra Nevada Mountains, California (Fig. 1), collected from 1993 to 2013. Although the study areas were not randomly selected, they spanned the length of the contiguous CSO range in the Sierra Nevada and encompassed all habitat types known to be used by spotted owls in the Sierra Nevada. The Lassen (LAS) study area was in the southern Cascades, but it was included in the Sierra Nevada province by the U.S. Forest Service for management purposes (USDA Forest Service 2004). Most of the LAS and Sierra (SIE) study areas were located on public land managed by the U.S. Forest Service, whereas the Sequoia and Kings Canyon (SKC) study area was in two national parks. Franklin et al. (2004) and Blakesley et al. (2010) described in detail the LAS, SIE, and SKC study areas; however, SKC has changed from Blakesley et al. (2010). In 2006, the study area for SKC was reduced; all analyses and estimates were based on data from this reduced study area. Longterm mark-resight data were available for each study population. The number of marked birds varied each year, but generally ranged between 40 and 90, with 2-17 new birds marked each year.
LAS, SIE, and SKC had territorial survey areas, where survey effort varied among years, and a subset core study area, which was surveyed consistently for estimation of λ t based on mark-resight data. In these core areas, survey coverage was complete, including areas between known territories. For all empirical comparisons, v www.esajournals.org CONNER ET AL.
we used data only from the core areas to meet the assumptions of the Pradel and occupancy models.

Field data collection
Although mark-resight data were collected prior to 1993, survey coverage of SKC was not adequate to meet the Pradel model assumption of constant study area size prior to 1991. Following Franklin et al. (2004), Blakesley et al. (2010), and Conner et al. (2013), λ t for the first two annual intervals were discarded in part to avoid potential bias; thus, the initial year for the Pradel Δ t estimate for SKC was 1993. We used the same initial year for all study areas and the same time periods for occupancy and Pradel analyses, which were 1993for LAS and SIE and 1993 for SKC (the SKC study ended in 2012).
Field methods for spotted owl population dynamics studies have been well described (Forsman 1983, Franklin et al. 1996, Pradel 1996, Anthony et al. 2006. In particular, Franklin et al. (2004) detailed the field methods and protocols used to collect the data in our study. Although the studies were designed to collect data to estimate survival and reproduction of individuals and pairs, the structure of the surveys made it possible to create data sets that were appropriate for modeling site occupancy as well. Because field methods have been well described, we present only a brief summary of survey methods that are relevant to the repeated site visits for generating encounter histories for the occupancy model, as well as for generating encounter histories used in the Pradel model.
Surveys to find and locate CSOs were conducted each year from 1 April to 31 August on LAS and from 1 March to 30 September on SIE and SKC. Surveys were conducted by vocally imitating spotted owl vocalizations. Most surveys were conducted at night, while some surveys were conducted during daylight hours to find roosts or nests, search for mates and fledglings, identify banded owls, and capture unmarked owls. For Pradel model λ t to represent changes in the number of owls in a population of interest, the study area size and boundary must remain unchanged through time Nichols 2002, Franklin et al. 2004) unless the model is structured correctly to accommodate a study area expansion, as in Blakesley et al. (2010). Consequently, we used data only from the core areas to interpret λ t as the annual rate of population change.
Following Tempel and Gutiérrez (2013), occupancy territories were defined as areas where we recorded diurnal (starting at civil twilight in the morning and ending at civil twilight in the evening) observations of owls in at least 3 yr from 1993 to 2012. We used only diurnal observations because CSOs are central-place foragers (Orians and Pearson 1979), so that the roost/nest area may be considered as a rough approximation to the territory center. Also, nocturnal observations may have been owls from neighboring territories or nonresident owls traveling through an area. We used the 3-yr requirement to define occupancy territories in order to screen out areas where nonterritorial owls may have been observed. This was particularly important for SIE and SKC data where March and September diurnal observations may have had a greater probability of being nonresident owls than the April through August observations. We averaged the diurnal locations (from 1993 to 2012) of each territory to define occupancy territory centers. We used half the mean nearest neighbor distance among occupancy territory centers for the radius of a circle around each center to define the occupancy territory.
Data from the field surveys were summarized differently for occupancy and Pradel models.
For the Pradel models, we used data from all known resident banded owls within the core areas to create encounter histories. An owl was considered resident on a territory in a given year if it was a member of a pair, twice in a season with observations at least 7 d apart, or if it was identified by color bands as the previous year's resident, due to high site fidelity in spotted owls (Blakesley et al. 2006, Seamans and Gutiérrez 2006, Seamans and Gutierrez 2007. If a resident marked owl was captured or resighted during any of the surveys in a year (i.e., it was identified as alive and on the study area), a 1 was entered into the encounter history; otherwise, a 0 was entered. For the occupancy model, we constructed an encounter history for each territory. Following Tempel and Gutiérrez (2013), we divided each sampling season (i.e., primary sampling period), which was approximately the breeding season, into 10 (LAS) or 14 (SIE and SKC) semimonthly v www.esajournals.org CONNER ET AL. periods (i.e., secondary sampling periods). For each semimonthly period, we assigned a 0 if no owls were detected and a 1 if at least one owl was detected during any survey. However, a 2 was assigned if fledglings were detected so that we could examine how reproductive success may have influenced detectability. We included only surveys that contributed to ≥30 min of survey effort within the occupancy territory without the detection of an owl, or surveys in which an owl was detected within the occupancy territory. If an occupancy territory was not surveyed during a semimonthly period or received <30 min of survey effort without detection of an owl, we assigned a missing value. One study area, SKC, had a year of missed data collection (2005), which resulted in fewer estimates of λ t and Δ t for both the Pradel and occupancy models, as described below.

Occupancy modeling
We used robust-design occupancy models to estimate annual territory occupancy dynamics (MacKenzie et al. 2003(MacKenzie et al. , 2006. We used ψ 1 γ ε p structure of the robust occupancy model, where ψ 1 = proportion of sites occupied during the first time period, γ = probability a site unoccupied in year t is occupied in the next year (t + 1; local colonization probability), ε = probability a site occupied in year t is unoccupied the next year (t + 1; local extinction probability), p = probability of detection at an occupied site, with year = primary sampling occasion (i.e., year), and visit = secondary sampling occasion. Estimated annual proportion of sites occupied (ψ t ), which we refer to as occupancy, is estimated from annual estimates of ε and γ (MacKenzie et al. 2006). We used Program MARK (White and Burnham 1999) for all analyses. Because data were not collected on SKC in 2005, ψ t could not be estimated for 2005. Consequently, Δ t = 2005 = ψ 2005 /ψ 1993 also could not be estimated.
Our model development procedure was a sequential process. We followed a common process of first modeling p while leaving the structure of other model parameters general (i.e., ψ 1 ε t γ t , where t = year modeled as different in each year; Nichols et al. 1997). We modeled p first with a general model that allowed p to vary categorically among all year by visit (j) combinations (p t × j ). We then constrained within-year variation following the strategy of Tempel and Gutiérrez (2013); that is, we used month as a categorical covariate to represent different states of the owl's breeding cycle, and fledging season (0 = prefledging, 1 March-31 May; 1 = postfledging, 1 June-31 September) to represent behaviors that may influence detection probability that are associated with the reproductive cycle. We also replaced visit with time trends (linear [T], loglinear [lnT], quadratic [TT], constant [.]) based on previous findings for detection probability (Seamans and Gutierrez 2007). Essentially, recapture models with T, lnT, and TT represent our hypothesis that observer proficiency increased as the studies progressed due to experience finding owls (T structure), but may level as there is a upper limit to detection (lnT and TT structures). Finally, we used reproductive status (0 = nonreproductive, 1 = reproductive) and initial detection (0 = before and up to first detection, 1 = after first detection) for each site (like individual covariates) as categorical variables to replace visit; see Tempel and Gutiérrez (2013) for details.
Using the top model for within-year p, we then constructed less parameterized models for ε, γ, and p in which year was modeled with time trends (linear [T], log-linear [lnT], quadratic [TT], constant [.]). We considered all possible combinations of temporal trends for ε, γ, and p. We did not use additional covariates (e.g., habitat variables within a territory) because our goal was to compare annual population dynamics using occupancy and mark-resight (Pradel's temporal symmetry) models.
We used ψ t to calculate occupancy Δ t as: Thus, overall change, which is estimated on the last year (k) of the study period, was estimated as: Δ k = ψ k ∕ψ 1 .
For the maximum-likelihood estimates from the model-averaging approach, we estimated the variance of Δ k using the delta method (Seber 1982). For the Bayesian MCMC approach, estimates and their distributions were from occupancy model ε t γ t p(best structure). We used MCMC sampling implemented in Program v www.esajournals.org CONNER ET AL.
MARK to estimate posterior distributions of ψ 1 , ε t , and γ t , from which we calculated the posterior distribution of ψ t MCMC . The posterior distributions of ψ t MCMC were used to estimate the posterior distribution of Δ k (Δ k MCMC ) as well as median Δ k and the 2.5th and 97.5th percentiles.
For all MCMC simulations, we used 4000 tuning samples, 1000 burn-in samples, and 20,000 realizations after thinning. We used vague priors for all parameters included in the model. Because all model estimates were logit-transformed parameters, we used a normal prior distribution with mean of 0 and a standard deviation of 1.75, which is a vague prior when back transformed to the real scale (2.5th and 97.5th percentiles of approximately 0.02 and 0.98, with a uniform distribution between those percentiles). We determined whether the Markov chains converged using the Gelman-Rubin statistic, R-hat (Gelman et al. 2004). For each parameter, we used 10 chains of 1000 each and used a threshold of Rhat < 1.1 to indicate adequate sampling of the posterior distribution of ψ 1 , ε t , and γ t . Before beginning the MCMC simulations, we used a routine in Program MARK's MCMC sampler to determine the number to thin by averaging the maximum number to thin for each parameter in order to achieve first-order Markovian independence. We estimated the median and 2.5th and 97.5th percentiles for the distribution of Δ k MCMC for each study area.

Mark-resight modeling
We used each owl's mark-resight encounter history to estimate the annual rate of population change in territorial owls (λ t ) using Pradel's temporal symmetry model (Pradel 1996, Franklin et al. 2004 in Program MARK (White and Burnham 1999). We used the (ϕ f p) structure of the Pradel model, where ϕ is apparent survival (probability that an owl alive in year t survived to the next year [t + 1] and remained on the study area [i.e., available for recapture/resight]), f is recruitment (number of territorial owls at t + 1 divided by number of territorial owls at t) and p is the recapture (by capture or resight) probability. For this form, λ is a derived estimate, calculated as λ = ϕ + f. For each study area, we estimated the overdispersion parameter (ĉ) using the median ĉ procedure in Program MARK under the Cormack- where s is an owl's sex. When ĉ was >1, we used ĉ to inflate variances of parameter estimates (Burnham and Anderson 2002).
Both LAS and SIE had expansion areas, areas in which surveying began several years after the start of the study. If these were not accounted for, new owls found in these areas would enter the Pradel model as new recruits, and result in a positive bias in estimates of λ t . Therefore, we accounted for these areas by grouping them separately and excluding estimates of λ t from the year of and year following the expansion, and estimates of p for the year of expansion (see Appendix G: Blakesley et al. 2010 for a detailed explanation of the procedure). In addition, as mentioned above, SKC had a year of missed data collection (2005). We used an unequal time interval to account for this in Program MARK, which resulted in average survival and recruitment over the 2-yr interval.
Similar to the occupancy modeling, we used a sequential approach for model development. We left the structure of ϕ and f general (i.e., ϕ t f t ) and first modeled p. Note that for the Pradel model, p is the probability of capturing an unmarked individual or of resighting or recapturing a marked individual that was still alive during a given year, where p in the occupancy modeling is the probability of detecting an owl(s) at an occupied territory during a survey. We did not include sex in models of ϕ and f to parallel the occupancy modeling, in which we did not consider the sex of a territory holder. In addition, previous analyses for CSO on the three study areas did not find sex effects on ϕ and f (Appendix H: Blakesley et al. 2010). We ran a series of models hypothesized to affect p. For p, we included sex (Blakesley et al. 2010), and modeled year with time trends (linear [s + T], log-linear [s + lnT], quadratic [s + TT], and constant [s]). The time trends represent our hypothesis that observer proficiency increased as the studies progressed due to experience finding owls (T structure), but may level as there is a upper limit to detection (lnT and TT structures). We also ran the time trend models without sex. For the next phase, we used the best structure on p and constructed models in which ϕ t and f t were reduced. Similar to the occupancy modeling, we also modeled year with time trends (linear [T)], log-linear [lnT], quadratic [TT], and constant [.]), using all possible combinations of temporal trends in ϕ t and f t . We used v www.esajournals.org CONNER ET AL. the derived parameter option in Program MARK to generate annual model-averaged estimates and variances of λ t (i.e., ϕ t + f t ).
We used estimates of λ t (λ t ) to estimate Δ t on each study area, which for the Pradel model is the proportional change in estimated population size, relative to the initial population size (Franklin et al. 2004 the last interval (k) for the study period, we calculated i Δ k = λ 1993 ×λ 1994 ×…×λ 2012 for LAS and SIE, and through ̂λ 2011 for SKC. For the modelaveraging approach, we estimated the variance of Δ k using the delta method (Seber 1982) and the variance-covariance of the model-averaged estimates. For the Bayesian MCMC approach, estimates and their distributions were from Pradel model ϕ t f t p(best structure). We used MCMC sampling implemented in Program MARK to estimate posterior distributions of ϕ t , and f t , from which we calculated the posterior distribution of λ t MCMC . The posterior distributions of λ t MCMC were used to estimate the posterior distribution of Δ k (Δ k MCMC ) as well as median Δ k and the 2.5 and 97.5 percentiles. We did not use the average estimate from the 2 yr (2004-2005 and 2005-2006) that arose from the year of missing data collection on SKC to avoid biasing the variance for estimates of Δ t and ultimately Δ k .
For all MCMC simulations, we used 4000 tuning samples, 1000 burn-in samples, and 20,000 realizations after thinning. We used vague priors for all parameters included in the model. For estimates of ϕ, which were logit-transformed parameters, we used a normal prior distribution with mean of 0 and a standard deviation of 1.75, as described for occupancy parameters above.
For estimates of f, which were log-transformed parameters, we used a normal prior distribution with mean of 0 and a standard deviation of 1.75, which is a vague prior when back transformed to the real scale (2.5th and 97.5th percentiles of approximately 0.04 and 30.44). We determined whether the Markov chains converged using the Gelman-Rubin statistic, R-hat (Gelman et al. 2004). For each parameter ϕ t , and f t , we used 10 chains of 1000 each and used a threshold of R-hat < 1.1 to indicate adequate sampling of the posterior distribution. Before beginning the MCMC simulations, we used a routine in Program MARK's MCMC sampler to determine the number to thin by averaging the maximum number to thin for each parameter in order to achieve firstorder Markovian independence. We estimated the median and 2.5th and 97.5th percentiles for the distribution of Δ k MCMC for each study area.

Occupancy
We located between 37 and 74 owl territories in the three study areas during 1993-2013 (Table 1). Of these, reproduction occurred on 84-93% at least once during the study, which indicates that we had identified biologically relevant territories. The top structure on p for all three study areas had month or visit as an additive effect with reproduction or initial detection effects (Table 2). Whether an owl was reproductive had a strong positive effect on detection probability for LAS and SKC, while whether an owl was initially detected had a strong positive effect on subsequent detection probability for SIE (Table 2).
There were different patterns in ε and γ among the study areas. LAS had an increasing ε and a relatively stationary γ, while SIE and SKC had v www.esajournals.org CONNER ET AL. Note: Only models with mass (w i ) > 0.001 are shown. † Key to model notation: K = no. of parameters; QAIC c = quasi Akaike's information criteria corrected for small sample size and lack of model fit; ΔQAIC c = difference between the model listed and the QAIC c of the best model; w i = model mass based on model QAIC c compared to all other model QAIC c values; T = linear time trend; lnT = logarithmic time trend; TT = quadratic time trend; . = constant time; yr = year; visit = semimonthly period during reproductive period, modeled as a categorical variable; month = month during reproductive period, modeled as a categorical variable; repro = categorical site covariate indicating whether a juvenile(s) were detected at the nest; initial = site covariate that separated recapture probability from initial capture probability; ε = local extinction probability; γ = local colonization probability; p = probability of detection.

Mark-resight
To model rate of population change, we used encounter histories for 331, 334, and 171 individual owls, for LAS, SIE, and SKC, respectively, with an average of 49-73 marked owls/ year located on the study areas (Table 1). Median ĉ ranged from 1.00 to 1.16, suggesting some overdispersion but no serious lack of fit (Table 3). We adjusted all variance estimates by site-specific estimates of c.
There were similar patterns in ϕ for all study areas; that is, ϕ was slightly increasing (Fig. 3). Patterns in f were different, with f slightly decreasing on LAS, remaining stable on SIE, and   Note: Models shown are those with w i > 0.001. † Key to model notation: K = no. of parameters; QAIC c = quasi Akaike's information criteria corrected for small sample size and lack of model fit; ΔQAIC c = difference between the model listed and the QAIC c of the best model; w i = model mass based on model QAIC c compared to all other model QAIC c values; s = sex; t = year effect as a categorical variable; T = linear time trend; lnT = logarithmic time trend; TT = quadratic time trend; . = constant time; ϕ = apparent survival; f = recruitment of territorial owls; p = probability of detection.
‡ The number of parameters was adjusted for study area expansion, which added three parameters. showing some annual variability on SKC (Fig. 3).
Top models for all three study areas had temporal trends for ϕ (Table 3). Top models for f were mixed; that is, there were either temporal trends or f was constant. Patterns in λ followed from ϕ and f; for LAS and SIE, λ slightly increased during the study period, while for SKC λ varied from year to year (Fig. 3).

Occupancy vs. mark-resight
For LAS and SIE, the trajectories for occupancy and Pradel Δ t showed declines, while for SKC the Pradel model showed an increase while occupancy showed relatively little change (Fig. 4). By the end of the study period, MCMC-based and model-averaged estimates of Δ t were similar (all ≤5% different, Fig. 5); for simplicity, we compare occupancy and Pradel MCMC-based estimates. The Pradel model showed more dramatic patterns of overall change than the occupancy model. The Δ k for the Pradel model at the end of the study period was 29% lower for LAS, 21% lower for SIE, and 26% higher for SKC relative to the occupancy model ( Fig. 4 and Table 4). The variance for the Pradel estimates of Δ t was greater than for the occupancy estimates; the 95% CI lengths of overall Δ t were 1.6-3.4× longer for the Pradel estimates.
Posterior distributions of Δ k were used to estimate probabilities of declines and increases in territory occupancy and population (Table 5). Different risks were revealed by the two models. For LAS and SIE, the biggest differences were for the ≥30% decline; while territory occupancy models showed low probabilities (0.08 and 0.01) for this larger decline, the Pradel model showed substantial-to-moderate probabilities (0.89 and 0.53; Table 5) for population change. However, the probability of any decline (≥0%) was similar between the models for LAS and SIE. In contrast to LAS and SIE, for SKC the probability of any decline of territory occupancy was much higher (0.56) compared to the Pradel model (0.13; Table 5). Note that for SKC, this comparison is not biologically meaningful because ψ 1 was approximately 1 (0.99); thus, the only real way its occupancy could move was down.

Bayesian MCMC-based vs. model-averaged estimates
The MCMC chains for all occupancy model parameters and areas converged (R-hat < 1.1).
Model-averaged and Bayesian MCMC estimates had similar trajectories; correlations between Δ t ranged from 0.77 to 0.89 (Fig. 5). However, because of the model used for the Bayesian MCMC approach [e.g., ε(t) γ(t) p(best structure)], the estimates showed much greater annual variation than model-averaged estimates (Fig. 5). Estimates from the MCMC posterior distribution closely reflected raw occupancy patterns (not shown), except that MCMC-based estimates were higher due to the incorporation of detection probability. The model-averaged estimates sometimes diverged from their MCMC counterparts, but the magnitude of the divergence varied from area to area (Fig. 5). However, by the end of the study period MCMC and model-averaged estimates of Δ k were all ≤5% of each other. Variances of the MCMC Δ t estimates were greater than for their modelaveraged counterparts. For the Pradel model, the 95% CI lengths were 21-24% longer for MCMC-based estimates relative to modelaveraged estimates, while for the occupancy model they were 28-96% longer.

DISCUSSION
Although correlations for trajectories of Δ t between territory occupancy and Pradel models were moderate for MCMC-based estimates and high for model-averaged estimates, the magnitude of the trajectories was different with the Pradel model reporting greater magnitude of changes in the population; that is, Δ k from the Pradel model was approximately 20-30% lower than for the territory occupancy model for LAS and SIE and 25% higher for SKC. In terms of declines, territory occupancy declined by 13-21%, while the Pradel model indicated the population declined by 31-44%. Although both methods report declines, estimates for the Pradel model were of greater magnitude which could lead to different management responses. For example, the International Union for Conservation of Nature and Natural Resources (IUCN) guidelines use a decline of >30% over a 10-yr time period as one of its criteria for categorizing threatened populations (IUCN 2001). Thus, results from the Pradel model may lead to conservation actions, while results from the territory occupancy model may result in status quo management. The difference in the magnitude of estimated occupancy model. For example, if two owls share declines can arise simply because of the different a territory, the loss of one may be reflected in resolution between the approaches. In general, survival and, consequently in λ t and Δ t , by the there can be decreases or increases in population Pradel model, but because the territory remains numbers that will not be captured by the territory occupied it is not reflected by the occupancy model. Tempel et al. (2014) noted a similar discrepancy between an integrated population model (IPM) and occupancy model, which they attributed to an increasing number of territories occupied by single owls during the study period. A similar mechanism can result in the differences in estimated increases between Pradel and occupancy methods. To wit, a second owl could move into a territory occupied by a single owl; this change can be captured by the Pradel model via recruitment, but not by the occupancy model. However, a multistate occupancy model that Estimates are from occupancy model ε(t) γ(t) p(best structure) and Pradel model ϕ(t) f(t) p(best structure). Data are for California spotted owls on three study areas in the southern Cascades and Sierra Nevada, California, 1993-2013. Overall realized population change is the proportion of the initial population size or number of territories remaining at the end of the time period † Estimates are the median and 2.5th and 97.5th percentiles from the posterior distribution.
includes the transition probability of a site moving between two occupancy states of pair vs. single bird may be an appropriate way to address this issue. And from a management perspective, a multistate occupancy model might be much more informative, because trends in the occupancy of breeding pairs might be more meaningful to population dynamics than trends in the number of territories with ≥1 owl. Future simulation efforts are needed to evaluate whether this approach would yield trajectories that are more similar in magnitude to the Pradel model. However, there is another mechanism that could give rise to differences between the occupancy and Pradel model for an increasing population; that is, the territory occupancy model inherently has a density-dependent structure, in that no additional territories can be colonized if: (1) the study area is at carrying capacity; or (2) the number of sample territories is fixed at the start of the study and all are occupied with no chance for detecting additional newly colonized sites in the study area. The SKC study area showed increases, with the probability of a ≥10% increase estimated as 0.07 for the occupancy model and 0.74 for the Pradel model (Table 5). Because of the retrospective nature of our study where we assigned occupancy sampling sites after the end of the sampling period, and because we surveyed areas between territories, if the population increase in SKC Pradel estimates led to additional territories, those territories would have been incorporated in our sample. Our occupancy estimate of 0.99 suggests that the number of owl territories on our study area may have reached carrying capacity, while the higher estimate for the Pradel model may indicate an increase in the number of pairs of owls at occupied sites over time.
The choice of territory size (and concomitantly number of occupancy sites for a study area) Estimates are from occupancy model ε(t) γ(t) p(best structure) and Pradel model ϕ(t) f(t) p(best structure). Data are for California spotted owls on three study areas in the southern Cascades andSierra Nevada, California, 1993-2013. Probabilities are based on a posterior distribution of overall realized population change (Δ t ), which is the proportion of the initial occupied territories remaining occupied at the end of the monitoring time period. v www.esajournals.org CONNER ET AL.
could also lead to differences between territory occupancy and Pradel estimates of population change. If occupancy sample sites were smaller than territories, there would be a negative bias in ψ in a declining population because the net loss of one territory holder can result in the loss of two occupied sites. Conversely, if sites were large, such that one site covered >1 territory, then losing one territory holder would not result in an unoccupied site and declines that exist may be missed. This problem has been noted for occupancy studies in continuous habitat by Efford and Dawson (2012), who conclude that occupancy is an inadequate metric for population monitoring when it is confounded with home-range (akin to territory in this context) size. For best results, in terms of tracking population changes through time, the choice of site size should match the territory size.
In addition to differences in magnitudes of the trajectories of Δ t , variance also differed between the methods with smaller variances for territory occupancy estimates compared to Pradel estimates. This may often be the case because detection probabilities are typically greater than capture or resighting probabilities because only a single individual has to be observed and confirmed to species for a detection to occur. In contrast, for a capture or recapture/resighting to occur, an individual has to be physically captured or observed well enough that its identity can be ascertained. This is a much higher bar for most species, and lower capture or detection probabilities almost always result in decreased precision in the parameters of interest. While the larger variance for the Pradel model may be because individual capture and resight probabilities are often lower than detection probabilities, this was not the case here as resight rates were high (90% were ≥0.85). Tempel and Gutiérrez (2013) found a similar result, which they speculated may be because occupancy models have larger effective sample sizes due to multiple surveys at a territory each year, while the Pradel model uses a single point. There is debate over what should be the effective sample size for robust-design closed capture-recapture models (Kendall and Bjorkland 2001), the category of models to which occupancy is most closely related (Barbraud et al. 2003). However, when we changed the effective sample size to reflect the most conservative approach for estimating effective sample size, in this case the number of sites sampled, there was very little change in the variance (and no change in the order of top models). It may also be that this difference is partly a theoretical variance calculation issue. The delta method variance estimate for the Pradel model, in which Δ t is the product of all annual λ t , involves summing logtransformed values of all variances and covariance of each multiplied estimate. Meanwhile, Δ t for occupancy is simply a division for each year that does not include a cumulative effect of added variances. Further investigations comparing occupancy model and Pradel model results should include studies using simulated data to help more fully understand why the occupancy model yields more precise estimates than the Pradel model.
For both occupancy and Pradel models, correlations were high for trajectories of Δ t between Bayesian MCMC-based and model-averaged estimates, and the magnitude of the trajectories was similar. However, there was a difference in the temporal pattern of Δ t between the two approaches. There was more annual variation in the MCMC-based estimates; thus, estimates of Δ t for individual years could differ from modelaveraged estimates (see LAS, Fig. 5). This is not due to an inherent difference between maximum likelihood and Bayesian inference, rather it is because top models for all study areas, and for both occupancy and Pradel models, were less parameterized time trend models. We note that this may not be the case for other studies; modelaveraged values could show high variation if the top models were time varying. Because we were comparing approaches, we were not concerned with comparing specific models for the model-averaging approach. Rather, we used model averaging as the practitioner uses it (e.g., Anthony et al. 2006, Doherty et al. 2010, Tempel and Gutiérrez 2013. However, for the Bayesian MCMC approach, we used a model in which annual parameters were unconstrained: ε(t) γ(t) p(best structure) and ϕ(t) f(t) p(best structure). We chose this highly parameterized model to avoid any bias in annual estimates of ψ or λ, and expected the variance to be reasonable given that we modeled the parameters as random effects (i.e., shrinkage estimates or empirical Bayes v www.esajournals.org CONNER ET AL. estimates), which removes sampling variation. Not surprisingly, estimates from the MCMCbased approach show annual variation better, but their variances were greater than for their model-averaged counterparts. Perhaps a future comparison should be made between likelihoodbased and Bayesian model-averaged estimates.
Although the higher variances associated with Bayesian MCMC-based estimates may be disadvantageous, they have a major advantage over model-averaged estimates; namely, that the probability of decline can be estimated for a retrospective analysis using the posterior distribution of Δ t . In general, evaluating whether a population declined typically relies on hypothesis tests with a null hypothesis of no decline. However, this approach does not answer the real question of interest; "What is the probability the population declined at or more than a specific amount?" (Reckhow 1990, Crome et al. 1996, Wade 2000, King et al. 2010. Using Bayesian MCMC methods to generate a posterior distribution of Δ t is a valuable tool for robustly estimating risk, or probability of declines, for retrospective analysis of monitored populations. Spatial relationships between occupancy and abundance have been an active topic in ecology, particularly at the landscape scale (e.g., Gaston et al. 2000, Webb et al. 2012, Roney et al. 2015, but also at smaller scales (Estrada andArroyo 2012, Guarino et al. 2012). The governing theory is that there will be a positive inter-and intraspecific relationship between regional occupancy and local abundance (Zuckerberg et al. 2009). In addition, other methods use occupancy data to estimate abundance directly, by either using additional information (e.g., counts; Royle 2004 or by exploiting densitydriven heterogeneity in detection probability (Royle andNichols 2003, Stanley and. Rather than simply using occupancy models, these methods may be intermediate approaches, in terms of cost and effort, for the detection of trends in populations. However, none of these studies provide insight as to how temporal trends in occupancy relate to trends in population abundance. Tempel and Gutiérrez (2013) conducted the only study we are aware of in which population dynamics based on occupancy and abundance were directly compared for the same population.
Based on model-averaged estimates, Tempel and Gutiérrez (2013) found estimates of Δ t from occupancy and abundance closely matched in pattern and magnitude. While LAS and SIE show similar patterns and magnitude of decline to the Tempel and Gutiérrez (2013) CSO study site, also in the Sierra Nevada, our results differ strikingly in the magnitude of Δ t between the occupancy and Pradel models. If we compare our declining populations, LAS and SIE, to theirs, we found the overall Δ t was ~20-30% lower for the Pradel model compared to the occupancy model, while their estimates were almost identical (3% higher for the Pradel model). However, more recently, Tempel et al. (2014) found a large divergence in the magnitude of Δ t from an IPM compared to the occupancy model. Here, the discrepancy was similar to what we found between the occupancy and Pradel model, with the IPM estimate approximately 30% lower than the occupancy estimate. They attribute part of the difference to the inclusion of earlier years (which had high recruitment but could not be included in occupancy because it did not meet study design assumptions) and/or an increase in the number of territories occupied by single owls over the study period (Tempel et al. 2014). However, because of the inclusion of the additional years, there is no way to diagnose why their results varied so dramatically from ours. Further field evaluation is necessary to understand why, when, and where Δ t will be aligned between territory occupancy and Pradel estimates.

CONCLUSIONS
Occupancy monitoring has seen increased use over the past decade, to some extent as shrinking budgets have motivated the use of less expensive approaches to monitoring. Occupancy models have successfully been used to rapidly and inexpensively assess population status over large landscapes (Finley et al. 2005, Hristienko and McDonald 2007, Budy et al. 2015. Further, for territorial species, occupancy data have extended beyond describing distributional dynamics and relative abundance to estimating reproduction (Nichols et al. 2007, MacKenzie et al. 2010) and survival (Roth and Amrhein 2010). In addition to meeting budget constraints, developing methods that allow for the use of occupancy data can v www.esajournals.org CONNER ET AL. make large amounts of data available for research into changes in population distributions, as well as for monitoring. Many agencies and institutions have storehouses of historical presenceabsence data that could provide a boon for researchers with the appropriate analytical methods (e.g., Hagler et al. 2011). Although some data may not be applicable because of failure to meet requirements for occupancy modeling (Boshoff and Kerley 2010), it is likely there is a great deal of data that can provide a valuable retrospective evaluation of large-scale and long-term population changes.
While territory occupancy has been successfully used to estimate other vital rates, results from our study suggest it may provide a less sensitive metric for estimating a population's trajectory. From our study, it appears the Pradel model is more sensitive to detecting population changes. However, a multistate occupancy model might match more closely and provide greater insight to population dynamics than trends in the number of territories with ≥1 owl. But, until the relationship between multistate occupancy and population dynamics is evaluated, careful consideration of the trade-offs between cost, the spatial scale of inferential interest, and risk tolerance should be exercised in the development of a monitoring program. If one is mainly concerned with detecting declines in a population, and the population is not at great risk, the financial savings of less expensive occupancy monitoring may be worth the loss of more sensitive population information. On the other hand, if monitoring species of concern, which often are characterized by life history strategies where longevity has been selected at the expense of reproduction (Cardillo et al. 2008), then declines may be missed or underestimated by occupancy monitoring. Smaller populations may decline, undetected by occupancy monitoring, to low levels where the probability of extinction increases to unacceptably high levels due to demographic variability and stochastic environmental events (White 2000a, Lande 2001, Morris and Doak 2002. Although both methods reported declines, the Pradel model estimated greater magnitudes of decline that could lead to different management responses. Namely, results from the Pradel model may lead to conservation actions necessary to avoid high extinction or extirpation risk for small populations, while results from the territory occupancy model may result in status quo management. In addition, we found both Bayesian MCMC-based and modelaveraged estimates of Δ t a robust way to evaluate population trends. However, we recommend Bayesian MCMC-based estimates of occupancy for estimating risk, or probability of declines, for retrospective analyses of monitored populations.