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
Abundance model
We modeled abundance and occupancy using separate models. We describe the abundance model here and subsequently describe the occupancy model.
State process
The abundance model includes a sub-model quantifying occupancy of grid cells, upon which abundance is conditioned (i.e., we estimate abundance within occupied grid cells).
\[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 abundance of independent clusters of individuals within \(d_{[max,high]}\) of points within an occupied grid cell as
\[N_{ijt} \sim Poisson(λ_{ijt}×z_{ijt}) \tag{D.3}\]
\[log(λ_{ijt})=α_{it}+ε_{jt}. \tag{D.4}\]
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.5}\]
\[ε_{it} \sim Normal(0,SD_{[ε]}^{2}). \tag{D.6}\]
We estimated \(SD_{[α],t}\) and \(SD_{[ε]}\) with gamma priors. For \(SD_{[α],t}\), we used a completely uninformed \(gamma(1,1)\) prior. For \(SD_{[ε]}\), we used a somewhat informed \(gamma(1.015117, 3.807330)\) prior based on the distribution of posterior estimates of this parameter across species from earlier analysis. We derive density, D, as
\[D_{ijt}= \frac{N_{ijt}× \hat{\bar{s}}}{K×∑_{d=1}^{d_{[max]}}A_{d}}, \tag{D.7}\]
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.8}\]
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}\)) as
\[n_{ijt} \sim Binomial\biggl(N_{ijt},p_{[a],tl}×p_{[p],t}×\biggl(\frac{K_{[s],ijt}}{K}\biggr)\biggr), \tag{D.9}\]
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, and \(K_{[s],ijt}⁄K\) is survey effort. Detection components can be expressed as the sum of probability vectors:
\[p_{[a],tl}=∑_{r=1}^{R} π_{[a],tlr} \tag{D.10}\]
\[p_{[p],tl}=∑_{d=1}^{d_{max}} π_{[p],tld}, \tag{D.11}\]
where \(π_{[a],tlr}\) and \(π_{[p],tld}\) are detection probabilities within time removal and distance bins (hereafter bin probabilities), respectively. We model the distribution of counted individuals among time removal and distance bins, respectively, as
\[\mathbf{y_{[a],ijt}} \sim Multinomial \biggl(n_{ijt},\mathbf{π_{[a],tl}^{[c]}} \biggr) \tag{D.12}\]
\[\mathbf{y_{[p],ijt}} \sim Multinomial \biggl(n_{ijt},\mathbf{π_{[p],tl}^{[c]}} \biggr), \tag{D.13}\]
where \(\mathbf{y_{[a],ijt}}\) and \(\mathbf{y_{[p],ijt}}\) are vectors of counts within time removal and distance bins. Bin probabilities for availability assume constant availability across the survey period,
\[π_{[a],trl}=p_{[ar],tl} (1-p_{[ar],tl})^{r-1}, \tag{D.14}\]
where \(p_{[ar],tl}\) is the availability for a single 2-min time removal bin. Bin 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.15}\]
where \(p_{[pd],td}\) is the perceptibility within distance bin \(d\) in year \(t\), and \(A_d\) is the area of distance bin \(d\). Multinomial cell probabilites for modeling \(\mathbf{y_{[a],ijt}}\) and \(\mathbf{y_{[p],ijt}}\) represent bin probabilites scaled to sum to 1:
\[π_{[a],trl}^{[c]}=π_{[a],trl}⁄p_{[a],tl}, \tag{D.16}\]
\[π_{[p],td}^{[c]}=π_{[p],td}⁄p_{[p],t}. \tag{D.17}\]
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}\). 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.18}\]
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.19}\]
where \(h_t\) and \(b\) are estimated parameters for the hazard function. For perceptibility models (Equation D.18 and Equation D.19), we draw year-specific shape parameters from hyperdistributions:
\[log(σ_t) \sim Normal(μ_{[σ]},SD_{[σ]}^2), \tag{D.20}\]
\[log(h_t) \sim Normal(μ_{[h]},SD_{[h]}^2). \tag{D.21}\]
We allowed availability to vary by year and BCR:
\[logit(p_{[ar],tl}) \sim Normal(μ_{[p_{ar}]},SD_{[p_{ar}]}^2). \tag{D.22}\]
Priors for annual variation in detectability components are
\[μ_{[σ]} \sim Normal(4.5,1.5), \tag{D.23}\]
\[μ_{[h]} \sim Normal(4.5,1.5), \tag{D.24}\]
\[μ_{[p_{ar}]} \sim Normal(0,1), \tag{D.25}\]
\[SD_{[σ]},SD_{[h]},SD_{[p_{ar}]} \sim gamma(1,1). \tag{D.26}\]
Occupancy model
We describe here the occupancy model, which estimates occupancy of points conditioned on occupancy of grid cells.
State process
We model species occupancy of grid cells in the same manner as for the abundance model (Equation D.1, Equation D.2). We model species occupancy of points within an occupied grid cell as
\[u_{ijt}|z_{ijt} \sim Binomial(K,θ_{it}×z_{ijt}) \tag{D.27}\]
\[logit(θ_{it}) \sim Normal(μ_{[θ],t},SD_{[θ],t}^{2}), \tag{D.28}\]
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.
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.29}\]
with \(p_{tl}\) is a year- and BCR-specific detection probability and \(K_{[s],ijt}⁄K\) is survey effort. Overall species detectability for a point survey, \(p_{tl}\), summarizes the probability vector
\[p_{tl}=∑_{r=1}^{R} \mathbf{π_{[R],tl}^{[c]}} \tag{D.30}\]
consisting of multinomial cell probabilities for
\[\mathbf{y_{[R],ijt}} \sim Multinomial \biggl(y_{ijtl}, \mathbf{π_{[R],tl}^{[c]}} \biggr), \tag{D.31}\]
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.32}\]
and are scaled to sum to 1,
\[π_{[R],trl}^{[c]}=π_{[R],trl}⁄p_{tl}, \tag{D.33}\]
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 abundance model (\(p_{[ar]}\) in Equation D.22, Equation D.25, and Equation D.26), including the shorter 5-min survey period in 2008-2009 (see explanation following Equation D.16).
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.34}\]
\[logit(\hat{ψ}_{itm}×\hat{θ}_{it}) \sim α_{[ψθ],im}+r_{[ψθ],im}×(t-1). \tag{D.35}\]
For occupancy, we multiplied \(ψ×θ\) to derive the unconditional probability of a point being occupied by the species, and then derived trends. We retained the maximum likelihood estimate (MLE) of the slope, \(\tilde{r}_{im}\), which is equivalent to the least squares trend. We also retained the MLE of the intercept \(\tilde{α}_{im}\). 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.11). The plot size for abundance estimation technically corresponds with truncation distance, but because density is abundance scaled by area, interpretation of density does not depend on plot size.
For occupancy, we truncated detections for most species at 125 m from survey points (i.e., half the distance between neighboring points). Ninety-five of the 358 species analyzed are typically detected at distances beyond 125 m such that truncating at this distance would entail excluding most detections. For these species, we truncated at 998 m. 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 125m radius plots, neighboring plots within a grid cell do not overlap, so estimation of area occupied is readily interpretable as the species range within the stratum. For species with 998m radius plots, neighboring plots can overlap, so a single individual can potentially occupy > 1 plot. For these species, occupancy less readily interpretable as the species range and instead more clearly quantifies space use (Latif et al. 2016). Occupancy estimates are not entirely comparable across species with different plot sizes. Nevertheless, large differences in occupancy should be ecologically meaningful. 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 alternative structures for the observation component of abundance and occupancy models in order to balance model fit, complexity, and computational efficiency. 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 for abundance models.
For abundance models, 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. We first considered the maximally complex abundance model (hereafter global abundance model), which included annual variation in both perceptibility and availability, and variability in availability among Bird Conservation Regions (Equation D.20 - Equation D.22). If the global abundance model did not converge after a specified number of attempts (for details on model fitting, convergence, and MCMC sampling criteria, see Section D.1.7 below), we then considered a constant perceptibility abundance model, which replaced either Equation D.20 and Equation D.23 or Equation D.21 and Equation D.24 with
\[log(σ) \sim Normal(4.5,1.5) \tag{D.36}\]
\[log(h) \sim Normal(4.5,1.5), \tag{D.37}\]
respectively. We successfully fitted models allowing variation in availability for all species, so we had no need of a constant availability model, even though we did initially consider it.
For occupancy models, we occasionally failed to achieve convergence when allowing detectability to vary among years and Bird Conservation Regions, in which case we fitted a constant-detectability model. Since detectability for occupancy models is analogous to the availability component in abundance models, the constant-detectability occupancy model replaces the equivalent of Equation D.22, Equation D.25, and Equation D.26 with
\[logit(p) \sim Normal(0,1). \tag{D.38}\]
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 20\), 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 50 max attempts for global abundance and occupancy models and no max for simpler alternatives.