Estimating occupancy dynamics for large‐scale monitoring networks: amphibian breeding occupancy across protected areas in the northeast United States

Abstract Regional monitoring strategies frequently employ a nested sampling design where a finite set of study areas from throughout a region are selected and intensive sampling occurs within a subset of sites within the individual study areas. This sampling protocol naturally lends itself to a hierarchical analysis to account for dependence among subsamples. Implementing such an analysis using a classic likelihood framework is computationally challenging when accounting for detection errors in species occurrence models. Bayesian methods offer an alternative approach for fitting models that readily allows for spatial structure to be incorporated. We demonstrate a general approach for estimating occupancy when data come from a nested sampling design. We analyzed data from a regional monitoring program of wood frogs (Lithobates sylvaticus) and spotted salamanders (Ambystoma maculatum) in vernal pools using static and dynamic occupancy models. We analyzed observations from 2004 to 2013 that were collected within 14 protected areas located throughout the northeast United States. We use the data set to estimate trends in occupancy at both the regional and individual protected area levels. We show that occupancy at the regional level was relatively stable for both species. However, substantial variation occurred among study areas, with some populations declining and some increasing for both species. In addition, When the hierarchical study design is not accounted for, one would conclude stronger support for latitudinal gradient in trends than when using our approach that accounts for the nested design. In contrast to the model that does not account for nesting, the nested model did not include an effect of latitude in the 95% credible interval. These results shed light on the range‐level population status of these pond‐breeding amphibians, and our approach provides a framework that can be used to examine drivers of local and regional occurrence dynamics.


Introduction
Substantial evidence suggests that amphibian declines are occurring worldwide (Wake 1991;Stuart et al. 2004;Adams et al. 2013). As is the case for many taxa, this assessment has been built on the aggregate of observations from a series of separate local studies of population dynamics and expert assessments used to infer large-scale population status. Our objective here is to describe a more synthetic alternative for studying large-scale occurrence dynamics, which allow for data from multiple study areas to be analyzed under a single unified framework. We show how this approach can be used to strengthen conclusions about both regional and local population dynamics.
Frequently, conservation scientists need to make conclusions about changes in population status across large scales, such as the entire range over which that species occurs. Observations for a single site or for a few closely occurring sites have limited value for making inference outside of the local area at which they are conducted. At the same time, it is often prohibitively costly to randomly sample individuals or points across regional scales. As an alternative, many large-scale monitoring programs rely on hierarchical sampling designs where a series of smaller areas are selected for study, and sites are randomly chosen and intensively subsampled from within these areas (Fig. 1). Although it is possible to analyze results from each study area separately, stronger inferences are possible if information is aggregated across individual study areas to simultaneously examine both local and regional variation in population dynamics.
Hierarchical study designs are commonly used to address ecological questions. For example, traditional analyses, such as nested ANOVA with subsampling, have a long history of dealing with the hierarchical nature of data collection. More recent developments in mixed effect modeling have increased the utility of analyzing hierarchical data for a range of applications (Bolker et al. 2009). Although now commonplace in many ecological fields, analyses that account for hierarchical sampling designs have been rare for many demographic studies (Schwarz 2002). Examples include methods used for survival analysis, capture-mark-recapture analysis, and occupancy analysis rely on complex likelihoods that do not easily allow for inclusion of random effects or other approaches for apportioning variance components in nested designs (Schwarz 2002;Williams et al. 2002). The computational challenges associated with including random effects in likelihood-based demographic methods has meant they have lagged behind other analyses in addressing the hierarchical way in which much ecological data are collected.
Recent advances in fitting hierarchical models using Markov chain Monte Carlo (MCMC) have opened the door for analyses of demographic data that account for nested structures (Gelman and Hill 2006;Royle and Dorazio 2008;King et al. 2010;K ery and Schaub 2012) and allow the inclusion of random effects to apportion variance components in nested sampling designs. Implementation of random effects in demographic analyses has most often focused on accounting for individual-level variation in parameters (e.g., King et al. 2006;Ford et al. 2012;Miller et al. 2012). Individual heterogeneity is a significant issue in many cases, and these advances have been important in addressing this source of variability (e.g., Cam et al. 2002). In the case of occupancy models, recent attention has also focused on dealing with the nested nature of occurrence, where occurrence at one level is dependent on occurrence at a higher level (e.g., Nichols et al. 2008;McClintock et al. 2010;Mordecai et al. 2011;Miller et al. 2012). Less effort has been devoted to dealing with the variance and covariance that occurs within and among groups when using spatially nested sampling designs. Perhaps the best examples of accounting for random effects in among-group variation in the population analysis literature are the approaches developed for modeling occupancy of multispecies communities (Dorazio and Royle 2005;Dorazio et al. 2006;K ery and Royle 2008). In these models, among-species Site A Visit 1 Visit 2 Year 1 Year 2 Year 3

