INTRODUCTION

Early studies of population dynamics focused primarily on temporal variation in demographic parameters (Errington 1945, Lack 1964), but a basic principle of modern population ecology is that survival, recruitment, and movement rates vary over both time and space (Holmes et al. 1994, Tilman and Kareiva 1997, Hanski 1999). Understanding the factors that influence spatio-temporal variation in demographic parameters has become a central objective of basic ecological research because this information is needed to answer fundamental questions about the dynamics of species distributions, the mechanisms governing range shifts, and the role of density dependence in population regulation and synchrony (Bjørnstad et al. 1999, Paradis et al. 1999, Pagel and Schurr 2012, Gurevitch et al. 2016). In applied settings, insights about spatial variation in demographic parameters can be used for purposes such as assessing management actions, identifying high-quality habitat, and guiding reserve design (Van Horne 1983, Murphy and Noon 1992, Sanderlin et al. 2012).

Although information about spatial demography is clearly needed, acquiring sufficient data for inference at broad spatial scales has proved difficult. For example, data on survival and recruitment typically come from mark–recapture studies that require substantial effort and financial resources (Williams et al. 2002, Saracco et al. 2010). Consequently, most long-term demographic studies have been conducted over small spatial extents, often at just a few plots (Perrins 1979, Rodenhouse et al. 2003, Sæther et al. 2016). Tremendous amounts of information about temporal dynamics have resulted from these studies, but they provide limited insights into the ecological processes governing spatial variation in abundance and distribution.

In contrast to long-term demographic studies conducted at a small number of plots, efforts to study and monitor the dynamics of populations at broad spatial scales have relied on count-based surveys such as the North American Breeding Bird Survey and the British Breeding Bird Survey (Freeman et al. 2007, Sauer et al. 2017). These studies have been valuable for determining when and where population changes occur, but they provide little information about the demographic processes that contribute to the observed dynamics. This has led to substantial debate about the causes of population declines (Rappole and McDonald 1994, Latta and Baltz 1997). Recent efforts have sought to enhance the value of count data for inference on population dynamics by developing hierarchical models with latent demographic processes (Newman et al. 2006, Chandler and King 2011, Dail and Madsen 2011). This framework makes it possible to study processes such as structured population dynamics, instead of simple temporal trends, which are often the focus of conventional analysis of count data (Buckland et al. 2007, Zipkin et al. 2014a, 2014b). Although these methods represent a major improvement over conventional approaches, not all demographic parameters can be estimated from count data, and it is difficult to account for individual heterogeneity in vital rates and detection parameters.

Integrated population models (IPMs) have been developed out of the recognition of the strengths and weaknesses of both mark–recapture data and survey data in studies of population dynamics (Besbeas et al. 2002, 2003, Brooks et al. 2004, Gauthier et al. 2007, Schaub and Abadi 2011). Conventional IPMs can be described as state-space models, in which a time series model or a matrix population model is used to describe the latent population dynamics, with information about the dynamics coming from more than one type of data, which are typically mark–recapture and survey data (Newman et al. 2014). Combining the 2 types of data often makes it possible to learn about processes, such as immigration, that could not be studied with either dataset in isolation (Schaub et al. 2007, Abadi et al. 2010b). Integrated population models can also be used to assess the demographic contributions to population growth and to answer important questions about the factors affecting species of conservation concern (Koons et al. 2017).

Spatially explicit IPMs were recently developed to extend the scope of conventional IPMs (Chandler and Clark 2014). This framework replaces the matrix population model of IPMs with a latent individual-based model (Diggle 2013, González et al. 2016) describing spatio-temporal variation in abundance and distribution. Although numerous types of data could be accommodated by spatial IPMs, efforts so far have focused on spatial capture–recapture (SCR) data (Efford 2004, Borchers and Efford 2008, Royle et al. 2014) and spatially referenced survey data (Chandler and Clark 2014). The primary advantage of spatial IPMs is the ability to model spatial, temporal, and individual-level variation in survival, recruitment, and dispersal. Other benefits include the potential to model all datasets conditional on the same state process, avoiding the assumption that the datasets are statistically independent (Abadi et al. 2010a). In addition, by retaining the spatial information inherent in the data, spatial IPMs can account for the ubiquitous source of variation in detection rates attributable to the distance between animals and sampling locations.

The objectives of the present study are to demonstrate how spatial IPMs can be fitted to data on avian population dynamics and to expand the spatial IPM framework to directly model spatio-temporal variation in recruitment. To demonstrate, we analyzed 4 yr of constant-effort mist-net data and distance-sampling data collected on Canada Warblers (Cardellina canadensis) over an elevation gradient near the southern limit of the species' breeding range. Populations near low-latitude range limits are predicted to decline and shift upward in elevation in response to climate change (Hampe and Petit 2005, Sekercioglu et al. 2008, Conroy et al. 2011), and as part of a preliminary investigation of this topic, we sought to determine how survival, recruitment, and population growth rates varied over the elevation gradient in our study area.

