Heterogeneity in song rates in the Collared Flycatcher (Ficedula albicollis) explained with the availability parameter in generalized N-mixture models: Its importance for abundance estimates in avian aural counts

1.   Introduction

The ease of obtaining data originating from repeated counts of animals and the possibility to get abundance estimates unaffected by imperfect detection makes N-mixture models (Royle, 2004) increasingly popular (Kéry and Royle, 2016). The paper, which introduced the model (Royle, 2004) has been cited 928 times based on CrossRef accessed via the Biometrics journal website (as of December 7, 2022), which indicates a high importance of this type of model in studies of abundance. However, their sensitivity to violations of assumptions has been repeatedly stressed (Hunt et al., 2012; Warren et al., 2013; Barker et al., 2018; Link et al., 2018) but their validations were reported equally often (Bötsch et al., 2019; Costa et al., 2019, 2021; Kéry, 2018; Neubauer et al., 2022). Among a few important assumptions, these models assume that detections are independent among individuals: detecting one individual at a site is assumed not to affect detecting another one (Royle, 2004). This assumption can be violated when song rates increase with abundance N_i at a site i (Kéry and Royle, 2016), where territorial males motivate each other to vocalize as demonstrated by Warren et al. (2013) and Neubauer et al. (2022). Within N-mixture models, unmodelled heterogeneity among individual birds in detection probability p should lead to underestimating abundance, so if this relationship is present, but unaccounted for, abundance will be biased low and thus, conclusions remain questionable. Simulations by Kéry and Royle (2016, pages 248–250), however, did not support this expectation: both parameters, detection probability p and abundance λ were estimated approximately correctly in the presence of individual heterogeneity in p, even if this was ignored in the model. One of the approaches suggested to overcome this type of heterogeneity in p is to model it as a function of abundance (if available) at a site N_i or a measure of it (like the maximum per-site count), despite some form of circularity probably being involved; few (if any) studies addressed this issue (Kéry and Royle, 2016). In this paper we expand this exercise, by taking advantage of generalized binomial N-mixture models, where detection is decomposed into availability probability and actual detection probability, provided an individual is available (Chandler et al., 2011). Song rates can then be directly captured by the availability parameter φ, which is separate from detection probability p. In contrast, both these processes are involved and summarized in a single 'detection probability' parameter in the simpler, standard binomial N-mixture model (Royle, 2004). We fit generalized binomial N-mixture models to the repeated counts, collected with a robust design protocol, of a forest passerine, the Collared Flycatcher (Ficedula albicollis) (Neubauer and Sikora, 2013). We tested whether (1) song rates increase with observed abundance, (2) adding an explanatory variable for φ in the form of observed count improves model fit and inference and (3) is there any effect on abundance estimates if this heterogeneity is not accounted for (i.e., a misspecified model is fitted). We also fitted a standard, binomial N-mixture model (Royle, 2004), which ignores the multiscale nature of these data to compare abundance estimates when availability is conflated with detection. We then explored the results to strengthen understanding and outline profits resulting from modelling density-dependence in hierarchical models for abundance.

2.   Methods

2.1   Study site and field methods

We used a dataset from Collared Flycatcher repeated counts in Las Warmiński nature reserve (N Poland), originally analysed with dynamic occupancy models (Neubauer and Sikora, 2013), except that binary detection-nondetection data were converted back to the original abundance data (per-site and per-visit number of singing males). The reserve, established in 1982, covers 1803 ​ha, including 1550 ​ha of forests (Fig. 1). This forest is mostly coniferous (with Scots Pine Pinus sylvestris dominating over Norway Spruce Picea abies), but isolated patches of deciduous stands, mostly with Pendunculate Oak (Quercus robur) (up to 240 years old), European Hornbeam (up to 180 years old) and Black Alder (Alnus glutinosa) along the waterbodies and rivers provide suitable breeding habitat spots for the Collared Flycatchers. Seventy-eight sites, being State Forests sections, with boundaries frequently drawn along forest roads, streams or lake shores, were highly variable in shape and size (0.5–11 ​ha). They shared suitable habitat (dominant tree age of 100+ years, deciduous trees prevailing) and were predefined prior to the study. Observers (the authors) surveyed them with standardized effort (6–8 ​min spent on counting singing birds per hectare) to detect the number of the Collared Flycatcher males present. Four individual visits (hereafter counts) to each site were performed between late April and late May 2012. The robust design approach was followed: counts were repeated close in time to form 'surveys', over which the population is assumed closed. In this study, two counts formed survey 1 (all 78 sites visited with 1–3 day intervals: during 27–30 April and again, during 30 April–3 May), and after ca. two weeks, the next two counts formed survey 2 (all 78 sites visited during 18–20 May and again, during 20–23 May, keeping intervals between visits to the same sites within the range of 1–4 days). Surveys 1 and 2 were separated by longer intervals (ca. two weeks), over which the population was assumed open. Sites were visited to count flycatchers between sunrise (05:00 during survey 1, 04:30 during survey 2) until afternoon, and usually ended around 13:00–14:00. For further details on fieldwork, see Neubauer and Sikora (2013).

