Point Count Data

The Boreal Avian Modelling Project (BAM) was initiated to facilitate the conservation of boreal birds. As part of these efforts, BAM has collated data from georeferenced point count surveys conducted in boreal, hemiboreal, and sub-Arctic biomes by independent inventory, monitoring, and research studies (Cumming et al. 2010, Barker et al. 2015). The BAM database (version 5) contains 250,822 off-road point count surveys conducted at 145,289 point count locations by 131 separate studies between 1991 and 2014. Of these surveys, 3% were 3 min, 60% were 5 min, and 37% were 10 min in total point count duration. We applied removal models to a subset of 89,304 off-road point counts at 25,754 locations conducted by 71 projects in boreal and hemiboreal regions (hereafter, ‘boreal') in which observers recorded when they first detected each bird relative to ≥2 time intervals (hereafter referred to as the ‘full dataset'; Figure 1). We excluded roadside surveys because of possible differences in singing behavior between birds breeding along roadsides vs. in off-road areas (Kociolek et al. 2011). The number of time intervals varied from 2 to 10, and there were 11 distinct combinations of count duration and the number and length of the time intervals (Figure 1). We restricted our analyses to detections of vocalizing gallinaceous gamebirds, shorebirds, woodpeckers, and songbirds, which included species that are diurnally active and are first detected primarily by their auditory displays during breeding surveys. We analyzed species detected aurally during a minimum of 75 point count surveys in the full dataset (median = 1,240 point count surveys, maximum = 31,334). We present results for the 152 species for which we could successfully estimate the constant singing rate parameter of the conventional removal model (Appendix Tables 2 and 3).

Count duration varied from 3 min to 20 min, but we omitted time intervals beyond 10 min to minimize the inclusion of movements by birds during counts, which could violate the closed population assumption of the removal model. For example, individuals that move toward the observer are more likely to be detected than ones that move away, and would inflate counts (Glennie et al. 2015). Point count radii in our dataset were at least 50 m, with 100-m radius (24%) and unlimited distance (72%) counts making up the majority. The removal model captures availability, while perceptibility is often estimated via distance sampling (Diefenbach et al. 2007). The conditional model formulation uses the cumulative counts over the time intervals relative to the total count, which makes the estimation of availability independent of the observer and distance-related perceptibility component and of the true underlying population abundance. By conditioning on the total counts, the abundance and other nuisance variables (such as perceptibility) are canceled out of the model likelihood (Sólymos et al. 2013).

As a result of this conditional modeling approach, we treated within- and between-year repeat visits to the same locations as independent observations. We checked that our assumption of conditional independence was valid, and we found neither systematic bias in our parameter estimates nor shrinking of standard errors as a result of treating data from revisits as independent (Appendix Figure 6). Using data from within- and between-year revisits to the same locations allowed us to better estimate date- and time-related patterns in our dataset, because revisits increased our sampling coverage of different dates and times, and therefore strengthened our inference regarding longitudinal variation. The diurnal and seasonal timing of the point counts varied widely, so we standardized the timing of each point count relative to local sunrise and spring, respectively, using the National Oceanic and Atmospheric Administration's Sunrise/Sunset Calculator ( https://www.esrl.noaa.gov/gmd/grad/solcalc/sunrise.html, Meeus 1999), implemented in the R package maptools (Bivand and Lewin-Koh 2016). We quantified the ordinal day of each point count, and also a standardized version of the ordinal day of each point count relative to the 30-yr average of the local start of the growing season (McKenney et al. 2006). Temporal covariates were standardized by dividing the observed values by the maximum possible values: time since local sunrise (SR; mean = 2.2 hr, range = −5.7 to 12.5 hr, maximum possible = 24 hr); ordinal day of year (DY; mean = 166, range = 128 to 200, maximum possible = 365); and days since local spring (LS; mean = 48.5, range = −8 to 94, maximum possible = 365).

Continuous Time-removal Models