METHODS

We begin by describing the general spatial IPM framework and then present the application to the Canada Warbler data. Spatial IPMs are hierarchical models with an ecological state model and at least 2 observation models describing how the data arise conditional on the latent state variables (Royle and Dorazio 2008, Chandler and Clark 2014). The ecological state model describes how abundance and spatial variation in density (i.e. distribution) change over time as functions of survival, recruitment, and movement. A spatio-temporal point process model is used as the state model because it provides a framework for modeling individual-level variation in vital rates. An overview of spatio-temporal point process models is beyond the scope of the present study, but briefly, these models are designed for inference about the processes influencing the number and spatial distribution of points over time (Diggle 2013, González et al. 2016). They can therefore be used as individual-based models of population dynamics in which the points are animal locations that enter, exit, and move within the spatial region via recruitment, mortality, and dispersal. The observation models in spatial IPMs describe how the capture–recapture and survey data (and possibly other types of data) arise conditional on the abundance and distribution of individuals during each sampling occasion. The observation models would not be necessary if it were possible to monitor each individual in the population over its lifetime. Spatial sampling and imperfect detection make this impossible in most studies, but hierarchical models allow for inference on the underlying state process even though it is only partially observed (Royle and Dorazio 2008).

We use the following notation throughout. Probability density (and mass) functions are represented using the p(·) notation. For example, if the random variable x is Poisson distributed with rate θ, we would write p(x) = Pois(θ). Conditional probability densities will be denoted by p(·|·). Continuing with the previous example, if the random variable y is binomial and depends on the random variable x, we would write p(y|x) = Bin(x,p). Note that the binomial probability here is p, which should not be confused with the density function p(·).

Initial Abundance and Distribution

The first time period is modeled differently than subsequent periods because there is no information about the survival and recruitment processes that gave rise to the initial population. We therefore model the initial period using a spatial point process that is independent of the other periods. Spatial point process models are characterized by an intensity function μ(s) that describes the expected number of points in an infinitesimally small area located at s. In other words, μ₁(s) is the density surface at t = 1 that we wish to estimate. Here, s, without subscripts, represents an arbitrary point in the two-dimensional spatial region 𝒮 ⊂ ℝ². If density is uniform throughout 𝒮, the state space of s_it, the point process is said to be homogeneous with the expected number of individuals given by E(N₁) = μ₁(s)|𝒮| where |𝒮| is the area of the state space. If density varies throughout 𝒮, the point process is inhomogeneous. Spatial variation in density can be modeled as a function of covariates, for example with a log-linear model: log (μ₁(s,β)) = v′(s)β, where v(s) is the set of spatially referenced environmental variables, which typically are formatted as “raster” data, and β is the vector of coefficients of the log-linear model. Although the intensity function (μ₁(s)) depends on the β coefficients, we will suppress them to be concise.

When the point process is inhomogeneous, the expected value of initial abundance is found by integrating the intensity function over the spatial region:

Equation 1 is also used in the data augmentation scheme to compute the probability that an individual is a member of the initial population: p(z_i,1) = Bern (ψ₁ = Λ₁ / M). The distribution of activity centers is proportional to the intensity function, and the probability density function for a single point is found by normalizing: p(s_i,1) = μ₁(s_i,1)/Λ₁. If the points are mutually independent, after accounting for covariate effects, the joint density of the point pattern is given by the product of M such terms. Otherwise, the spatial dependence in point locations could be modeled with a Markov point process (Reich and Gardner 2014).

Survival, Recruitment, and Dispersal

In subsequent time periods (t > 1), the probability of being alive is modeled conditional on the individual's previous state. If it was alive in the previous period, it survives with probability φ. If the individual was not alive and had not been previously recruited, then it becomes recruited with probability γ˜_t. These processes are modeled as

where

The variable a simply indicates whether an individual is available to be recruited. This formulation is now commonly used in individual-based extensions of the Jolly-Seber model (Royle and Dorazio 2008, Royle 2009). However, unlike other open-population capture–recapture models and spatial IPMs, we model the probability that an individual is recruited (γ˜_t) as a function of a point process for the number and locations of recruits:

where γ_t(s) is the per capita recruitment rate at location s. As with the point process for initial abundance and distribution, the recruitment intensity function can include spatial covariates. The intensity function could also be expanded to include density dependence.

In the absence of dispersal, spatial variation in density at t > 1 is determined by spatial variation in survival and recruitment: μ_t(s) = φ_t−1(s)μ_t−1(s) + γ_t−1(s)μ_t−1(s); and the expected population growth rate is λ_t(s) = μ_t(s)/μ_t−1(s) = φ_t−1(s) + γ_t−1(s). The probability that an individual is recruited can be computed by dividing the expected number of recruits by the number of individuals that are available to be recruited: γ˜_t = Γ_t/A_t where A_t = .