Figure  1.  Schematic map of the study area, Las Warmiński nature reserve. Surveyed sites are depicted by grey polygons (their boundaries are approximate), while the distribution and maximum per-site counts over four visits (C_maxi) of Collared Flycatcher males is shown with circles varying in size for each site. Łańsk (a resort) and Cyranka, marked with dashed, red lines are the only permanent settlements within the reserve. Forests in green, white depicts deforested areas (mainly marshes and meadows). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

DownLoad: Full-Size Img PowerPoint

2.2   Analysis of count data: models and parameters

We used the generalized, binomial N-mixture model developed to account for geographical openness by considering temporary emigration and estimating it with a parameter called availability probability φ (φ ​= ​1 – temporary emigration; Chandler et al., 2011) in addition to abundance λ and detection probability p. The temporary emigration can be any process which makes an individual periodically unavailable for detection. The φ parameter emerges from decomposing the processes responsible for detecting an individual into the probability of being available (p_a, sensu Nichols et al., 2009) and the "true" probability of detecting an individual (p_d) provided it is available, i.e., it sings (so that the probability of detecting an individual is a product p_a ​× ​p_d, see Nichols et al., 2009). The availability parameter can adequately address either of the two types of temporary emigration leading to the individual being unavailable for detection: its movements in space resulting in leaving the site (but not necessarily the territory) physically or being inactive (like staying silent or being hidden) despite being present at a site. These two processes were named "spatial" and "random", respectively, by Kéry and Royle (2016, page 485), and we follow these names here. Despite being relatively intuitive, depending on which of the two processes it addresses in a given case, availability might or might not be incorporated into density estimates, the quantity for which hierarchical models were developed (Kéry and Royle, 2016). If temporary emigration is due to spatial movements, then the density under the model is D = (λ ​× ​φ)/area. Not accounting for temporary emigration leads to overestimation of abundance, since entire territories are 'included' into the area surveyed (a 'site'), despite their parts extending (well) beyond the boundaries of sampling sites (see Chandler et al., 2011). This overestimation can be particularly strong with point counts, but including availability or applying distance sampling protocols or both help to solve this problem (Neubauer and Sikora, 2020). With individuals staying silent (like not vocalising in birds, hiding or burrowing in other animal groups), density is D ​= ​λ/area, since individuals go undetected due to being unavailable for detection, despite being present. However, if availability is included into density estimations in such cases, the estimated quantity is the number of available individuals (e.g., visible or vocally active, that may potentially be detected). In the simpler binomial N-mixture model (Royle, 2004), a parameter named detection probability summarizes both p_a and p_d.

Generalized, binomial N-mixture models have a hierarchical structure, composed of three submodels: for estimating abundance, availability probability and detection probability, given an individual is available for detection (Chandler et al., 2011). We fitted 16 such models (see Table 1) to count data of the Collared Flycatcher using the gpcount function. The abundance submodel had its parameter (abundance) constant but included log (site area) as an offset, which makes the model effectively estimating density per area unit (ha). The submodels for availability differed, to make possible direct comparisons of models considering density dependence and ignoring it. We used observed, maximum count C_max at a site as a covariate for availability, being a factor with four levels: no male observed, a single male observed, two males observed and three to five males observed, with the latter category pooling the very few sites where 3 or 5 males were detected (no sites held 4 males). The availability probability was thus either (1) kept constant (a model which does not account for density-dependent probability of vocalising), (2) modelled as survey-dependent (as above, but survey dependence allows to address within-season variation, i.e., differences in availability between surveys 1 and 2), (3) modelled as a function of maximum count C_max at a site (a model which addresses density-dependence in the probability of vocalising) or (4) modelled as both survey- and maximum count-dependent (a model which addresses variation from (2) and (3)). The detection probability (p_d sensu Nichols et al., 2009, provided that an individual is available) in our models was modelled with another set of covariates, which included observer effect (a factor with two levels: G.N. or A.S.) and time of day (expressed in minutes since sunrise, see above; the latter being a survey-specific covariate and cannot be treated as a covariate for availability). Therefore, the models forming our set could be ranked from the simplest to the most complex, given the structure of submodels for availability and detection. The most complex model (model (4) in Table 1) can be described as:

Table  1.  Generalized binomial N-mixture models fitted to repeated count data of Collared Flycatcher (Ficedula albicollis) in Las Warmiński nature reserve, N Poland, spring 2012. The three Greek letters denote submodels for abundance (λ), availability probability (φ) and detection probability (p). Heterogeneity in availability probability (i.e., probability of vocalising or song rate) was either unmodelled (with availability probability kept constant across sites and unrelated to abundance, denoted φ) or modelled as a function of maximum count at a site (denoted φ_Cmax) or survey (φ_survey) or both. The remaining submodels were either the same across models (abundance, λ, modelled as constant across sites but with site area included as an offset) or modelled as a function of observer (p_obs), time of the day, (p_time) or both (detection probability, p). Models are ranked according to AIC, from top-supported (the lowest AIC value) to least supported. N_p is the number of parameters in the model, '+' denotes models with weights < 0.01.

Model structure	np	AIC	ΔAIC	ωAIC	Cumulative ω
(1) { λ, φsurvey + Cmax, p }	7	499.83	0.00	0.37	0.37
(2) { λ, φsurvey + Cmax, ptime }	8	500.80	0.97	0.23	0.60
(3) { λ, φsurvey + Cmax, pobs }	8	501.83	1.99	0.14	0.74
(4) { λ, φsurvey + Cmax, pobs + time }	9	502.62	2.79	0.09	0.83
(5) { λ, φCmax, p }	6	503.19	3.35	0.07	0.90
(6) { λ, φCmax, ptime }	7	503.57	3.74	0.05	0.95
(7) { λ, φCmax, pobs }	7	505.18	5.35	0.03	0.95
(8) { λ, φCmax, pobs + time }	8	505.47	5.64	0.02	0.97
(9) { λ, φsurvey, p }	4	554.78	54.95	+	1.00
(10) { λ, φsurvey, ptime }	5	554.89	55.05	+	1.00
(11) { λ, φsurvey, pobs }	5	556.77	56.94	+	1.00
(12) { λ, φsurvey, pobs + time }	6	556.81	56.97	+	1.00
(13) { λ, φ, ptime }	4	558.26	58.43	+	1.00
(14) { λ, φ, p }	3	558.82	58.99	+	1.00
(15) { λ, φ, pobs + time }	5	560.12	60.29	+	1.00
(16) { λ, φ, pobs }	4	560.79	60.96	+	1.00

 | Show Table

DownLoad: CSV

where M_i denotes the super population size, i.e., the total number of individuals potentially to be detected, N_i denotes the number of individuals available for detection given φ and C_i is the number of individuals actually detected, given p. Any simpler model can be obtained by removing the effects to be omitted.

We focused on paired comparisons of fitted models within our setting by comparing two models each time: one with density dependence in availability accounted for and the other one with this form of heterogeneity unmodelled and identical structure of the remaining submodels. This was based on information criteria (AIC). Model fitting and model selection was performed with the R (version 3.6.3, R Core Team, 2021) unmarked package (Fiske and Chandler, 2011) using the gpcount (Chandler et al., 2011) and pcount (for standard, binomial N-mixture model) functions. Models were ranked with the modSel function applied to an unmarked FitList object.

To get the estimates of density (derived quantity, computed as λ ​× ​φ) along with its variance, we applied a nonparboot function from the unmarked package. It uses non-parametric bootstrapping for variance estimation under single models. With this technique, the original dataset is resampled with replacement many times, a given model fitted to each resample and parameters then estimated. We repeated bootstrapping 100 times for the three top-supported models to get approximate distributions of parameters; means with SD and confidence intervals (Table 2, Figs. 1 and 2) were computed from these distributions.

Table  2.  Parameter estimates from the generalized N-mixture models fitted to count data of Collared Flycatcher (Ficedula albicollis) in Las Warmiński nature reserve, 2012 (N Poland). Models (1), (2) and (3) are the top-supported models from Table 1. All the values are on a normal scale with their 95% confidence intervals in parentheses. Availability probability is shown for the sites occupied by 1, 2 and 3 to 5 males. Estimates of the total population size (summed abundances across all sites surveyed) for the surveys 1 and 2 are also shown.

Model	Abundance	Availability			Detection	Population size
		1 male	2 males	3–5 males
Model (1)
{ λ, φsurvey + Cmax, p }	1.696 (0.730–2.34)		Survey 1:		0.327 (0.196–0.462)	survey 1:
		0.190 (0.091–0.535)	0.306 (0.140–0.892)	0.394 (0.169–0.910)		101 (65–185)
		Survey 2:				survey 2:
		0.279 (0.157–0.569)	0.421 (0.218–0.999)	0.529 (0.277–1.000)		148 (93–222)
