Appendix D — Data Analysis
D.1 Statistical Models
Indices
The following indices are used in all model descriptions:
\(i\) - strata with max \(I\)
\(j\) - grid cells with max \(J\)
\(k\) - points within a grid cell with max \(K=16\), and \(K_{[s]}\) is the number of points surveyed in grid cell \(j\)
\(l\) - Bird Conservation Regions with max \(L\); note that strata are nested with Bird Conservation Regions, so an index \(i\) implies indexing by \(l\).
\(d\) - distance bins with max \(d_{[max]}=10\)
\(d_{[mid]}\), \(d_{[low]}\), and \(d_{[high]}\) are the distance mid-point, lower bound, and upper bound for distance bin \(d\). \(d_{[max,high]}\) is the upper bound for the furthest distance bin (i.e., the radius of the entire point-level plot).
\(r\) – 2-min time removal bins within the 6-min survey with max \(R=3\)
\(t\) – years with max \(T\)
\(m\) – posterior samples
Hurdle model (unified abundance & occupancy)
For most species, we modeled abundance and occupancy using a unified hurdle model. The “hurdle” refers to the point-level portion of the state process whereby we assume at least one individual present at each occupied point, and we model the number of additional individuals present over and above the first. Thus, with a hurdle structure zeros can only arise from the occupancy portion of the model, in contrast with zero inflation wherein zeros can arise from either the occupancy or abundance sub-models. By removing the ambiguity of how zeros arise, the hurdle model allows us to more clearly quantify point-level occupancy distinctly from abundance using the same model. Because abundance and occupancy models largely represent subsets of the sub-models contained in the hurdle model, we first describe the hurdle model in complete detail, and then more briefly describe the abundance and occupancy models, focusing on how the latter differ from the former.
State process
We model species occupancy of grid cells as
\[z_{ijt} \sim Bernoulli(ψ_{it}) \tag{D.1}\]
\[logit(ψ_{it}) \sim Normal(μ_{[ψ],t},SD_{[ψ],t}^{2}), \tag{D.2}\]
where \(ψ_{it}\) is the probability of species occupancy of a grid cell in stratum \(i\) in year \(t\), \(z_{ijt}\) is the latent occupancy state of grid cell \(j\) in stratum \(i\) during year \(t\), and \(μ_{[ψ],t}\) and \(SD_{[ψ],t}^{2}\) are year-specific mean and variance hyper-parameters. We exclude data for all strata in which the species was never detected (effectively conditioning occupancy on being within the species range). We model species occupancy of points within an occupied grid cell as
\[u_{ijt}|z_{ijt} \sim Binomial(K,θ_{it}×z_{ijt}) \tag{D.3}\]
\[logit(θ_{it}) \sim Normal(μ_{[θ],t},SD_{[θ],t}^{2}), \tag{D.4}\]
where \(θ_{it}\) is the probability of species occupancy of a point within an occupied grid cell in stratum \(i\) and year \(t\), \(u_{ijt}\) is the latent occupancy state describing the number of points occupied within grid cell \(j\), stratum \(i\), and year \(t\), and \(μ_{[θ],t}\) and \(SD_{[θ],t}^{2}\) are year-specific mean and variance hyper-parameters. We model abundance of independent clusters of individuals within \(d_{[max,high]}\) of points within a grid cell as
\[N_{ijt} \sim Poisson(λ_{ijt}×u_{ijt})+u_{ijt} \tag{D.5}\]
\[log(λ_{ijt})=α_{it}+ε_{jt}. \tag{D.6}\]
where \(α_{it}\) and \(ε_{jt}\) are stratum- and grid-cell-specific random effects, respectively, and \(α_{it}\) is also subject to a fixed effect of year:
\[α_{it} \sim Normal(μ_{[α],t},SD_{[α],t}^{2}) \tag{D.7}\]
\[ε_{it} \sim Normal(0,SD_{[ε]}^{2}). \tag{D.8}\]
We derive density, D, as
\[D_{ijt}= \frac{N_{ijt}× \hat{\bar{s}}}{K×∑_{d=1}^{d_{[max]}}A_{d}}, \tag{D.9}\]
where \(K×∑_{d=1}^{d_{[max]}}A_{d}\) is the effective area (\(km^{2}\)) of estimation for a grid cell and \(\hat{\bar{s}}\) is the estimated mean size for independent clusters, which is sampled from the distribution
\[\hat{\bar{s}} \sim Normal(\hat{\mu_{s}},\hat{SE_s}), \tag{D.10}\]
where \(\hat{\mu_s}=Mean(s)\) and \(\hat{SE_s}=SD(s)/\sqrt{n}\), in which \(s\) represents observed cluster sizes and \(n\) is the number of independent clusters detected. We derive stratum-level density estimates by averaging all \(D_{ijt}\) within each stratum, and we take the area-weighted average of strata estimates to obtain superstrata density and occupancy estimates.
Observation process
We model the sum of observed counts (of independent clusters of individuals) across all points surveyed within a grid cell (\(n_{ijt}\)) and the number of points where the species was detected (\(y_{ijt}\)) as
\[n_{ijt} \sim Binomial\biggl(N_{ijt},p_{[a],tl}×p_{[p],t}×\biggl(\frac{K_{[s],ijt}}{K}\biggr)\biggr), \tag{D.11}\]
\[y_{ijt} \sim Binomial\biggl(u_{ijt},p_{[s],ijt}×\biggl(\frac{K_{[s],ijt}}{K}\biggr)\biggr),~\mathrm{and} \tag{D.12}\]
\[p_{[s],ijt}=1-(1-p_{[a],tl}×p_{[p],t})^{N_{ijt}}, \tag{D.13}\]
where \(p_{[a],tl}\) is a year- and BCR-specific availability component of detection probability, \(p_{[p],t}\) is a year-specific perceptibility component of detection probability, \(K_{[s],ijt}⁄K\) is survey effort, and \(p_{[s],tl}\) is the probability of detecting the species at an occupied point. Detection components summarize probability vectors,
\[p_{[a],tl}=∑_{r=1}^{R} \mathbf{π_{[a],tl}} \tag{D.14}\]
\[p_{[p],tl}=∑_{d=1}^{d_{max}} \mathbf{π_{[p],tl}}, \tag{D.15}\]
consisting of multinomial cell probabilities for
\[\mathbf{y_{[a],ijt}} \sim Multinomial \biggl(n_{ijt},\mathbf{π_{[a],tl}^{[c]}} \biggr) \tag{D.16}\]
\[\mathbf{y_{[p],ijt}} \sim Multinomial \biggl(n_{ijt},\mathbf{π_{[p],tl}^{[c]}} \biggr), \tag{D.17}\]
where \(\mathbf{y_{[a],ijt}}\) and \(\mathbf{y_{[p],ijt}}\) are vectors of counts within time removal and distance bins, respectively. Cell probabilities for availability assume constant availability across the survey period,
\[π_{[a],trl}=p_{[ar],tl} (1-p_{[ar],tl})^{r-1}, \tag{D.18}\]
and are scaled to sum to 1,
\[π_{[a],trl}^{[c]}=π_{[a],trl}⁄p_{[a],tl}, \tag{D.19}\]
where \(p_{[ar],tl}\) is the availability for a single 2-min time removal bin. Surveys in 2008-2009 were 5 min instead of 6 min in length, so for these years, we calculate the availability for a 1-min interval, \(p_{[a1],tl}=1-\biggl((1-p_{[ar],tl})^3\biggr)^{1/6}\), and then calculate the unscaled cell probability for the third (1 min) period as \(π_{[a],(t≤2)(r=3)l}=p_{[a1],tl}×(1-p_{[ar],tl})^{r-1}\). Cell probabilities for perceptibility account for area of bins with increasing distance:
\[π_{[p],td}=\frac{p_{[pd],td}A_d}{∑_{d=1}^{d_{[max]}}A_d},~\mathrm{and} \tag{D.20}\]
\[π_{[p],td}^{[c]}=π_{[p],td}⁄p_{[p],t}, \tag{D.21}\]
where \(p_{[pd],td}\) is the perceptibility within distance bin \(d\) in year \(t\), and \(A_d\) is the area of distance bin \(d\). We model \(p_{[pd],td}\) using one of four possible perceptibility models: 1) a half-normal constant model, 2) a hazard rate constant model, 3) a half-normal year model, and 4) a hazard rate year model. For half-normal functions, we calculated $p_{[pd],td} for each distance class as:
\[p_{[pd],td}=\frac{2\pi\int_{c=d_{[low]}}^{c=d_{[high]}}exp\biggl(-\biggl(\frac{c^2}{2\sigma_{t^2}}\biggr)\biggr)c~dc}{A_d}, \tag{D.22}\]
where \(d_{[low]}\) and \(d_{[high]}\) are the bounds for distance bin \(d\) and \(σ_t\) is the half-normal shape parameter. Because of the lack of an analytical solution to the integral of the hazard rate function, we calculated \(p_{[pd],td}\) at the midpoint, \(d_{[mid]}\), of each distance class:
\[p_{[pd],td}=1-exp\biggl(-\biggl(\frac{d_{[mid]}}{h_t}\biggr)^b\biggr), \tag{D.23}\]
where \(h_t\) and \(b\) are estimated parameters for the hazard function. For perceptibility models (Equation D.22 and Equation D.23), we draw year-specific shape parameters from hyperdistributions:
\[log(σ_t) \sim Normal(μ_{[σ]},SD_{[σ]}^2), \tag{D.24}\]
\[log(h_t) \sim Normal(μ_{[h]},SD_{[h]}^2). \tag{D.25}\]
We also consider models allowing availability to vary by year and BCR:
\[logit(p_{[ar],tl}) \sim Normal(μ_{[p_{ar}]},SD_{[p_{ar}]}^2). \tag{D.26}\]
Priors for annual variation in detectability components are
\[μ_{[σ]} \sim Normal(4.5,1.5), \tag{D.27}\]
\[μ_{[h]} \sim Normal(4.5,1.5), \tag{D.28}\]
\[μ_{[p_{ar}]} \sim Normal(0,1), \tag{D.29}\]
\[SD_{[σ]},SD_{[h]},SD_{[p_{ar}]} \sim gamma(1,1). \tag{D.30}\]
Abundance-only model
For species that we modeled abundance separately from occupancy, we describe here the abundance model, focusing on the model components that differed from the hurdle model.
State process
We model species occupancy of grid cells in the same manner as for the hurdle model (Equation D.1 and Equation D.2). Because we model occupancy separately (see Section D.1.4 below), we omit Equation D.3 and Equation D.4 from the abundance model. Without point occupancy, we model abundance at occupied grid cells (in contrast with Equation D.5) more simply as
\[N_{ijt} \sim Poisson(λ_{ijt}×z_{ijt}), \tag{D.31}\]
with variability among strata and grid cells modeled the same (Equation D.6, Equation D.7, and Equation D.8). Derivation of density, D, also does not change (Equation D.9, Equation D.10).
Observation process
We model the sum of observed counts in the same manner as the hurdle model (Equation D.11), but we omit observation models for point occupancy (i.e., omit Equation D.12 and Equation D.13). The remaining observation sub-models exactly match those of the hurdle model (Equation D.14 - Equation D.30).
Occupancy-only model
For species that we modeled occupancy separately from abundance, we describe here the occupancy model, focusing on the model components that differed from the hurdle model.
State process
We model species occupancy of grid cells and points in the same manner as for the hurdle model (Equation D.1 - Equation D.4). The remainder of the ecological portion of the hurdle model (Equation D.5 - Equation D.10) are ommitted from the occupancy model.
Observation process
We model species detections at points using
\[y_{ijt} \sim Binomial\biggl(u_{ijt},p_{tl}×\biggl(\frac{K_{[s],ijt}}{K}\biggr)\biggr), \tag{D.32}\]
with \(p_{tl}\) here being equivalent to \(p_{[s],ijt}\) in Equation D.12. Overall species detectability for a point survey, \(p_{tl}\), summarizes the probability vector
\[p_{tl}=∑_{r=1}^{R} \mathbf{π_{[R],tl}^{[c]}} \tag{D.33}\]
consisting of multinomial cell probabilities for
\[\mathbf{y_{[R],ijt}} \sim Multinomial \biggl(y_{ijtl}, \mathbf{π_{[R],tl}^{[c]}} \biggr), \tag{D.34}\]
where \(\mathbf{y_{[R],ijt}}\) is a vector of representing the frequency of points where the species was detected across time removal bins when first detected. Cell probabilities for the time removal model assume constant detectability across the survey period,
\[π_{[R],trl}=p_{[R],tl} (1-p_{[R],tl})^{r-1}, \tag{D.35}\]
and are scaled to sum to 1,
\[π_{[R],trl}^{[c]}=π_{[R],trl}⁄p_{tl}, \tag{D.36}\]
where \(p_{[R],tl}\) is species detectability for a single 2-min time removal bin. We model variability in species detectability here (\(p_{[R],tl}\)) in the same manner as the availability component of the hurdle model (\(p_{[ar]}\) in Equation D.26, Equation D.29, and Equation D.30), including the shorter 5-min survey period in 2008-2009 (see explanation following Equation D.19).
Trend estimation
We derived posterior estimates of occupancy and abundance trends for any strata or superstrata for desired time periods from strata or (rolled-up) superstrata estimates of occupancy or density, respectively. For stratum or superstratum, i, we fitted a generalized linear model (GLM) to each posterior sample, m, of density or occupancy estimates:
\[log(\hat{D}_{itm}) \sim α_{[D],im}+r_{[D],im}×(t-1)~\mathrm{or} \tag{D.37}\]
\[logit(\hat{ψ}_{itm}×\hat{θ}_{it}) \sim α_{[ψθ],im}+r_{[ψθ],im}×(t-1). \tag{D.38}\]
For occupancy, we multiplied \(ψ×θ\) to derive the unconditional probability of a point being occupied by the species, and then derived trends. To correctly propagate GLM estimation uncertainty, we sampled a random intercept \(\tilde{α}_{im}\) and slope \(\tilde{r}_{im}\) from the mean and standard error estimated for these parameters for each posterior sample. We then summarized across \(m\) to derive the posterior estimates for \(r\) and predicted \(D\) and \(ψ×θ\).
Plot size
Following established practice for distance sampling data (Buckland et al. 2001:15, 151), we truncated detections for each species beyond the 90% quantile of recorded distances from the surveyor. The purpose of this practice is to improve the fit of the perceptibility curve by minimizing the influence of outliers on estimation (e.g., \(p_{[p],tl}\) in Equation D.15). In our unified framework, truncation sets the plot size for estimation of density and occupancy. Because density is abundance scaled by area, interpretation of density does not depend on plot size, but interpretation of occupancy requires consideration of plot size depending on desired interpretation.
Occupancy describes the proportion of plots where at least one member of the species occurs (i.e., presence-absence of the species). As such, occupancy can be multipled by stratum area to estimate the area occupied, and plot size is the resolution for this quantity. For species with a plot size < 125m radius (i.e., the maximum plot size whereby neighboring plots within a grid cell do not overlap), estimation of area occupied is readily interpretable as the species range within the stratum. For plot sizes > 125m radius, neighboring plots can overlap, so a single individual can potentially occupy > 1 plot. Thus, for species with plot radii > 125m, area occupied is less readily interpretable as the species range and instead more clearly quantifies space use (Latif et al. 2016). Because of the difference in plot size among species, occupancy estimates are not entirely comparable across species. Nevertheless, large differences in occupancy should be ecologically meaningful. Moreover, application of estimates for quantifying the proportion or area of a stratum occupied would be comparable across species with plot sizes < 125m radius and roughly comparable across species with larger plot sizes. Trends in occupancy will also be comparable across species regardless of plot size. Additionally, trends in occupancy can be more clearly interpretable as changes in species range regardless of plot size, whereas trends in abundance may primarily reflect changes in population size.
Model selection
For each species, we considered a series of alternative model structures in order to balance model fit, complexity, and computational efficiency. Reflecting the estimation objectives for the analysis, the ecological component of all models shared the same structure, but varied in complexity of the observation process model. Because accounting for heterogeneity in bird observability is critical for accurate population estimation, we aimed to fit maximally complex observation models supported by the data. Moreover, because corrections to count data generated from distance sampling are relatively robust to heterogeneity (see pooling robustness in BUckland et al. 2001:41-42), we prioritized allowing heterogeneity in availability over perceptibility.
We first considered two potential functions for describing the decline in perceptibility of birds with distance from the surveyor: a half-normal model and a hazard model. We preliminarily fitted detection-only versions of these models to the raw distance data (i.e., distances not indexed by sampling unit or year). We used the Watanabe-Akaike Information Criterion (WAIC; Watanabe 2010, Hooten and Hobbs 2015) to select the most parsimonious perceptibility functional form and then retained that function for each species for modeling distance-based perceptibility in the primary analysis.
For primary analysis of data for each species, we first considered the maximally complex hurdle model (hereafter global hurdle model), which included annual variation in both perceptibility and availability, and variability in availability among Bird Conservation Regions (Equation D.24 - Equation D.26). If the global hurdle model did not converge after a specified number of attempts (for details on model fitting, convergence, and MCMC sampling criteria, see Section D.1.8 below), we then considered a constant perceptibility hurdle model, which replaced either Equation D.24 and Equation D.27 or Equation D.25 and Equation D.28 with
\[log(σ) \sim Normal(4.5,1.5) \tag{D.39}\]
\[log(h) \sim Normal(4.5,1.5), \tag{D.40}\]
respectively. If a constant perceptibility model also failed to converge, we switched to a global abundance model. We allowed for simpler models (i.e., a constant-perceptiblity abundance model and constant-availability hurdle and abundance models), but we only resorted to the global abundance model for 2% of species, and did not need anything simpler to estimate abundance for any species.
When estimating abundance with the abundance model, we estimated occupancy with occupancy models. As with abundance models, we were always able to meet convergence and sampling criteria with a global occupancy model (i.e., detectability varied among years and Bird Conservation Regions as in Equation D.26), so we never needed a constant detectability model.
Model fitting
We performed all data processing, analysis, and derivation of estimates in R (R Core Team, 2023) and model fitting in Nimble (de Valpine et al. 2017, 2023). All scripts for the annual analysis are maintained in a git repository and available upon request.
The R code automated combining strata-level population and trend estimates to derive superstrata estimates at multiple spatial scales, including quantification of estimation uncertainty. Additionally, we automated periodic checking of convergence and switching to alternative models in cases where more complex models failed to converge within a specified number of attempts.
Model fitting with Nimble entailed sampling model parameter space with 2 independent Markov Chain Monte Carlo (MCMC) chains. For each species, the automated process began sampling with chains of length = 10000 iterations, discarding the first 5000 as burn-in and retaining every 50th sample, then checked convergence of chains using the \(\hat{R}\) metric and effective sample size using \(n_{effective}\) (Gelman et al. 2004). We required \(\hat{R}<1.1\) and \(n_{effective}\geq 10\), and continued sampling for another 10000 iterations until we either achieved these criteria or reached the maximum number of attempts, with an attempt being 10000 iterations and checks of \(\hat{R}\) and \(n_{effective}\) occurring after each attempt. We specified a lower max attempts for the global hurdle and abundance models and a higher max attempts for the constant-perceptibility hurdle model. Initially, these levels were 10 and 20 attempts, respectively. For subsequent analyses implemented to troubleshoot errors for some species, we set max attempts higher to 20 and 50, respectively. Going forward, we may further tune the number of attempts, max \(\hat{R}\), and min \(n_{effective}\) to balance estimation rigor with computation speed and efficiency.