Time-removal models are based on a removal experiment wherein animals are trapped and thereby removed from the closed population of animals being sampled (Zippin 1958). When applying removal models to avian point count surveys, the counts of singing birds (Y_ij, …, Y_iJ) within a given point count i (i = 1, …, n) are tallied relative to when each bird is first detected in multiple and consecutive time intervals, with the survey start time t_i0 = 0, the end times of the time intervals t_ij (j = 1, 2, …, J), and the total count duration of the survey t_iJ. Each individual bird is counted once, so that individuals are ‘mentally removed' from a closed population of undetected birds by the surveyor. In the continuous time-removal model (hereafter, ‘removal model'), singing events by individual birds are assumed to follow a Poisson process. We can use the rate parameter of the Poisson process (ϕ) to estimate the singing rates of birds during a point count (Alldredge et al. 2007a).

In the conventional removal model (M₀), all individuals of a single species at a given location and time are assumed to be homogeneous in their singing rates. The time to first detection follows the exponential distribution f(t_ij) = ϕ exp(−t_ij ϕ), and the cumulative density function of times to first detection in time interval (0, t_iJ) gives us the probability that a bird sings at least once during the point count as p(t_iJ) = 1 − exp(−t_iJ ϕ) (Barker and Sauer 1995, Alldredge et al. 2007a, 2007b, Sólymos et al. 2013). Although the term ‘detection' implies both the perceptibility and detectability processes, the perceptibility component is not part of our models due to the conditional likelihood approach. The perceptibility component should be addressed in a distance sampling framework independent of the removal models (Handel et al. 2009, Sólymos et al. 2013).

In the continuous time formulation of the finite mixture (or 2-point mixture) removal model, the cumulative density function during a point count is given by p(t_iJ) = (1 − c) 1 + c [1 − exp(−t_iJ ϕ)] = 1 − c exp(−t_iJ ϕ), where ϕ is the singing rate for the group of infrequently singing birds, and c is the proportion of birds during the point count that are infrequent singers (see Jewell 1982 for a discussion of mixtures of exponential distributions). The remaining proportions (1 − c; the intercept of the cumulative density function) of the frequent singers are assumed to be detected instantaneously at the start of the first time interval (Alldredge et al. 2007b). In the simplest form of the finite mixture model, the proportion and singing rate of birds that sing infrequently is homogeneous across all times and locations (model M_f).

The singing rate in the conventional removal model (ϕ_i) may also vary with respect to survey-specific covariates. We hereafter refer to these models collectively as time-varying conventional removal models (M₀^ϕ). Previously, researchers have applied covariate effects to the parameter ϕ_i of the finite mixture model, assuming that the parameter c is constant irrespective of time and location (i.e. only the infrequent singer group changes its singing behavior; Farnsworth et al. 2002, Aldredge et al. 2007a, Etterson et al. 2009, Reidy et al. 2011). An alternative parameterization is for c_i rather than ϕ to be the time-varying parameter, allowing individuals to switch between the frequent and infrequent singer group depending on covariates. We defined 2 time-varying finite mixture model classes: (1) M_f^ϕ, with covariate dependent rate parameter p(t_iJ) = 1 − c exp(−t_iJ ϕ_i); and (2) M_f^c, with the proportion of infrequent singers dependent on covariates: p(t_iJ) = 1 − c_i exp(−t_iJ ϕ).

In all 3 classes of time-varying models (M₀^ϕ, M_f^ϕ, and M_f^c), we modeled covariate effects as linear or quadratic effects of time since local sunrise (SR), ordinal day of year (DY), and days since local spring (LS) in 14 different candidate models (Table 1). We estimated model parameters using the conditional maximum-likelihood procedure described by Sólymos et al. (2013) and implemented in the R package detect (R Core Team 2017, Sólymos et al. 2018). This included using the cmulti function, with type argument ‘rem' for conventional (M₀, M₀^ϕ) models and type argument ‘fmix' and ‘mix' for the 2 mixture model types (M_f^ϕ and M_f^c, respectively; both options give identical results for M_f). Covariate effects for ϕ_i were modeled using the logarithmic link function to ensure that ϕ_i > 0, whereas covariate effects for c_i were modeled using the logistic link function to ensure that c_i was in the interval (0, 1).

Sample Size Requirements for Removal Models