Site B Visit 1 Visit 2
Year 1 Year 2 Year 3

Species 1
Site A Visit 1 Visit 2 Year 1 Year 2 Year 3

Site B Visit 1 Visit 2
Year 1 Year 2 Year 3 Species 2 Figure 1. As typical of many large-scale studies, we collected data using a hierarchical sampling design. We sampled for breeding evidence for two species of amphibians at randomly selected wetlands that occurred in representative sample of protected areas in the northeastern United States. In addition to this hierarchical spatial design, samples also occurred within a temporal nested design, with data collected across 8 years with multiple visits occurring within each year.
variation is accounted for using random effects. We demonstrate how that the same approach can and should be used to characterize among-group variation that emerges from hierarchical spatial sampling. We analyzed data on the occurrence of breeding wood frog (Lithobates sylvaticus) and spotted salamander (Ambystoma maculatum) populations from 2004 to 2013 in randomly selected vernal pools that occurred within 14 federally protected lands in the northeastern United States (Fig. 1). Vernal pools (small, isolated temporary wetlands) are important for maintaining biodiversity by providing breeding and foraging habitats and water sources for many wildlife species (Colburn 2004). In particular, vernal pools provide essential breeding and larval habitat for wood frogs and mole salamanders (Ambystoma species). Amphibians breeding in vernal pools may be at particular risk of population decline, as the breeding habitat's small size, short hydroperiod, and isolated nature make them more difficult to locate and protect. Both the wood frog and the spotted salamander breed in vernal pools in proximity to upland forest habitat and are widespread in the eastern United States. The southern extent of the wood frog range is found in the southern Appalachian Mountains of Georgia, while the spotted salamander occurs throughout all of the southeast except Florida (Lannoo 2005). Occupancy rates are lower in habitats in urban landscapes than in intact hardwood forests, and occurrence of both species is maximized in short-duration hydroperiod wetlands (Rubbo and Kiesecker 2005). These are among the first species to breed in the spring, and their vernal pool habitats are typically dry by mid-summer (Colburn 2004); thus, their breeding activity and success is dependent on precipitation and pool hydroperiod (Brooks 2004).
We fit hierarchical occupancy models (MacKenzie et al. 2006;K ery and Royle 2008;Royle and Dorazio 2008;Zipkin et al. 2012) using observations of breeding effort for each species. We demonstrate two approaches to the analysis. Our first goal (and the stated goal of many monitoring programs) was to simply estimate the observed trends in occupancy at area and regional levels. Estimating annual changes in occupancy using linear trends also allowed us to explore factors that may explain different occupancy trajectories among areas. Specifically, we hypothesized that a latitudinal gradient in occupancy trends was present. In the northeastern United States, latitude is a good proxy for two factors we hypothesized would have negative effects on occupancy: urbanization and the expected impact of climate change. In both cases, we expected a greater negative effect in the southern portion of the region where human population densities are much higher, and in the southern end of the species' ranges where warming and drying are likely to have a greater influence on habitat occupancy. Our second goal was to fit models using a dynamic occupancy approach, in which information was shared across refuges about the distribution of transition probabilities and annual differences in these probabilities. This allowed us to make inferences about year-to-year changes in occupancy at both the study area and regional level.