Model (2)
{ λ, φsurvey + Cmax, ptime }	1.537 (0.510–2.252)		Survey 1:		0.373 (0.212–0.628)a	survey 1:
		0.240 (0.109–0.719)	0.346 (0.149–0.906)	0.413 (0.197–0.967)		106 (60–180)
			Survey 2:			survey 2:
		0.334 (0.162–0.999)	0.460 (0.209–1.000)	0.564 (0.315–0.999)		148 (86–224)
Model (3)
{ λ, φsurvey + Cmax, pobs }	1.760 (0.720–2.284)		Survey 1:		A.S.: 0.344 (0.193–0.530) G.N.: 0.339 (0.208–0.482)b	survey 1:
		0.178 (0.083–0.477)	0.271 (0.115–0.702)	0.366 (0.189–0.798)		97 (65–148)
			Survey 2:			survey 2:
		0.270 (0.138–0.582)	0.386 (0.190–0.897)	0.501 (0.297–0.887)		148 (95–239)
a Only the intercept is shown, despite this model including insignificant time effect in detection. b Estimates for both observers (A.S. and G.N.) are shown.

 | Show Table

DownLoad: CSV

Figure  2.  Relationships between the song rate (i.e., probability of vocalising) measured with availability probability in the generalized binomial N-mixture models and maximum count (C_maxi) at a site for the Collared Flycatcher repeated surveys in the Las Warmiński reserve, N Poland, spring 2012. Means (points) ​± ​1 SD (whiskers) are shown. The three symbols for each C_maxi indicate estimates under the three top-supported models in Table 1 (models 1, 2 and 3, with AIC weights of 0.42, drawn with dark grey, 0.17, drawn with medium grey, and 0.16, drawn with pale grey). Estimates of availability probability from the models with identical structure except for ignoring relationships between availability and C_maxi (thus only a single availability value, equal for all the sites, regardless of maximum count is estimated under these models) are shown on the right with empty symbols. These three models (Table 1: models 9, 10 and 11), had AIC weights of zero.

DownLoad: Full-Size Img PowerPoint

3.   Results

There was overwhelming evidence that models with availability modelled as a function of observed, maximal per-site abundance C_max were much better supported than the ones where this heterogeneity was not accounted for. The cumulative AIC weight of the former group of models was 1 (Table 1). Models with identical structure as the three top-supported models, but missing this single effect in φ, were worse in terms of support by at least 50 AIC units. The top-supported models (with ΔAIC ≤ 2) agreed that the availability probability of the Collared Flycatcher males increased with maximum observed count (a measure of local abundance), suggesting that per capita song rates were higher where there were more males at a site. For sites with a single male, two males and 3–5 males, model-averaged estimates from the three top-supported models were 0.20 (95% confidence intervals: 0.13–0.40), 0.31 (0.19–0.56) and 0.40 (0.24–0.64), during survey 1 and 0.30 (0.20–0.52), 0.43 (0.27–0.65) and 0.54 (0.37–0.75) during survey 2, respectively (Fig. 2). If this relationship was ignored in a model, a single availability value for all the sites was estimated higher: mean values were 0.48 for survey 1 and 0.68 for survey 2 (Fig. 2).

Abundance estimates for the whole studied area (i.e., the total population size for all the sites surveyed) were very similar for surveys 1 and 2 under the three top-supported models (97–106 males during survey and 148 during survey 2), with wide confidence intervals (Table 2). Analoguous estimates from the three identical models, missing the effect of maximum count in φ but with otherwise identical structure were respectably lower, particularly for survey 2 (Fig. 3). This confirms that unmodelled heterogeneity in availability leads to an underestimation of abundance: estimated population sizes were lower by 22–27% under the models ignoring heterogeneity in availability and the total population size was estimated at 76–77 males during survey 1 and at 108–112 males during survey 2 (Fig. 3). The total population size estimated under the standard binomial N-mixture model fitted to all four counts (thus ignoring the multiscale nature of the data and with detection probability kept constant) was 142 (95% CI: 93–227) males, which is very similar to the survey 2 estimates under generalized models (Table 2).

Figure  3.  Total population sizes of the Collared Flycatcher in the Las Warmiński reserve, N Poland, spring 2012, estimated under generalized binomial N-mixture models. Means (points) ​± ​1 SD (whiskers) are shown. Models (1), (2) and (3) drawn with filled symbols are the top-supported models (Table 1: models 1, 2 and 3) and include the relationship between availability and maximum count at a site (C_maxi), while 'paired' models (9), (10) and (11) are shown with empty symbols and have identical structure as (1), (2) and (3) except that they ignore this relationship.