We first examined the minimum sample size required to reliably estimate the parameters of the time-invariant M₀ and M_f models. We stratified the full point count dataset by project, and then estimated the model parameters for 3,002 combinations of the 65 projects and 152 species that had a minimum of 2 nonzero total counts (hereafter referred to as the ‘project-level dataset'). We classified an estimation as a success if the Nelder-Mead optimization algorithm for the conditional maximum-likelihood estimator (MLE) was able to converge on a point estimate of ϕ for M₀ or on a joint estimate of (ϕ, c) for M_f, and if the algorithm could estimate the corresponding variances, Var[log(ϕ)] for M₀, and jointly Var[log(ϕ)] and Var[logit(c)] for M_f. We fit logistic regression models separately to the outcomes of estimates and variances for M₀ and M_f models to predict the probability of successfully estimating the model parameters as a function of log number of point counts with nonzero totals for the species in a given project, and calibrated the fitted models to predict the minimum number of nonzero counts required to have a 90% probability of successfully estimating removal model parameters and variances. The number of nonzero counts was highly correlated with the mean count (Spearman's r = 0.82, P < 0.001) in the full dataset; therefore, we did not separately assess relationships with mean counts.

Adjusting Point Counts for Variation in Count Duration

We used a subset of the full dataset for 133 species from 38,540 point count surveys with 0–3, 3–5, and 5–10 min intervals (hereafter referred to as the ‘validation dataset') to assess the bias and variance based on the time-invariant models (M₀ and M_f; estimates based on the validation data) when predicting the number of birds counted in the full 10-min point count, using the data from the 0–3 min and 0–5 min portions of the same point counts as inputs. The expected total number of individuals of a species counted across point counts of duration t can be written as E[Y_t] = N p_t, where p_t is the estimated availability probability for duration t from the removal model (M₀ or M_f), and N is the unobserved total number of individuals of that species present across the point counts. It follows that N = E[Y₃]/p₃ = E[Y₅]/p₅ = E[Y₁₀]/p₁₀ for intervals of the same point count, in which case N is assumed to be constant due to the closed population. We computed the bias of the removal model when adjusting the 3-min and 5-min point counts to a 10-min standard, respectively, as (ΣY₃/p₃) / (ΣY₁₀/p₁₀) and (ΣY₅/p₅) / (ΣY₁₀/p₁₀), with a value of 1 resulting when the adjustment was unbiased (hereafter referred to as the ‘corrected relative count'). For each species, we fit M₀ or M_f models and drew B = 1,000 replications from the asymptotic distribution of the parameter estimates to represent uncertainty, and we calculated bias (mean signed difference from the expected value of 1) and variance (mean squared deviation of the corrected relative counts from their sample mean across the B replications) based on corrected relative counts from the B replicates. We graphically displayed the distribution of bias and variance across species as a function of sample size, count duration (3 or 5 min), and model type (M₀ or M_f). Bias and variance are presented as proportions of the 10-min corrected count. We used the median and 5% and 95% quantiles to represent the distribution of bias and the 5–95% interval to represent the distribution of variance of the species whose sample sizes were equal to or greater than the sample size selected for comparison.

Relative Support for Models and Plotting Time-varying Effects

For each of the 152 species, we fit 44 candidate models that included the time-invariant (M₀ and M_f) and time-varying models with temporal covariate effects (M₀^ϕ, M_f^ϕ, and M_f^c; 14 models for each) based on the full dataset (Table 1, Appendix Table 3). Although our dataset was composed of many point counts, the sample size for conditional likelihood estimation is defined by the number of point counts with >0 total counts. As a result, sample sizes varied considerably across species. Therefore, we compared the relative support for each model within a species using Akaike's information criterion corrected for small sample sizes (AIC_c). We assessed model support across species by counting the number of species for which each model was best-supported (smallest AIC_c value). We only compared models for the 147 species for which all model parameters were estimated. We compiled information on each species' migratory strategy (winter resident vs. migrant) from Poole (2012), and used a chi-square test to see whether migratory strategy affected how often ordinal day (DY) and days since local spring (LS) were selected in the best-fitting models (Appendix Table 3).

We examined temporal patterns in availability probabilities across the species for which the best-fitting removal model included temporal covariates (SR, DY, or LS; Table 1, Appendix Table 3). We did so by predicting availability probabilities based on a 3-min point count duration for each species as a function of the observed covariate values from the full dataset (response curves). We conditioned the predictions on the mean of the other covariates. We overlaid the response curves of all species that shared the covariate effect, and described the general patterns in availability probabilities across species using a 2-dimensional kernel density estimate with a bivariate normal kernel (Wand and Jones 1995) using the R package KernSmooth (Wand 2015). We identified the contour for the 75% bivariate quantile to identify the highest density region, that is, where most of the species' response curves overlapped. We also identified 5%, 25%, 50%, 75%, and 95% contours based on the cumulative density conditional on the covariate values, that is, percentage of the species with the same or higher availability given the date and time. We show average responses for all 3 time-varying model classes (M₀^ϕ, M_f^ϕ, and M_f^c).