Data collection
We conducted repeat-visit surveys to assess the presence of breeding wood frogs and spotted salamanders in vernal pool wetlands during the spring of 2004-2013. Sites were selected using a nested sampling design with the goal of understanding local and regional dynamics of these species. We collected 2-10 years of data in 14 federally protected areas from throughout the region (Table 1). Protected areas were included in the sampling program based on interest of local resource managers. Although areas were not selected randomly, locations are spatially representative of the region and the variety of federally protected areas in the northeastern United States. Protected areas range from northern Maine south to West Virginia, Maryland, and Washington D.C. Within these protected areas, vernal pool wetlands were chosen for sampling by either (1) random selection from a complete sample of vernal pools (if a comprehensive wetland habitat map had been completed) or (2) sampling random points within the area and mapping clusters of vernal pools (Van Meter et al. 2008). Not all areas participated in all years, dependent on annual priorities of the resource manager and funding availability. Wetlands were sampled in spring to coincide with timing of breeding for the focal species. Timing of visits varied according to latitude to match differences in typical breeding phenologies across the region. Two to four visits were attempted each year, and during each visit, two observers searched each pool for egg masses of wood frogs (Lithobates sylvatica) and spotted salamanders (Ambystoma maculatum) by walking the perimeter. Each observer independently recorded whether or not they observed egg masses for each species during the visit. For analyses, we considered a detection to have occurred during a visit if at least one observer recorded the species. We recorded date for all visits to be used as a detection covariate to control for potential confounding due to variation in the timing of egg-laying and potential predation.