DownLoad: Full-Size Img PowerPoint

Detection probabilities were rather low, ranging between 0.3 and 0.4 under the three top-supported models (Table 2).

4.   Discussion

This study is most likely the first one based on real-world data, which explores the relationship between availability probability and local abundance with generalized, binomial N-mixture models. We show that (1) per capita song rates were dependent on a measure of local abundance – the maximum observed count C_max and that (2) modelling availability as a function of C_max helps to explain this heterogeneity and improves inference. Previously, male-male vocal interactions known as counter singing, were found to increase song rates and affect detection probability with structurally simpler, standard binomial N-mixture models for the Marsh Tit (Poecile palustris), another forest passerine (Neubauer et al., 2022). Similarly, Warren et al. (2013) demonstrated that per-capita song rates increased with abundance in Golden-cheeked Warblers (Setophaga chrysoparia), violating the assumption of independence in detection probability in binomial N-mixture models. This resulted in abundance estimates biased low, and the strength of this bias increased with local abundance. However, in standard binomial N-mixture models, detection probability includes availability and is a product of both processes, while in the generalized version we used here, these are separate parameters. It is therefore natural that if availability reflects the probability of vocalising, the independence assumption refers to being vocally active, not to detection itself (unlike in N-mixture models [Royle, 2004], where detection probability summarizes both availability and actual detection). In the present work, we showed that, despite some form of circularity being involved (as already stressed by Kéry and Royle, 2016), models accounting for this form of heterogeneity – i.e., song rates dependent on local abundance – were much better supported than models ignoring it. Therefore, using observed, maximum per site counts as a covariate for φ improves model fit and helps to explain this variation. If ignored, it introduces a negative bias in abundance estimates with binomial N-mixture models as well (Warren et al., 2013). At the other hand, structurally the simplest binomial N-mixture model, when fitted to our flycatcher data, produced abundance estimates similar as the generalized models that consider relationship between C_max and availability (Table 2). It remains to be studied whether this is just a case or a regular feature of binomial N-mixture models with a simple structure.

Kéry and Royle (2016) mention that in WinBUGS it is also possible to specify a binomial N-mixture model, where the latent abundance N_i (not a proxy for it, as we did in this study using C_max) is included as a covariate for detection. Some such models experimentally fitted to simulated data seemed estimable, but others – fitted to real-world data – failed to produce promising estimates (Kéry and Royle, 2016, page 310).

Previous analysis of the count data used here (Neubauer and Sikora, 2013), applied dynamic occupancy models (MacKenzie et al., 2003). In their simplest form, as used in Neubauer and Sikora (2013), these models can only estimate detection probability (as a product of being available and being detected), and we've used maximum count C_max as a covariate for detection probability, to show that the latter increases with abundance (more precisely, the number of neighbours, that is C_max – 1 in an attempt to reflect habitat quality). However, in occupancy models, the estimated detection probability refers to a species (i.e., the probability of detecting at least one individual at a site). With detection decomposed into availability probability and actual detection probability in the generalized N-mixture models (Chandler et al., 2011) we used here, one can directly address the probability of vocalising (song rate) via the former parameter, which makes these models likely more appropriate to estimate abundance than binomial N-mixture models, particularly when the field protocol is based on point counts (see Neubauer and Sikora, 2020).

In summary, the so far little explored relationship between local abundance (or a measure of it) and detection (in standard binomial N-mixture models, Royle, 2004) or availability probability (in generalized binomial N-mixture models, Chandler et al., 2011) clearly deserves further attention. Using local, observed counts, or density, as a covariate for availability or detection should improve model fit and resulting abundance estimates (Kéry and Royle, 2016, page 310) with any commonly applied field protocols, like point counts. It remains to be verified with both simulations and real-world data.

Author contributions

GN and AS formulated the question, designed field work, collected data and wrote the paper. GN analyzed the data. Both authors read and approved the final manuscript.

Ethics statement

The study was conducted in compliance with permissions from the Regional Directorate for Environmental Protection, Olsztyn, Poland.

Declaration of competing interest

The authors declare that they have no competing interest.

Acknowledgements

The authors wish to express their gratitude to Marc Kéry and an anonymous Reviewer for their constructive suggestions, extremely helpful in improving the manuscript. Our thanks go also to Daniel O'Connell for improving the language, and to Tomasz Chodkiewicz and Zenon Rohde for their help with map preparation.