Removal Models vs. PIF Approach to Adjusting for Incomplete Availability

We compared adjustments for incomplete availability of individual species based on removal models with those used in the Partners in Flight (PIF) Population Estimates Database (Blancher et al. 2013) to determine how PIF estimates of continental population sizes might be improved by using removal model–based estimates of availability. PIF uses a time adjustment (here referred to as T_adj; Blancher et al. 2013) to account for diurnal variation in availability, which is calculated as the ratio of the maximum count estimated from a polynomial regression to the mean count based on 50 stops in the BBS data (Rosenberg and Blancher 2005, Blancher et al. 2013), and assumes that the probability of availability = 1 at the time of the maximum count.

We defined our species adjustment factor as the reciprocal of average availability (U_inv = 1 / p_avg). We calculated the correlation between T_adj and U_inv and compared the distributions across 129 species for which PIF adjustments were available (Blancher et al. 2013; Appendix Table 3). We also assessed how maximum availabilities (p_max) compared with the PIF assumption of the maximum being 1. We calculated average (p_avg) and maximum (p_max) availability probabilities by using the best-supported model for each species (Table 1, Appendix Table 3) for predicting the probability of availability for each point count in the entire BAM dataset (Figure 1) given the covariate values, and then calculating the average and maximum of the predicted values for each individual species. We set count duration to be 3 min, to be consistent with the BBS methodology that T_adj was based on.

Sample Size Requirements for Time-invariant Removal Models