The recruitment point process describes the distribution of newly recruited individuals, but when dispersal is possible, the distribution of individuals that were recruited prior to t requires an additional model. Any dispersal process could be considered, but here we focus on dispersal kernels of the form κ_t(s,s_{i, t−1}, τ_t), which could be Gaussian, negative exponential, etc. The kernel specifies the relative probability of dispersing to location s at time t, conditional on being previously located at s_i,t−1, with τ being the vector of parameters that govern the shape of the kernel. The probability density function for s_it is therefore a mixture that depends on whether the individual is available to be recruited:

This probability density is proportional to either the dispersal kernel or the intensity function, meaning that the actual density requires computing a normalizing constant found by integrating either function over the spatial region 𝒮.

The recruitment point process and the Markovian survival process determine the expected value of abundance at time t > 1 according to E(N_t) = E(R_t) + E(S_t). The expected number of recruits was defined previously, and, in the case that survival probability is constant among individuals, the expected number of survivors is E(S_t) = E(N_t−1)φ_t−1. Although the expected value of abundance is useful for prediction, and the observed data can be modeled conditional on the expected values, interest will often be on realized abundance—the number of individuals actually alive in 𝒮 at time t—which is given by . Realized and expected values of abundance can also be computed for any region within the state space. For example, the realized number of individuals in region ℬ ⊂ 𝒮 is , where I(·) is the indicator function returning 1 if its argument is true and 0 otherwise.

In nonspatial open-population capture–recapture models, when permanent emigration is possible, φ must be interpreted as “apparent survival,” defined as the probability of surviving and not permanently emigrating from the study area. However, as with other recently developed open-population SCR models (Ergon and Gardner 2014, Schaub and Royle 2014), our model provides an opportunity for estimating actual, instead of apparent, survival because movement and survival can be distinguished. However, if there aren't enough data to estimate dispersal, and hence the probability of permanent emigration, it may be necessary to assume that activity centers do not move among years. In this case, φ should be interpreted as apparent survival, with permanent emigration being the case where an individual permanently moves to an area where its encounter rate is negligible.

To complete the state model, we introduce one more partially observed latent variable, u_ikt, which denotes the location of individual i during secondary sampling occasion k in primary sampling occasion t. Distinguishing between primary and secondary sampling occasions corresponds to the “robust design,” in which it is assumed that mortality and recruitment are negligible during, but not among, primary periods (Pollock 1982). By collecting replicate observations within primary periods, it becomes possible to account for more sources of heterogeneity in capture probability. In our case, it also provides a way of modeling the within-season movement process that influences where individuals are detected during distance-sampling surveys. A simple model for u_ikt assumes that the locations are independent bivariate normal outcomes centered on the individual's activity center: , where σ is the scale parameter related to home-range size. Alternatively, a Markov movement model such as an Ornstein-Uhlenbeck movement model could be used to account for serial correlation and describe space use (Blackwell 1997, Hooten et al. 2017).

Observation Models

Capture–recapture data. Nonspatial, open-population, capture–recapture (CR) models are well-developed for “capture history” data in which y_ikt is a binary variable indicating whether individual i; i = 1, … , n0 was captured on secondary sampling occasion k; k = 1, … , K within primary sampling occasion t (Jolly 1965, Seber 1965, Pollock 1982). Such data do not retain information about the location of capture, and nonspatial CR models ignore the spatial region within which the population occurs. By ignoring space, these models make it difficult to account for variation in capture probability among individuals that arises from the distance between animals and traps. Of greater importance from the perspective of ecological research is the fact that nonspatial CR models do not allow for inference on spatial variation in density and other demographic parameters, except when space can be naturally subdivided into a small number of discrete units (Nichols and Kendall 1995).

In SCR models, the location of capture is an important component of the data structure. Although numerous types of data are suitable for SCR models, the data are commonly in the form of a four-dimensional array, in which each element, y_ijkt, is a binary variable indicating whether individual i was captured or encountered at trap j; j = 1, … , J on occasion k in year t. The location of each trap is stored in a J × 2 matrix, with x_j representing the two-dimensional Cartesian coordinates of trap j. The simplest nonspatial CR models treat capture probability p as a scalar, whereas most SCR models assume that capture probability is a function of the distance between individuals and traps (Efford 2004, Royle et al. 2014). The distance metric used in SCR models is not the distance between an individual's actual location (u_ikt) and a trap but is instead the distance between an individual's activity center (s_it) and the trap. The reason for this is that u_ikt is unknown on occasions when the individual is not captured, and accounting for this would require an explicit movement model, which would add an unnecessary degree of complexity in studies that are not focused on fine-scale movement behavior (Borchers 2012).