Statistical analysis
To properly account for the nested design under which the data were collected, we analyzed data with hierarchical models fit using a Bayesian framework. The statistical models we use are extensions of basic static and dynamic occupancy models (MacKenzie et al. 2002(MacKenzie et al. , 2003, fit in a Bayesian state-space framework (Royle and K ery 2007;Royle and Dorazio 2008), and extended to include among-area and among-year random effects for the primary parameters in the model. Analyses were analogous to methods used for community modeling of occupancy (K ery and Royle 2008; Zipkin et al. 2012), where amongspecies random effects were replaced by among-area effects.
We present two separate analyses, each formulated to address different monitoring objectives. The first example uses a static estimator of occupancy (MacKenzie et al. 2002), where information about the occupancy status in the current year is assumed to be independent of the status in the previous year. In this approach, two primary parameters are of interest: w ijthe probability a site in the ith area is occupied in the jth year; and p ijk the probability of detecting egg masses at a site in the ith area in the jth year during the kth visit given the site is actually occupied. This approach satisfies the desire to evaluate changes in occupancy probabilities, but does not explicitly address the dynamic processes that lead to annual changes in occupancy. The second example is based on a dynamic occupancy estimator (MacKenzie et al. 2003), where occupancy of a site in a given year is allowed to depend on the occupancy status in the previous year. The parameterization we used for the detection was identical for both examples, but a different parameterization was used to describe occurrence patterns. In the dynamic approach, only the initial occupancy (w ij ) for an area is directly estimated. Occupancy in subsequent years is a function of two transition parameters: c ijthe probability a site that was unoccupied in the ith area in year j was occupied in year j + 1 (colonization); and φ ijthe probability a site that was occupied in the ith area in year j was occupied in year j + 1 (persistence).
We use a common structure for both approaches, incorporating both fixed and random effects using a logit-link function (Bolker et al. 2009). We specified all random effects to be normally distributed on a logit scale. We used the same structure for detection in both examples where detection included an intercept (a), a fixed effect for date (b), and a random effect for differences among combinations of area and year (d). Values of date were centered to have a mean of 0 for each area so that the effect only accounted for within-area variation. The probability of detection for the kth visit to the ith area in the jth year is specified as follows: Other parameters are specific to the examples and are described in the following sections. We also needed to specify priors for the parameters. In all cases, we specified intercept terms for fixed effects (a) to be distributed Uniform (0, 1) on the real scale, fixed slope parameters (b) to be Normal (0, r = 100) on a logit scale, and variance terms for random effects (r p ) to be Gamma (0.1, 0.1).
We fit all models using Markov chain Monte Carlo methods to estimate the posterior distribution under each of our models. MCMC simulations were fit using JAGS v 3.2.0 (Plummer 2012) and run using the statistical package runjags (Denwood 2008) in R v. 2.14.0 (R Core Team 2012). We fit three chains of 20,000 samples after an initial burn in period of 5000 samples for each model. Models were checked for convergence based on trace plots and Gelman-Rubin convergence statistics. We include code for each of the models in the bugs language in the supplementary materials.

Model 1 -Occupancy trend analysis
We use an implicit approach where occupancy in each year is not conditional on the occupancy state in the previous year. This approach is useful when we are interested in directly modeling how occupancy changes across time rather than the values for the transition parameters (persistence and colonization) governing these annual changes. Our model in this case needs to describe occupancy for the ith area in the jth year: The model includes an area-specific intercept, a i , and describes a linear trend across years, where Y j is the continuous covariate year and has mean value of 0 across all years for each area. A random effect d i describes the trend for the ith area. Hyperparameters for the random trend effect are the mean trend for all areas, l, and the standard deviation for the among-area variation in trends, r.
We included a fixed effect to account for the effect of a predictor variable for the area-specific trend parameters: L i is the latitude of the ith area transformed to have a mean of 0 for all areas and b is the effect of latitude on trend. We used a similar set of prior probabilities to those used in the first model: logit(a)~Uniform(0, 1), bÑ ormal(0, 100), l~Normal(0, 100), and r~Gamma (0.1, 0.1).
We examined how including the random effect structure affected the mean and credible intervals of the parameter estimates and the implications of ignoring nonindependence of wetlands within areas. To do this, we ran an additional model without d i and compared trend estimates and credible intervals to the model with this term.

Model 2 -Annual changes in occupancy
We were also interested in making inferences regarding the dynamic parameters governing annual changes in occupancy. To do this, we used a model where the probability of site occupancy is conditional on the occupancy status in the previous year. The persistence probability (φ) is the probability that a previously occupied site in time tÀ1 stays occupied (i.e., the species does not go locally extinct) in t and the colonization probability (c) is the probability an unoccupied site in time tÀ1 becomes occupied in t. These transition probabilities were allowed to vary by year and among areas as random effects. This formulation allows functional flexibility in annual changes while taking advantage of information provided by the previous occupancy status of sites. Rather than separately estimate these values for each study area using a fixed effects design, we assumed these probabilities came from a common distribution with a mean value for all areas and years (logit[a]) and used hyperparameters to describe variation among areas and years (d i and h j , respectively). Persistence for the ith area in the jth year φ ij is equal to For both species, we found that on average southern areas were declining more quickly. Trend estimates for individual areas were generally imprecise except for those with the greatest sampling effort.
Colonization was formulated the same way. Noninformative priors were chosen for both so that logit(a)Ũ niform(0, 1) and r d~G amma(0.1, 0.1) and r h~G amma(0.1, 0.1). In addition, we specified the starting occupancy probability w i1 in the initial year to vary independently among areas (i.e., no hyperparameters) with a prior specified to be w ij~U niform(0, 1).
This approach allowed us to simultaneously make inferences at both the regional and area levels. We were able to estimate a mean occupancy within the region for each of the years, conditional on the assumption that our sample of areas was representative of the public lands in . We examined how including the random effect structure influenced posterior distributions for the effect of a one-degree increase in latitude on trend (b from Model 1 in the text). For both species, we found that when the random effect was included for among-area variation in trend (A and B), the posterior distribution was wider than when the random effect was not included (C and D). Not properly accounting for the nested structure led to overestimation of confidence in estimates. the region as a whole. The random effects structure treats the wetlands within areas as repeated measures of occupancy status within areas and properly propogates sampling error for areas into regional estimates. At the same time, simultaneously estimating occupancy for all areas can improve estimates for individual areas. To the extent that shared factors affect dynamics for all areas, this will be picked up by the random effect structure used in models. This provides added information when estimating effects for individual areas. We used k-fold cross-validation to test whether this was actually the case (Hooten and Hobbs 2014). We compared predictions from the full hierarchical model to a model where hyperparameters were no longer included and, instead, all parameters were treated as fixed effects and estimated separately for each study area. We used 5-fold validation and deviance of the predicted observations as our measure of prediction accuracy. This approach allowed us to determine whether our approach truly improved predictions for data that were held out from the primary data set.

Model 1
Estimated trends for individual areas were imprecise in most cases (Fig. 2). Higher precision occurred for areas with a larger number of sites sampled and where sampling occurred annually across the whole study period. Despite this, we did find some support for a latitudinal gradient in trends across all areas (Fig. 3). The 95% credible interval for the estimated effect in each species did include 0 (À0.056 to 0.129 for wood frog and À0.030 to 0.123 for spotted salamanders). The posterior probability of a positive effect of latitude on trend (i.e., consistent with a positive latitudinal gradient in trends; Fig. 3) was 0.835 and 0.916 for wood frogs and spotted salamanders, respectively. While supportive of our hypothesis, the results are not conclusive because of significant overlap of the credible intervals with 0. However, we would have concluded that there was overwhelming support for the relationship if we had not included the random effect in the model and thus treated each wetland as independent despite the nested structure of data collection. Based on narrower credible intervals, we would have estimated a 0.994 and 0.992 probability of a positive relationship (Fig. 3).

Model 2
At the regional level, there appeared to be a slight decline in wood frog occupancy during the middle years of the study period, with the lowest occupancy occurring in 2008 and 2009 (Fig. 4). Mean occupancy for spotted salamanders was stable with perhaps a minor decline in the final year of the study. Interestingly, there does not seem to be a very strong correlation in annual changes between the two species at the regional level, indicating they may respond differently to annual environmental variation or to different drivers.
Overall turnover was higher for the shorter-lived wood frog than for the spotted salamander. Mean among-area and among-year colonization and persistence probabilities for wood frogs were 0.119 and 0.896, respectively, and the mean values for spotted salamanders were 0.085 and 0.933, respectively.
Although regional estimates of occupancy were fairly consistent, there was significant annual variation in occupancy within the individual protected areas (Fig. 5). In some cases, distinct declines occurred, while in others, occupancy increased across the study period. Areas also varied in their mean occupancy, ranging from having low occupancy in all years to those where occupancy was nearer to 100% in all years. Estimates for each of the areas also show characteristics of the overall regional pattern, reflecting the fact that information was shared across areas due to the random effect structure in the analysis. The strongest shrinkage to the regional mean occurred for areas where sample sizes were small. It is possible to generate estimates in years where sampling did not occur for an area. The characteristics of the area during the sampled years and information provided by patterns in other areas provide a basis for making inferences during these years. We distinguish these estimates in Figure 5 using dotted lines. In general, credible intervals are wider in unsampled years where loss of precision increases with the number of years elapsed since the area was last sampled. Interestingly, there are minimal effects on credible intervals for effects on occupancy in areas where only a year or two are missing in between years where an area was sampled.
(A) Figure 5. We used hierarchical dynamic occupancy models to estimate annual occupancy probabilities and 95% credible intervals for wood frogs (A) and spotted salamanders (B) for each of the participating areas. We only display data for years data were collected at an area using solid lines and for years when data were not collected using dashed line. Estimates are both influenced by the data collected at the specific study area and by the mean value for all areas because of the random effect structure. Spotted salamanders were never observed at Great Swamp National Wildlife Refuge, and therefore, not estimates are displayed for the species at this site. Estimates were generated using Model 2. We found strong support that the hierarchical structure (random effects model) improved predictions when compared to a model where parameters were estimated independently for each study area (fixed effects model). In the case of wood frogs, the deviance for predictions from the hold-out data sets (5-fold cross-validation) was 7213.86 for the random effects model compared to 7255.76 for the fixed effects model (lower deviance indicates better model fit). The improvement was even greater for the spotted salamander model where deviance was 6919.71 for the random effects model compared to 7309.10 for the fixed effects model.

Discussion
As is the case with many species, despite decades of research on both the spotted salamander and the wood frog, there is little information on their status at regional scales. Instead, our knowledge of their status has had to be extrapolated from smaller-scale studies. Existing locallevel information suggests that these species, especially wood frogs, exhibit cyclical dynamics in breeding effort. Petranka et al. (2007) found a decline in occupancy of wood frog and less pronounced decline in spotted salamander egg masses from 1996 to 2006 at both restoration and reference sites in a single mitigation site in North Carolina. Using calling survey data, Gibbs et al. (2005) found populations of wood frogs in NY were stable between surveys conducted from 1973 to 1980and 2001to 2002. In Wisconsin, Trenham et al. (2003 report stable populations of wood frogs over the period from 1981 to 1998. Occupancy of wetlands in the Chesapeake and Ohio Canal National Historic Park have been relatively stable for both species from 2004 to 2013, although other species have been declining . Our expectation was that populations of wood frogs and spotted salamanders occurring in protected areas would exhibit stable occupancy patterns (i.e., no trend in occupancy). Important to this expectation is that, as a metric for population change, we expect occupancy to be less sensitive to annual variability in breeding effort (or cyclic patterns in abundance) than annual abundance. While previous research finds variation in annual abundance, most do not find absence of breeding effort (Berven 1990;Pechmann et al. 1991). Our analyses suggest that on average at the regional level, occupancy probabilities were relatively stable from 2004 to 2013.
However, we did find support for among-area differences in the occupancy dynamics across the period of study (Figs. 2,5). Some areas exhibited substantial increases and decreases for both species. Similarly, trends were directionally consistent with our prediction of a north-south gradient in population trends across the region, with greater decreases in the southern end of the sampling region (Fig. 2). The estimated latitudinal gradient in trends (b) was greater for wood frogs, which is not surprising given that our study area encompasses the southernmost portion of their range (Lannoo 2005). Schwarz (2002) noted more than a decade ago that proper accounting for variance components is a significant and unmet challenge for many demographic analyses. Although there are still significant gaps, the development of Bayesian hierarchical and MCMC methods has opened the door to addressing this deficiency. Other examples of recent extensions in the demographic literature include methods to deal with individual variation (Ford et al. 2012); group effects (Zipkin et al. 2009); path analysis (Cubaynes et al. 2012;Gimenez et al. 2012); experimental design (Schwarz 2002); and spatial autocorrelation (Chelgren et al. 2011;Johnson et al. 2013). The strength of many of the demographic analysis methods has always been the ability to properly weight observed samples to account for nondetected (unsampled) individuals. This has often precluded the use of analytical frameworks in many other fields. The great potential of Bayesian methods for demographic analyses is the ability to address these limitations. Our methods make further progress in this area and are likely to be useful for a range of species and monitoring programs that rely on hierarchical sampling efforts.
By designing a monitoring program with multiple protected areas across the northeastern USA, and incorporating this design explicitly in our analysis, we were able to predict occupancy dynamics and characterize regional patterns in trends than by evaluating populations solely at a regional or local level. Hierarchical models offer a powerful approach for analyzing occurrence data collected using nested designs, both properly accounting for the sampling structure and allowing for inference at and across multiple scales. Further, the utility of our nested monitoring design has clear implications for resource management. Concern is often triggered by population status at large scales, while management generally seeks to improve conditions at the local level . Understanding conditions at the local scale and how they relate to regional patterns is important not only in determining where action is most needed but also to tailor efforts to address local variation. A key step for doing this in our case will be to extend inference to better understand the factors that lead to within-site and among-year variation in occupancy dynamics.