We were able to estimate the singing rate parameter (ϕ) for all 152 species using the conventional removal model (M₀) fit to the full dataset from the boreal study area (Appendix Table 2). The mean probability of availability increased from 0.51 ± 0.15 SD for 3-min point count duration to 0.68 ± 0.15 SD for 5-min and 0.88 ± 0.12 SD for 10-min count duration (Figure 2). A logistic regression model (β₀ = −1.16 ± 0.10 SE, β_{log(# nonzero counts)} = 0.92 ± 0.04 SE) predicted that a species needed to be detected during 45 ± 6 SE point counts for the M₀ model to have a 90% probability of successfully converging on a project-specific estimate of ϕ. Sample size requirements for Var[log(ϕ)] were identical to that of the point estimates.

FIGURE 2

Estimates of availability probabilities for landbird species in northern North America from conventional (M₀, 152 species) and finite mixture (M_f, 147 species) removal models relative to point count duration (min). Dotted vertical lines indicate common point count durations (3, 5, and 10 min).

We jointly estimated the parameters for the M_f models (ϕ, c) for 147 species using the full dataset. Mean probabilities of availability based on the M_f models across the 147 species were 9–16% lower than for the M₀ models, and increased from 0.46 ± 0.17 SD for 3-min, to 0.56 ± 0.17 SD for 5-min, to 0.74 ± 0.17 SD for 10-min count duration (Figure 2). A logistic regression model (β₀ = −1.16 ± 0.09 SE, β_{log(# nonzero counts)} = 0.66 ± 0.03 SE) predicted that a species needed to be detected during 158 ± 24 SE point counts for the M_f model to have a 90% probability of successfully converging on project-level estimates of ϕ and c. The sample size requirement for a 90% probability of successfully estimating the associated variance was substantially larger, 963 ± 177 SE (β₀ = −1.93 ± 0.09 SE, β_{log(# nonzero counts)} = 0.60 ± 0.03). The number of nonzero counts for the 5 species for which we were unable to estimate M_f model parameters varied between 88 and 419 (Appendix Table 2).

Adjusting Point Counts for Variation in Count Duration

The bias when correcting 3- and 5-min duration counts to 10-min counts was affected by count duration. A 5-min count duration decreased median bias and tightened the 90% quantile range of bias closer to 0 (no bias) across species for both the M₀ and M_f models, when compared with the distribution of 3-min based biases. The median bias with respect to 10-min counts for the M₀ model was 0.05 for 3-min and −0.03 for 5-min counts, and did not improve further with increasing sample size (Figure 3). Although the 90% quantile range for M₀ models shortened with increasing sample size, it stayed centered around the biased median. Median bias with respect to 10-min counts for the M_f model was small (absolute value <0.01) regardless of sample size; the 90% quantile range shortened with sample size. The high end of the 90% quantile range was close to 0 (no bias), whereas the lower limit spread wider into the negative direction, changing from −0.27 at a sample size of 20 to −0.12 at a sample size of 200 (Figure 3).

FIGURE 3

Sample size effects on (A) bias and (B) variance in removal model predictions of the number of individuals of a species counted during a 10-min point count when using the counts of individuals from 3-min (gray) and 5-min (black) point counts. Data from 133 landbird species with varying numbers of nonzero counts (horizontal axis) were used. The distributions represented by each box are from species whose sample size was equal to or greater than the sample size noted on the horizontal axis. Model types (conventional [M₀] and finite mixture [M_f]) are shown side by side; boxes show the 90% quantile range, and lines (for bias) represent the median.

Variance was also greatly affected by count duration. Increasing the count duration from 3 to 5 min resulted in reducing the 90% limit for the M_f models from 0.06 to 0.03. Variance of the corrected relative counts under the M₀ model was <0.01 irrespective of duration and sample size, but nevertheless showing similar relative patterns. The 90% limit for variance of corrected relative counts under the M_f models decreased to 0.03 and 0.01 with 200 nonzero counts for the 3- and 5-min counts, respectively (Figure 3).

Probability of availability and corrected relative counts based on the time-invariant removal models (M₀ and M_f) varied among species with different sample sizes (Figure 4). The Rusty Blackbird (Euphagus carolinus), with 254 nonzero counts, showed high uncertainty in both probabilities and corrected relative counts, but bias with respect to 10-min counts was small for both models because the finite mixture model intercept (probability at time 0, 1 − c) was close to 0 (c = 0.96 ± 0.18 SE), and thus the singing rates were similar (M₀: ϕ = 0.16 ± 0.00 SE; M_f: ϕ = 0.14 ± 0.05 SE). The high uncertainty for the M_f model was a result of imprecise c estimates and high (0.6) correlation between c and ϕ (high rate with low intercept (1 − c), and vice versa). The Western Wood-Pewee (Contopus sordidulus) had 493 nonzero counts in the validation set, and, as a result, its M_f estimates, especially c, were more precise, with low correlation (ϕ = 0.11 ± 0.05 SE, c = 0.73 ± 0.03 SE, correlation = 0.15). As the intercept (1 − c) was higher, the singing rates differed more between models (M₀: ϕ = 0.27 ± 0.00 SE). This resulted in larger differences between the 2 models in corrected relative counts, with M_f being unbiased, although uncertainties (ranges of confidence intervals) were comparable. The Connecticut Warbler (Oporornis agilis) had 1,162 nonzero counts, and had the most precise parameter estimates of the 3 example species and the highest intercept (M₀: ϕ = 0.40 ± 0.00 SE; M_f: ϕ = 0.17 ± 0.04 SE, c = 0.54 ± 0.03 SE, correlation = 0.56), resulting in low bias and variance compared with the other 2 species.

FIGURE 4

Availability probabilities (left) estimated from time-invariant conventional (M₀) and finite mixture (M_f) removal models for 3 example species, the Rusty Blackbird (Euphagus carolinus; RUBL), Western Wood-Pewee (Contopus sordidulus; WEWP), and Connecticut Warbler (Oporornis agilis; CONW), with 95% confidence intervals (gray shading). Corrected relative counts (right) were calculated as the corrected count divided by the 10-min-based corrected count. Ratio = 1 indicates unbiasedness; continuous vertical lines are 95% confidence intervals. Dotted vertical lines indicate common point count durations (3 and 5 min).

Relative Support for Models and Plotting Time-varying Effects

We fit 44 conventional and finite mixture removal models to 147 species detected during at least 75 point counts in the full dataset (fitting finite mixture models was unsuccessful for 5 species; Appendix Tables 2 and 3). Finite mixture models were best-supported for 92% of the 147 species (M_f^ϕ best-supported for 26%, M_f^c best-supported for 66%; Table 1). The only species for which the time-invariant (M₀) model was best-supported based on AIC_c was the Yellow-billed Cuckoo (Coccyzus americanus; n = 196). The best-supported models included time since local sunrise (SR) for 75%, ordinal day (DY) for 39%, and days since local spring (LS) for 37% of the species. Best models included both time since sunrise and either ordinal day or days since spring for 52% of the species (Table 1).

Availability probabilities generally declined across the observed ranges of times since sunrise, and were consistently high between 0 and 4 hr after sunrise. Responses were more variable before and after this 4-hr period. Species' responses to time since sunrise were predominantly nonlinear (77 species with a quadratic and 33 with a linear relationship). Availability declined slowly with respect to date, and was highest across species between ordinal days 155 and 175 and 30–60 days after local spring. Species' responses to dates were predominantly linear for both ordinal day (40 species with a linear relationship and 18 with a quadratic relationship) and days since spring (35 species with a linear and 19 with a quadratic relationship). The temporal patterns of predicted availabilities across all species were very similar between the different classes of time-varying removal models (M₀^ϕ, M_f^ϕ, and M_f^c; Figure 5).

FIGURE 5

Availability probabilities estimated from 3 classes of time-varying removal models—(A) M₀^ϕ (conventional with time-varying rate parameter), (B) M_f^ϕ (finite mixture with time-varying rate parameter), and (C) M_f^c (finite mixture with time-varying proportion of infrequent singers)—with covariate effects of hour since local sunrise, ordinal day, and day relative to local spring overlaid for all terrestrial landbird species in northern North America for which the respective variable was part of the best-fit model (see Table 1; other variables were conditioned at their means). Availability probabilities were based on 3-min point count duration. Darker colors in the 2-dimensional kernel density plots indicate regions in which time-varying response curves of multiple species overlap; the continuous contour line indicates the highest density region (75% quantile of the bivariate density). Dashed contours indicate cumulative proportions of species with availability probabilities less than or equal to the value in the contour label conditional on the value of the focal covariate on the horizontal axis.

The selection frequency for ordinal day (DY) was higher for migratory species (41%, 51 out of 124) than for winter residents (30%, 7 out of 23), but the difference was not significant (χ² = 0.54, P = 0.47). The selection frequency for days since local spring (LS) was somewhat lower for migrants (35%, 44 out of 124) than for resident species (43%, 10 out of 23), also with a nonsignificant (χ² = 0.93, P = 0.34) difference (Table 1).

Removal Models vs. PIF Approach to Adjusting for Incomplete Availability

Across 129 species of landbird, the PIF adjustment factor (T_adj) for BBS counts had median value of 1.27 and ranged from 1.03 to 25.88. The mean of the average availability (p_avg) across these species was 0.43 ± 0.17 SD; the mean of the maximum availability (p_max) was 0.66 ± 0.23 SD. Our inverse availability based adjustment factor (U_inv) derived from fitted removal models was 1.7 times higher than T_adj (U_inv median = 2.57, range = 1.25–89.69); U_inv and T_adj were not correlated (Spearman's r = 0.005, P = 0.95). Average availability was very low, and, as a result, U_inv was very high for the Eastern Bluebird (n = 99, U_inv = 6,788,866.26) and Rusty Blackbird (n = 420, U_inv = 10,022.05) based on their best-supported finite mixture models (Appendix Tables 2 and 3). We show these results to demonstrate the problems that arise when numbers of detections are too small for the mixture model. These 2 species were excluded from the comparisons.

Conventional vs. Finite Mixture Models

We evaluated alternate formulations of the time-removal model (a conventional removal model and a finite mixture removal model), with and without covariates, for 152 bird species, using a compilation of point count data spanning northern North America. We found that predicted availability probabilities under conventional and finite mixture models were very similar in their ranges of probability values and shapes of their response curves to predictor variables. However, finite mixture models were supported for 92% of species with sufficient nonzero counts, in accordance with previous studies that found that mixture models provided better fit (Farnsworth et al. 2002, Etterson et al. 2009, Reidy et al. 2011). By incorporating variation in singing behavior within a survey event, the finite mixture model relaxes the homogeneous singing rate assumption of the conventional removal models, thus reducing bias when using availability probabilities to adjust survey counts for detectability.

Count duration had a profound effect on average counts. In our data, the first 3-min and first 5-min counts contained, on average, 59% and 73% of the full 10-min counts from the same surveys. Both conventional and finite mixture models reduced this bias considerably. Finite mixture models adjusted the shorter-duration counts with low bias (<1%), whereas conventional models tended to underpredict counts by 5% for 3-min counts and overpredict counts by 3% for 5-min counts when compared with the 10-min counts from the same surveys. The corrected counts using the conventional models were characterized by much smaller variances (<1%) than those from finite mixture models (<6%). However, the variance of finite mixture model–based corrected counts reached acceptable levels for species with at least 200 nonzero counts, the sample size for which bias reduction was greatest.

We evaluated the number of nonzero counts needed to reliably estimate the parameters of the removal models, and found that conventional models required a minimum of 50 nonzero counts, whereas finite mixture models required at least 160 nonzero counts. This corroborates the substantial reduction in variance of the adjusted counts both within and among species at 200 nonzero counts. The low precision of the finite mixture model parameters can considerably inflate the variance of corrected relative counts at low sample sizes (Figures 3 and 4), and at least 1,000 nonzero counts were required to reliably estimate the variance–covariance matrix for finite mixture models. Because the number of nonzero counts was correlated with mean counts (relative abundances) of the study species, it follows that less abundant species will require larger samples for reasonable removal model estimates.

The 1,000 nonzero counts per species required for mixture models (Amundson et al. 2014) is somewhat greater than what is required by other models of detection probability, such as distance sampling (Buckland et al. 2001:240). However, when sample sizes were sufficient, the mixture models were nearly unbiased for most species when adjusting point counts for protocol variation. Because associated variances were often high, even with much larger sample sizes, we recommend that researchers compare results from conventional and finite mixture models, preferring finite mixture models when the estimates and their variances are not outlying values. Otherwise, the benefits of potential bias reduction might be obscured, especially at low detection rates.

Implications for Survey Design and Modeling

Counts of longer duration increased the numbers of detections and were therefore effective for reducing bias and variance in the estimates of availability and corrected relative counts. We found that increasing the count duration from 3 to 5 min reduced the range of variation in bias (with respect to 10-min counts) across species, and reduced variance by at least 50% for most species, using both conventional and finite mixture models. Thus, when designing field surveys, lengthening the count duration to 5 or 10 min is an important consideration to increase the accuracy and precision of population estimates by sampling design (Barker et al. 1993, Ralph et al. 1995, Matsuoka et al. 2014). We recognize that shorter count durations (3 min) are favored for logistical reasons, but the choice of count duration must be based on a careful consideration of the tradeoff between survey costs and study objectives.

Longer count durations increase the chances that individuals may move in or out of the surveyed area or will be double counted, which violates assumptions of the removal model and most other models of detectability applied to point counts. Double counting was found to comprise 0.9–3.4% of detections during experimental point counts (Simons et al. 2009), and may also have been low in our study. We excluded count intervals beyond 10 min to minimize the effect of movement-related biases on our estimates. Although longer-duration counts might have been inflated compared with counts in nested subintervals, this would not have changed our conclusions in terms of relative bias among different models.

We found that, among the 135 species for which finite mixture models were found to fit better than conventional models, for 72% there was more support for the models with a time-varying proportion of infrequent singers, when compared with the models with time-varying singing rates. This indicates that our proposed mechanism of individuals switching membership between the groups of frequent and infrequent singers is more realistic than the alternative of fixed proportions as used by previous studies (Farnsworth et al. 2002, Aldredge et al. 2007a, Etterson et al. 2009, Reidy et al. 2011). The overall patterns regarding responses of singing rates to dates and times were very similar, regardless of support across the model classes. The periods when we found detectability to be high (0–4 hr after sunrise and June 4–24 [ordinal day 155–175]; see Results) agreed well with long-standing protocol recommendations (Ralph et al. 1995). Conducting point counts during these times of day and year will make surveys much more efficient for accumulating detections and will thereby improve the resulting accuracy and precision of abundance estimates.

Besides these general guidelines, local and regional context and study goals should also shape survey timing. Although sunrise time already reflects geographic differences, weather and resulting local spring date can vary considerably across years, thus influencing the timing of migration for nonresidents, and the timing of high vocal activity for most species. Regional differences can also introduce heterogeneity into bird availability through geographic advancement of breeding through the growing season. We captured this breeding progression with the days since spring variable (vs. ordinal days). The fact that the 2 predictors related to time of year were supported by different species suggests differences among species' plasticity of response to interannual variations in weather conditions. We found no significant differences between winter resident and migratory species with respect to selection frequencies of time-of-year covariates. This indicates that, although ordinal day and multidecadal averages (our days since spring variable) are important predictors of availability, time-of-year covariates are inadequate for establishing linkages between breeding phenology and weather events for migratory species (Both and Visser 2001, Both et al. 2010). Estimating year-specific effects with respect to days since spring vs. ordinal day might be indicative of changes in breeding phenology with a changing climate, which will be important to account for in programs monitoring long-term changes in bird population sizes. Accounting for year-to-year variation in covariates is important given the large interannual variation in spring green-up (15–20 days in the boreal region; Delbart et al. 2006) and the strong indications of advancing phenology in boreal and hemiboreal regions with climate change (Linkosalo et al. 2009, Beaubien and Hamann 2011).

Conclusions

Our study demonstrates the importance of accounting for temporal heterogeneity in bird availability when analyzing avian point count data. This is especially true when combining data across studies with different locations, sampling timeframes, or protocols. Location influences the expected time of the breeding season, and the timing of surveys relative to the local breeding season affects availability, for example, through time-varying behaviors. It is difficult to retrospectively control for these sources of heterogeneity introduced by sampling design, and assuming constant availability in space and time is unrealistic, especially when analyzing heterogeneous and/or spatially extensive datasets. Thus, a model-based approach to analysis is required. We argue that the continuous time-removal modeling framework outlined in this manuscript is sufficient. Easily obtainable covariates, such as the ones used in this study, can help to better correct availability bias of counts.

Comparing our average availability-based corrections with the time adjustment factors developed by Partners in Flight (PIF; Rosenberg and Blancher 2005, Blancher et al. 2013), we found our removal model–based adjustment factors to be 1.6 times higher. This was largely due to the PIF assumption of 100% availability during the times of day with the maximum counts, whereas our maximum availability estimates for a 3-min point count averaged 66%. As a result, PIF population size estimates are most likely to be underestimates for most species. This availability-related underestimation is relatively small compared with the 5-fold difference in population sizes when adjusting point counts for perceptibility probability using distance sampling versus PIF's expert-derived adjustment for maximum detection distances (Matsuoka et al. 2012). Similarly, both field tests (Confer et al. 2008) and sensitivity analyses (Thogmartin 2010) have indicated that PIF population size estimates are more strongly influenced by the choice of adjustment for perceptibility (maximum detection distance) than availability (time adjustment). We feel that availability probabilities from removal models could replace PIF's time adjustment (Rosenberg and Blancher 2005, Blancher et al. 2013) to modestly improve population size estimates, particularly if covariates were included in the removal models to account for date- and time-related heterogeneity in singing rates (Thogmartin et al. 2006). We did not find a correlation between correction factors based on removal models and PIF's time adjustments, even after adjusting for maximum availability. This is most likely due to geographic differences (mostly nonboreal vs. boreal regions, respectively) and the on- vs. off-road nature of the BBS and BAM datasets, which are complexities to be tackled by future studies.

The probability of availability, given the presence of an individual, can be efficiently estimated by conventional or finite mixture type removal models, with or without covariates for bird behavior. In this paper, we have outlined a flexible and continuous-time parameterization of these model types implemented in R package detect (Sólymos et al. 2018). We have extended removal model methodology to flexibly model time-varying finite mixture models using covariates for either the singing rate (ϕ) or the proportion of infrequent singers (c) parameter. We found that the choice of model had an impact on the estimability of unknown model parameters, and that this choice also affected the bias and variance of the corrected relative counts. In general, a conventional removal model is recommended for smaller sample sizes. Model selection approaches can be used to compare support for different parameterizations, but only at much higher sample sizes (>200 nonzero counts, but preferably >1,000). Lengthening the count duration from 3 min to 5 or 10 min increased sample sizes of detections and improved the accuracy and precision of removal model estimates. Well-informed survey design combined with various forms of removal sampling are useful for accounting for availability bias in point counts, thereby improving population size estimates and allowing for better integration of disparate studies at large spatial scales.