Spatial survival analysis accounts for female-biased breeding dispersal and provides realistic estimates of true annual survival in migratory warblers

Breeding dispersal — between-season change in breeding location — is usually female-biased in birds and creates problems in accurately estimating annual survival, as conventional Cormack-Jolly-Seber (CJS) survival models cannot discriminate between mortality and undetected emigration. Recently, spatial CJS (s-CJS) models have been developed that use data on breeding dispersal within a population to account for undetected emigration and provide corrected estimates of true annual survival, a development that promises to advance avian conservation initiatives that require accurate estimates of annual survival. Using a 14-year dataset on a color-banded population of Hooded Warblers (Setophaga citrina) in northwest Pennsylvania, I examined female-biased breeding dispersal and performance of a s-CJS model in estimating true annual survival of females and males. I also compared my findings to published literature on other migratory North American warblers, a group with many species of high conservation concern. Breeding dispersal in the Hooded Warbler study population is strongly female-biased, with median dispersal distances of 151 m for females (n = 227) and 51 m (n = 336) for males. Although most individuals disperse short distances, the observed pattern of breeding dispersal within the study site was best modeled using a heavy-tailed Cauchy dispersal kernel, a model that indicates the presence of a substantial tail of undetected long-distance dispersal, particularly in females. Using the Cauchy model, s-CJS analysis yielded realistic estimates of Hooded Warbler true annual survival, 0.61 for both sexes, and resolved ambiguities evident in much lower estimates of apparent annual survival, 0.45 for females and 0.54 for males, derived from conventional CJS analysis. Because long-distance breeding dispersal is widespread in migratory warblers and especially frequent in females, analyses of warbler survival should employ s-CJS methods whenever possible, as estimates of apparent annual survival derived from conventional CJS methods will be poor under-approximations of true annual survival.

Methods
Field Methods, Mapping, and Breeding Dispersal Distance

The study was conducted at Hemlock Hill Field Station in northwest Pennsylvania, USA (41.8°N, 79.9°W). Hemlock Hill consists of 173 ha of primarily beech-maple-hemlock forest located within a rural landscape of forest fragments, agricultural land, and abandoned fields; it has been the site of a long-term field study of Hooded Warblers since 2010 (Mumme 2018, Lignac and Mumme 2023, Mumme 2023). Intensive efforts are made each year to capture and color band all resident breeding adults and to census the entire study population. Data for the current study come from 340 females and 333 males that were color-banded on the study site 2010–2022 and re-sighted or recaptured there 2011–2023. 

Each year I used GPS to map the location, recorded as decimal latitude and longitude, of each resident breeding adult, based on the mean of its nest locations and the location that I judged to be the approximate center of its breeding territory based on resighting locations; including the approximate center of the breeding territory was necessary for mapping the locations of unmated males and a small number of mated pairs for which no nests were found. Seven females changed mates and breeding territories in the middle of a given nesting season, and for these females I mapped their location using their initial breeding territory and nests initiated with their original mate. Total sample size was 1236 location-years, 567 for females and 669 for males (Figure 1).

For color-banded breeding females and males that returned to the study site in a subsequent year, I calculated their breeding dispersal distance as the distance in meters between their locations in year x and year x + 1 (214 females, 326 males), or in a few cases where individuals were absent 1–2 years before returning, year x + 2 (12 females, 9 males) or year x + 3 (1 female, 1 male); inclusion of individuals absent 1–2 years is unlikely to alter the results, as their breeding dispersal distances did not differ from those of birds returning in consecutive years, for both females (t = 0.3, df = 225, P = 0.73) and males (t = 1.5, df = 334, P = 0.13). The total sample size available for analysis was 227 breeding dispersal events for females and 336 for males. I tested for sex differences in breeding dispersal distance using a general linear mixed model (GLMM) with an identity link function in which sex was considered a fixed effect and an individual’s unique identification code was entered as a random effect to account for multiple measures of breeding dispersal distance obtained from the same individual. To normalize the data, breeding dispersal distances were log-transformed prior to analysis. Statistical analysis was conducted using JMP Pro 17.0 (SAS Institute, Cary, North Carolina, USA).