Numerous distance-based capture probability functions can be considered, but a common choice is based on the kernel of the Gaussian distribution: , where is the Euclidean distance between the activity center and the trap. The parameter p₀ is capture probability at zero distance (i.e. capture probability when an individual's activity center is coincident with a trap). The scale parameter σ determines the rate at which capture probability decreases with distance. This is the same scale parameter that can be used in the bivariate normal movement model described previously. As with all CR models, additional temporal and behavioral effects on capture probability could be considered.

If the encounter history data are binary, they could be modeled as Bernoulli outcomes: . The inclusion of ensures that an animal can only be captured if it is alive during primary period . The Bernoulli model is one of several possible encounter models, and it assumes that an animal can be detected at multiple traps during each secondary sampling occasion. This implies that “traps” need not be standard live traps in which individuals are physically restrained. Instead, SCR models can be used with data from camera traps, hair snares, acoustic recorders, and other types of passive detectors, for which the Bernoulli assumption is not always ideal. In the case of mist-net data, individuals are physically restrained and, hence, they may not be captured at more than one net per occasion. For this type of data, a categorical observation model can be used in place of the Bernoulli model. Data for the categorical model are formatted such that y_ikt indicates the trap index at which individual i was captured on secondary occasion k in primary period t. Or, if an individual was not captured on a particular occasion, y_ikt is recorded as J + 1. The model then becomes p(y_ikt|s_it, z_it) = Categorical(π_i,1,t, … , π_i,J+1,t), with cell probabilities constructed using the multinomial logit (Royle et al. 2014:254–256). Although the categorical model is often recommended for mist-net data, the Bernoulli model may be more appropriate if a secondary sampling occasion consists of many hours of netting, such that an individual could be captured at more than one location within an occasion. If these “same-day recapture” events are independent, then the Bernoulli model could be used. Otherwise, the same-day recaptures could be discarded before implementing the categorical model.

Distance-sampling data. In distance sampling, individuals are surveyed from either points or line transects, and detection probability is assumed to decline with radial or perpendicular distance from the transect (Buckland et al. 2001). Conventional distance-sampling models are nonspatial, with the observed data being the distances to the subset of detected individuals. Spatial distance-sampling models are different, in that detection is modeled conditional on individual locations, which makes it possible to directly model the spatial distribution of individuals (Royle et al. 2014, Borchers et al. 2015). Below, we describe how either location or distance data could be modeled conditional on the point process model.

Let x˙_l reference the lth survey location in a collection of L point-count sites. The probability of detecting individual i is based on a monotonically decreasing function of distance between the observer and the animal, such as , where the dot notation distinguishes distance-sampling parameters from their SCR counterparts. The probability of detecting an animal at zero distance, p˙₀, is usually taken to be unity because animals can often be assumed to be detected with certainty at the survey point. If the actual locations of individuals are recorded when performing distance sampling, the observed data would be the subset of u_ikt detected from the L points. However, reconciling these data with the SCR data cannot be done with certainty because the identity of the individuals detected during distance sampling will almost certainly be unknown, because marks (such as color bands) cannot easily be seen when conducting surveys. A model is therefore needed to account for the unknown identity of the detected individuals if the SCR and distance sampling are to be modeled conditional on the same realized point pattern.

To simplify the observation model, and to reflect the nature of many distance-sampling datasets, assume that the distance-sampling data are binned into B annuli within the radius of the point-count plot. The number of individuals in each annulus will be denoted by . The data in this case could be modeled as binomial counts: , with being the average detection probability within each annulus, computed by integrating over the distance cut points (c₁, … , c_B+1):

where r is the radial distance from the observer.

This approach to modeling the distance data conditional on s_it and u_ikt is extremely computationally challenging, even if u_ikt is marginalized as described by Royle and Dorazio (2008:238–241). A much more computationally appealing approach is to model the distance-sampling data conditional on the underlying density surface, μ_t(s), instead of on the realized point pattern. One option for doing this is to assume that density within each point-count plot is uniform and that the number of individuals in each plot is Poisson distributed. This leads to a Poisson model for the distance-sampling data in which the expected number of detections at point l in distance bin b is given by the expected number of individuals in bin b multiplied by (Royle et al. 2004, Sillett et al. 2012).

Canada Warbler Data and Model Specification

We collected spatial capture–recapture and distance-sampling data on Canada Warblers during May–July, 2014–2017, in Nantahala National Forest, Macon County, North Carolina, USA (35°4′35′′N, 83°28′42′′W). The study area was selected because of the pronounced elevation gradient (600–1,600 m) within a relatively small area (63 km²) spanning the range boundaries of many high-elevation bird species, including the Canada Warbler. In the southern portion of their range, Canada Warblers primarily occur above 1,000 m, where they establish territories in hardwood forests with a dense understory of woody vegetation (Reitsma et al. 2009, Becker et al. 2012). Sampling along the elevation gradient and across the range boundary was conducted as part of a broader effort to understand the factors influencing population dynamics near the southern limits of the breeding range.

We collected SCR data at 9 sites, 8 of which were sampled in more than 1 yr (Figure 1). At each site, 20 net locations were established and arranged in 4 rows, with the outer rows spaced by 50 m and the inner 2 rows spaced by 100 m. Five nets were placed in each row, and midpoints of nets within rows were separated by 25 m. In general, only 10 of the 20 nets were operated each day, for ∼6 hr day⁻¹ beginning 30 min before sunrise, from May 1 to July 1. The 4 days of sampling were usually consecutive, except when weather interfered. Nets were not operated during rain or high winds. To increase capture rates, we broadcast Canada Warbler vocalizations using speakers placed at 2 nets on days 2 and 4 of the netting session. We used nylon mist nets (32 mm mesh, 12 m long, 1.5 m tall) that were checked every 30–60 min, and each captured individual was marked with a USGS aluminum band and 3 color bands. Age, sex, and morphological measurements were also recorded, and birds were released within 30 min of being extracted.

FIGURE 1.

Location of mist-net plots and number of male Canada Warblers captured at each plot in Nantahala National Forest, Macon County, North Carolina, USA. Background colors represent elevation (m). Only 2 sites were sampled in 2014. Not shown is individual-level information about within- and among-season recaptures. UTM = Universal Transverse Mercator, zone 17N.

Distance-sampling data were collected at 70, 71, 109, and 109 points in 2014, 2015, 2016, and 2017, respectively. Survey locations were arranged on a 500 m grid covering the study area (Figure 2). Each sampled location was visited before 1100 hours once per year, between late April and early July. Surveys lasted 10 min and were divided into four 2.5 min intervals. For each detected individual, we recorded the species, detection cue (e.g., song, chip, visual), time intervals of detection, and location on a simple map with 10 annuli defined by radii with 10 m distance increments. Prior to data collection, observers were trained in distance estimation using range finders. Distance data were collapsed into 5 distance bins defined by 20 m increments from 0 to 100 m. By recording both distance and time to detection, we were able to estimate the effect of distance on detection probability, as well as the effects of other covariates on p˙₀, which in this case can be defined as the probability that an individual vocalized (Alldredge et al. 2007, Chandler et al. 2011, Sólymos et al. 2013).

FIGURE 2.

Location of point-count plots and number of male Canada Warblers detected at each plot in Nantahala National Forest, Macon County, North Carolina, USA. Background colors represent elevation (m). Thirty-eight new plots were added in 2016. Although only total counts are shown, time-to-detection and distance-sampling data were collected at each point. UTM = Universal Transverse Mercator, zone 17N.

We modeled initial density, survival, and recruitment as functions of elevation. Although numerous other abiotic factors and biotic interactions can influence survival and recruitment, many of them are likely to covary with elevation, and our intent here was not to tease apart these effects but rather to demonstrate how spatial IPMs can be used to draw inferences about spatial variation in demographic parameters. We used an adult, male-only model without dispersal because >90% of our point-count detections were of males and because we did not observe marked individuals moving >200 m between years (see below). We considered 3 models for the relationship between Canada Warbler demographic parameters and elevation. The first model treated each parameter (initial density, survival, and recruitment) as a linear function of elevation on the link scale. Log link functions were used for initial density and recruitment, and a logit link was used for survival. The second model allowed for quadratic effects of elevation on the link scale to allow for possible non-monotonic relationships between demographic parameters and elevation (Lichstein et al. 2002). The third model used a logistic equation for each demographic parameter to evaluate the possibility that demographic rates had upper asymptotes. For example, the logistic equation for initial density was , where is the asymptote and determines the rate at which density approaches the asymptote over the elevation gradient. We aggregated a 30-m-resolution digital elevation model to 180 m resolution, and we assumed that survival, recruitment rate, and density were constant within each grid cell. This allowed us to evaluate the integrals in Equations 1 and 4 using a Riemann sum.

We augmented the SCR data using M = 350, and we used a Bernoulli observation model because we had several same-day recaptures, which are incompatible with the categorical observation model. Variation in capture probability was modeled as a function of distance between activity centers and nets, and we also accounted for the playback effect by modeling a unique baseline capture probability, p₀, for passive and playback net occasions. To account for the fact that not all nets were operated on each occasion, we fixed capture probability to zero for net sites that were not operational on a particular occasion. The state space was created by placing a 400 m buffer around the trap locations.

The distance-sampling and time-of-detection data were modeled conditional on the expected value of abundance at each point-count plot, which was computed by area expansion. Specifically, we multiplied local density, μ(s), by the area of the 100 m radius point-count plots. We assumed that time of detection and detection distance were independent, which allowed us to model the data using 2 independent multinomials. The multinomial cell probabilities for the distance-sampling data were based on Equation 6. The time-to-detection cell probabilities were p˙₀, (1 − p˙₀)p˙₀, (1 − p˙₀)²p˙₀, (1 − p˙₀)³p˙₀, where p˙₀ is the probability that an individual vocalizes during a 2.5 min interval. We modeled variation in p˙₀ using a logit link with time of day and date (days since May 1) as covariates.

We fitted models using JAGS 4.2.0 run with the “rjags” package in R 3.4.1 (Plummer 2003, R Core Team 2017). Vague prior distributions were used for all parameters (Appendix 1). We created 3 Markov chains, each consisting of 10,000 iterations, and we discarded the first 5,000 iterations as burn-in. Convergence was assessed using visual inspections of the Markov chains and with the Gelman-Rubin diagnostic (Gelman and Rubin 1992). Model selection was based on posterior deviance. Code is provided in Appendix 1.

RESULTS

We captured 109 adult male Canada Warblers during the 4 yr study. Seventy-six individuals (69.7%) were captured just once, whereas 21 birds were captured twice, 6 were captured 3 times, 4 were captured 4 times, 1 was captured 5 times, and 1 was captured 6 times. Twelve individuals (11.0%) were captured in 2 consecutive years, and no individuals were captured in >2 yr. All within- and among-year recaptures were within 200 m of the original capture location, and no individuals were captured at more than 1 of the 9 sites, which suggests that site fidelity was high and adult dispersal low. The total number of individuals captured, averaged over the 4 yr, increased with elevation (Pearson's r = 0.68). No Canada Warblers were captured at the 2 lowest sites, and the greatest number of captures occurred at the highest site (Figure 1). We detected 30, 25, 62, and 45 Canada Warblers on our point-count plots during the 4 yr of the study. The proportions of sites with ≥1 detection were 0.43, 0.35, 0.57, and 0.41 from 2014 to 2017. Some of the variation among years was due to the addition of 38 new survey locations in 2016. Annual counts and observed occupancy tended to increase with elevation, but the relationship was weaker than in the mist-net data (Pearson's r = 0.41; Figure 2). Below 1,000 m, Canada Warblers were detected at only 2 of the 30 point-count plots. The proportion of point-count plots with detections was higher above 1,000 m (38 of 79 plots), yet in contrast to the increasing trend observed in the mist-net data, observed occupancy and average counts were relatively constant above 1,000 m (Figures 1 and 2). Detection frequencies in the five 20-m-wide annuli were 37, 53, 30, 30, and 12. Note that detection frequencies would have increased with distance if detection probability was 1, because the area of each annulus increases with distance. The frequencies of first detections during the four 2.5 min intervals were 109, 26, 18, and 9, indicating that most birds were first detected early in the 10 min surveys.

All 3 models converged, with Gelman-Rubin statistics <1.1 for each parameter. The choice of M = 350 was deemed sufficient because the estimated probability that the super-population size was greater than 300 was <0.001. The most supported model was the logistic model, in which each ecological process (initial density, apparent survival, and recruitment) was modeled with an upper asymptote and a lower boundary at zero (Table 1). The mean deviance for this model was 35 units less than the next most supported model (the linear model), and it included the same number of parameters. The quadratic model had 3 more parameters than the other 2 models, but these parameters did little to reduce the mean and variance of the posterior deviance.

TABLE 1.

Deviance statistics used for model selection. The logistic model was selected as the best model in the set because it had lower mean deviance than the other models, and it was tied for the lowest number of parameters. The number of parameters does not include the number of latent variables, which was constant among models.

Initial abundance increased rapidly with elevation, as indicated by a positive value of and a 95% credible interval (CI) that did not include zero (Table 2). The effect can be seen clearly in the 2014 density surface (Figure 3). Density was <0.03 males ha⁻¹ below 1,000 m and increased more than sixfold to 0.20 males ha⁻¹ at high elevations in the study area.

TABLE 2.

Posterior summary statistics for the most supported (logistic) model, including 95% credible intervals (CI).

FIGURE 3.

Estimated density surfaces for the 4 yr of the study. Color indicates number of male Canada Warblers per hectare; gray contour lines show elevation. In years 2015–2017, density was modeled as a function of spatial variation in apparent survival and recruitment. Data are from 109 point-count locations (+) and 9 mist-net sites (×) in Nantahala National Forest, Macon County, North Carolina, USA. UTM = Universal Transverse Mercator, zone 17N.

Apparent survival increased with elevation, from 0.040 (95% CI: 0.0032–0.12) at 900 m to 0.29 (95% CI: 0.16–0.42) at the highest elevation (Figure 4). The 95% CI for , the effect of elevation on apparent survival, did not include zero (Table 2). There was much less evidence that per capita recruitment varied over the elevation gradient. The 95% CI for included zero, and the effect size was small (Table 2 and Figure 4). Point estimates suggest that recruitment was not high enough to offset (apparent) mortality, as posterior means of λ_t(s) were <1 over the elevation gradient (Figures 4 and 5). However, 95% CIs include 1, indicating that there is not enough evidence to conclude that the population was declining.

FIGURE 4.

Estimated relationships between Canada Warbler vital rates and elevation over the range of elevations in our study area in Nantahala National Forest, Macon County, North Carolina, USA. Apparent survival (φ) is the probability of surviving and remaining in the study area. Recruitment (γ) is defined as the per capita rate at which new individuals enter the population. Population growth rate is λ = γ + φ. See text for details.

FIGURE 5.

Spatially explicit population growth rates (λ_t(s)) for Canada Warblers in Nantahala National Forest, Macon County, North Carolina, USA. Point estimates (posterior means) were <1 throughout the study area, but 95% credible intervals (CI) included 1. UTM = Universal Transverse Mercator, zone 17N.

Playback had a large effect on capture probability. Without playback, the baseline capture probability (p₀) was 0.029, compared to 0.156 at nets with playback (Table 2). The scale parameter (σ) of the capture probability model indicated that capture probability was negligible at nets that were >500 m from an individual's activity center (Figure 6). The probability that a bird vocalized during the 10 min point count was 0.96 (Figure 6) and was not strongly affected by either date or time of day (Table 2). Canada Warbler detection probability decreased rapidly with distance, being approximately zero at 100 m (Figure 6).

FIGURE 6.

Estimated encounter functions for the point-count data (A) and mist-net data (B) for Canada Warblers in Nantahala National Forest, Macon County, North Carolina, USA. “Distance” in A is the distance between an observer and a bird. “Distance” in B is the distance between a mist net and an activity center.

DISCUSSION

We described a spatially explicit, individual-based model that allows for population-level inference from capture–recapture and survey data. Unlike previous spatial IPMs that focused on temporal variation in recruitment (Chandler and Clark 2014), we developed an approach for modeling spatio-temporal variation by adopting a point process for the abundance and distribution of recruits in each year. The model can be used to assess the effects of environmental variables on demographic parameters in studies of population dynamics, and it should be useful for informing conservation decisions because it provides a means of identifying high-quality habitat, defined as environmental conditions where survival and recruitment are highest (Van Horne 1983).

Our model yields maps of spatio-temporal variation in density, and it can therefore be viewed as a type of hierarchical species distribution model (Hefley and Hooten 2016). However, most species distribution models are designed for presence-only data or survey data, whereas our spatial IPM accommodates data on marked individuals, allowing for insights into the demographic processes that contribute to changes in species distributions (Schurr et al. 2012, Normand et al. 2014). Although numerous new methods have been developed for using data on unmarked animals to study the effects of environmental variables on species distributions and population dynamics (Dail and Madsen 2011, Sollmann et al. 2015, Nadeem et al. 2016), data on unmarked animals provide much less direct information about demographic processes than capture–recapture data. For example, if a population is at equilibrium, count data collected on unmarked animals cannot be used to determine whether mortality is equal to recruitment or mortality and recruitment rates are low. By accommodating data on marked animals, we can distinguish between these possibilities and obtain precise estimates of vital rates (Zipkin 2014a). In particular, monitoring marked animals over multiple years provides much more information about survival than can be obtained from surveys of unmarked individuals. Although less informative than the data on marked individuals, the survey data do provide important information about spatial variation in apparent survival and recruitment. For example, if abundance is increasing in some regions of a study area but not in others, it must be the result of spatial variation in recruitment. Similarly, declines in abundance indicate that recruitment is not sufficient to offset losses due to mortality or emigration. Nonetheless, if resources were unlimited, the survey data would not be needed because SCR data provide all of the required information about distribution and demography. The primary reason for utilizing the survey data is that it is typically cost prohibitive to collect SCR data at a fine spatial resolution over a sufficiently large region to characterize spatial population dynamics. Combining capture–recapture with survey data, which can be collected with much less effort, is therefore desirable. This intuitive idea has been recognized by avian ecologists for many years, but only recently has it been made possible in a statistical inference framework (Rappole et al. 1998, Ahrestani et al. 2017).

The state model in our IPM is a spatio-temporal point process, in which the points (individuals) enter and exit the population via recruitment and mortality. Point process models have a long history in ecological research for datasets in which the points are directly observed (Stoyan 1982, Rathbun and Cressie 1994). However, in our case, the points are unobserved activity centers whose locations are inferred from the capture and detection locations. The use of a spatio-temporal point process model allows for inferences about individual, spatial, and temporal variation in demographic parameters, and the inferences are scalable from the individual to the population level and from small to large spatial regions. Spatial scaling is possible because the intensity functions vary continuously in space and can be used to predict demographic parameters at any point in the study area, yet they also can be averaged over larger regions that may be of interest to managers. Among the benefits of modeling individuals is that it alleviates problems such as the ecological fallacy—drawing inferences on individuals from group data—and the “modifiable aerial unit problem” in which inferences become strongly scale dependent when modeling aggregated spatial data (Robinson 1950, Openshaw and Taylor 1979, Clark et al. 2011). Point process models avoid these dilemmas by modeling the aggregated data (e.g., the observed counts) conditional on the individual-level latent variables. In addition to avoiding bias due to aggregation, these models provide an opportunity for understanding the spatial scale at which environmental variables most strongly influence ecological processes (Chandler and Hepinstall-Cymerman 2016). Moreover, point process models allow for statistical inference and therefore overcome many of the limitations of conventional individual-based models that are often criticized as being too ad hoc and complex to allow for generalizable insights from empirical data (Grimm 1999, Hooten and Wikle 2010).

Although spatial IPMs allow for individual-level variation in demographic parameters, the only source of individual heterogeneity that we focused on was that due to location. Location is clearly important because individuals in different environments often have different probabilities of surviving and reproducing, but we ignored other important sources of variation arising from differences in sex, age, and individual traits, which are the focus of most matrix population models and integral projection models (Easterling et al. 2000, Caswell 2001, Ellner and Rees 2006, Ghosh et al. 2012). Future work could accommodate these sources of variation by borrowing ideas from the nonspatial capture–recapture literature in which individual-level covariates are modeled as partially observed latent variables (King et al. 2008). For example, data on sex-specific detection rates could be used to estimate the sex ratio of the entire population. Distinguishing between males and females would make it possible to estimate fecundity, which would allow for the assessment of hypotheses about reproduction. In addition, a sex-structured model with dispersal would allow one to understand the relative contributions of fecundity and dispersal to recruitment (Ergon and Gardner 2014, Schaub and Royle 2014), which would be difficult to achieve with matrix or integral projection models (Merow et al. 2014).

One limitation of the proposed modeling framework is computation time. Excessive run times (>1 wk model⁻¹) led us to simplify our model by treating the SCR and distance-sampling data as conditionally independent. Specifically, we modeled both datasets conditional on the same intensity function, but not on the same realized point pattern. The latter option is straightforward but very computationally challenging with distance-sampling data because each detection is conditional on the unobserved activity centers and movement events of all individuals in the state space. By ignoring the identity of the individuals detected in our distance-sampling surveys, we avoided much of the computational burden, but we discarded some information about the locations of activity centers and the scale parameter , which is associated with territory size. Ignoring this information should have no effect on inference if the capture–recapture and survey locations are far enough apart to ensure that individuals captured in mist nets are not detected during distance sampling. Future work should attempt to develop computationally efficient methods for implementing the conditional modeling approach to account for this form of statistical dependence.

Our application of the spatial IPM to the Canada Warbler data provided several interesting insights. First, apparent survival was low compared to that of other long-distance passerines (Sillett and Holmes 2002), with <40% of individuals surviving and remaining within the study area. This low estimate resulted, in part, from the fact that we never captured an individual in >2 consecutive years. It is unlikely that permanent emigration contributed substantially to this low estimate because, like most passerines, Canada Warblers exhibit high site fidelity (Hallworth et al. 2008, Cline et al. 2013) and none of the marked birds in our sample were observed to move >200 m among years. It is more likely that events during the nonbreeding season and interactions among seasons contributed to the low apparent survival (Sillett and Holmes 2002, Small-Lorenz et al. 2013, Rockwell et al. 2017). Although apparent survival was low, and the point estimates for population growth suggested that the population was declining, the credible intervals for population growth rates included 1, indicating that recruitment may be high enough to offset low survival. However, with only 4 yr of data, more research is needed to determine whether this population is undergoing a long-term decline as predicted by climate-change models (Matthews et al. 2004).

We found evidence that the pronounced density gradient in our study area was the result of apparent survival increasing with elevation. Both the SCR and distance-sampling data indicated that individuals at lower elevations returned less often than individuals at higher elevations, either because they died or because they dispersed beyond our study area. As more years of data on dispersal become available, it should be possible to distinguish between losses due to mortality vs. emigration. In contrast to apparent survival, per capita recruitment did not exhibit substantial variation over the elevation range. Rather, the number of recruits in an area was determined primarily by local density in the previous year. For example, at low elevations where density was low, recruitment was also low, as can be seen by the consistent lack of detections at the point-count plots below 1,000 m.

Additional research is needed on the mechanisms that govern the relationships between elevation and apparent survival and recruitment. Clearly, it is not elevation per se that influences vital rates, and future work should attempt to determine how demographic parameters are affected by environmental variables that covary with elevation. For example, modeling the direct effects of weather variables could provide valuable insights, because recent research has demonstrated that birds do not always shift their distributions upward in elevation in response to climate change, but instead may be more likely to track precipitation and temperature patterns (Tingley et al. 2009, 2012). More research is also needed to assess the relative influences of dispersal and fecundity on spatial variation in recruitment. The modeling framework developed here provides a platform for addressing these research questions and advancing knowledge about the factors that influence spatial population dynamics.