Identification of Appropriate Dispersal Model

I used the simulation capabilities of the spatial CJS (s-CJS) model of Schaub and Royle (2014) to identify the most appropriate dispersal kernel with which to model Hooded Warbler breeding dispersal. Schaub and Royle (2014) model dispersal as change in location along both the x and y axes in a two-dimensional equidistant grid. For clarity, I refer to these changes as x-y dispersal distances to distinguish them from breeding dispersal distance described in the previous section; while breeding dispersal distances are always positive, x-y dispersal distances can include negative values, reflecting dispersal movements that result in decreases in either longitude (x), latitude (y), or both. Under most circumstances, the distributions of x-y dispersal distance along each axis would be expected to be symmetrical and centered around 0, and can be modeled using distributions with these characteristics (Schaub and Royle 2014).

I first represented the 173 ha of the study area within a 20 x 18 equidistant x-y grid where each grid unit = 100 m. I then replotted the 1236 bird location-years according to this grid system (Figure 1) so that the 563 observed breeding dispersal distances (227 for females, 336 for males) could also be expressed as x-y dispersal distances. I then used the simulation capabilities of the model of Schaub and Royle (2014) to simulate dispersal within and outside of this gridded study area (Figure 1) using three different theoretical dispersal kernels: (1) a random normal distribution, (2) a t distribution with df = 2, and (3) a Cauchy distribution. The three models can also be viewed as three different t distributions with decreasing degrees of freedom and increasingly heavy tails, as the normal and Cauchy distributions are equivalent to t distributions with df = ∞ and df = 1, respectively.  Although the model of Schaub and Royle (2014) allows dispersal kernels on the two axes to have different values of the scale parameter σ (σ1 and σ2 for the x and y axes, respectively), the observed distribution of x and y dispersal distances within the Hooded Warbler study area did not differ significantly for either females (Kolmogorov-Smirnov D = 0.070, P = 0.63) or males (Kolmogorov-Smirnov D = 0.065, P = 0.47). I therefore set σ1 = σ2 = σ in all dispersal simulations. 

Using the three dispersal kernel models, I then generated simulated x-y dispersal data that could be compared to observed x-y dispersal data for each sex, searching for the combination of kernel model and σ that would best fit the observed x-y dispersal data, as determined by the Kolmogorov-Smirnov two-sample test. Although each simulation was designed to generate about 6,600 total simulated dispersal events, the number of simulated x-y dispersal events that occurred within the study area — the ones that could be compared to observed values — varied with both dispersal kernel and σ. The simulations (n = 56) focused on regions of σ-space that minimized the Kolmogorov-Smirnov D statistic, and maximized the corresponding P-value, for each of the three tested dispersal kernels. 

Estimates of Detection Probability and Annual Survival

I estimated detection probability (p) and apparent annual survival (ɸ) via conventional CJS analysis using R version 4.3.2 (R Core Team 2023) and package marked version 1.28 (Laake et al. 2013). For this analysis, histories of individual females and males were combined into a single dataset and analyzed with a 4-parameter model with separate estimates of p and ɸ for each sex. I then used the Bayesian spatial CJS (s-CJS) model of Schaub and Royle (2014) to estimate the same parameters correcting for undetected emigration. The s-CJS model was run using JAGS version 4.3.2 (Plummer 2003) and the R package r2jags 0.7-1 (Su and Yajima 2021). For the s-CJS analysis I used the dispersal kernel model identified in the previous section and ran the sexes separately. Vague priors (random uniform distribution between 0 and 1) were used for p and ɸ. Because preliminary models encountered convergence issues in estimating the critical dispersal scale parameters σ1 and σ2, particularly in females, more specific priors — random uniform distribution between 0.001 and 2 — were employed for σ1 and σ2 in the final models; these more specific priors are justified, as the simulations described in the previous section indicated that their posterior estimated values would likely fall comfortably within that range. The Markov chain Monte Carlo (MCMC) settings specified three chains, each with 420,000 iterations with the first 20,000 discarded, for a total of 1.2 million saved iterations. No convergence issues were encountered in the final models, with all Rhat values less than 1.02 for females and 1.01 for males.