Trends in bird abundance differ among protected forests but not bird guilds
Jeffrey W. Doser1, 2, Aaron S. Weed3, Elise F. Zipkin2, 4, Kathryn M. Miller5,

Andrew O. Finley1, 2, 6

1Department of Forestry, Michigan State University, East Lansing, MI, 48824,
USA
2Ecology, Evolution, and Behavior Program, Michigan State University, East
Lansing, MI, 48824, USA
3Northeast Temperate Inventory and Monitoring Network, National Park Service,
Woodstock, VT 05091, USA
4Department of Integrative Biology, Michigan State University, East Lansing,
MI, 48824, USA
5Northeast Temperate Inventory and Monitoring Network, National Park Service,
Bar Harbor, ME, 04609, USA
6Department of Geography, Environment, and Spatial Sciences, Michigan State
University, East Lansing, MI, 48824, USA
Corresponding Author: Jeffrey W. Doser, telephone: (585) 683-4170; email:
<EMAIL_ADDRESS>ORCID ID: 0000-0002-8950-9895
Running Head: Trends in forest bird abundance

## Abstract

Improved monitoring and associated inferential tools to efficiently identify
declining bird populations, particularly of rare or sparsely distributed
species, is key to informed conservation and management across large spatio-
temporal regions. We assess abundance trends for 106 bird species in a network
of eight forested national parks located within the northeast U.S.A. from
2006-2019 using a novel hierarchical model. We develop a multi-species, multi-
region removal sampling model that shares information across species and parks
to enable inference on rare species and sparsely sampled parks and to evaluate
the effects of local forest structure. Trends in bird abundance over time
varied widely across parks, but species showed similar trends within parks.
Three parks (Acadia National Park and Marsh-Billings-Rockefeller and
Morristown National Historical Parks (NHP)) decreased in bird abundance across
all species, while three parks (Saratoga NHP and Roosevelt-Vanderbilt and
Weir-Farm National Historic Sites) increased in abundance. Bird abundance
peaked at medium levels of basal area and high levels of percent forest and
forest regeneration, with percent forest having the largest effect. Variation
in these effects across parks could be a result of differences in forest
structural stage and diversity. By sharing information across both communities
and parks, our novel hierarchical model enables uncertainty-quantified
estimates of abundance across multiple geographical (i.e., network, park) and
taxonomic (i.e., community, guild, species) levels over a large spatio-
temporal region. We found large variation in abundance trends across parks but
not across bird guilds, suggesting that local forest condition might have a
broad and consistent effect on the entire bird community within a given park.
Research should target the three parks with overall decreasing trends in bird
abundance to further identify what specific factors are driving observed
declines across the bird community. Understanding how bird communities respond
to local forest structure and other stressors (e.g., pest outbreaks, climate
change) is crucial for informed and lasting management.

Keywords: avian populations, Bayesian analysis, community dynamics,
hierarchical model, imperfect detection, National Park Service, removal
sampling, species richness

## Introduction

Developing sampling designs and efficient statistical methods to monitor
trends in species (e.g., Lany et al. 2020) and communities (e.g., Farr et al.
2019) is critical to inform management of landscapes, wildlife, and other
natural resources (Stem et al., 2005). Species occupancy and abundance trends
have been used extensively over the last decade to identify and prioritize
management tasks needed to protect wildlife (e.g., Guisan and Thuiller, 2005;
Zipkin et al., 2010; Vander Yacht et al., 2016). Species distribution models
and similar methods are key tools to evaluate changes in populations and
communities (e.g., Wilkinson et al. 2019), an increasingly important task as
climate change and habitat loss cause new challenges to conservation.
Monitoring of bird species distributions in particular has received much
attention in the ecological literature as a result of their popularity, wide
development of statistical methods to account for imperfect detection (e.g.,
Farnsworth et al. 2002; Royle 2004b), large-scale public science programs
(e.g., Breeding Bird Survey, eBird), and bird sensitivity to environmental
stressors (Bibby, 1999; Canterbury et al., 2000; Gregory et al., 2003).

While monitoring changes in species distributions over time is a common task,
assessments of abundance trends is less common and more difficult, but
critical as declines in abundance can lead to reductions in ecosystem services
(Inger et al., 2015; Rosenberg et al., 2019). Key studies have shown large
declines in bird abundance in the last several decades across extensive
spatio-temporal regions in North America (Stanton et al., 2018; Rosenberg et
al., 2019) and Europe (Inger et al., 2015). Such research, finding large
declines in both common and rare bird species, suggests the need for large-
scale monitoring of bird abundance.

Forest bird species in particular play an essential role in forest ecosystems
primarily via their foraging ecology that results in ecological services such
as scavenging carcasses, nutrient cycling, pest management, and seed dispersal
(Bruns, 1960; Whelan et al., 2008, 2010, 2015). Understanding changes in
forest bird abundance over time is a key task in maintaining ecological
integrity of forest ecosystems. Diversity of forest birds is strongly linked
to structural forest composition and complexity (MacArthur, 1958, 1964;
MacArthur and MacArthur, 1961; MacArthur et al., 1966; Karr and Roth, 1971;
Willson, 1974), but how that complexity and composition affects bird abundance
over time is unresolved.

Monitoring of bird abundances occurs primarily through the use of point-count
surveys, in which observers record all birds seen at a location for a set
amount of time. Sampling and accompanying modeling techniques to address
imperfect detection in count data have been widely developed over the last 50
years and include removal sampling (Farnsworth et al., 2002), distance
sampling (Buckland et al., 1993), repeated counts (Royle, 2004b), capture-
recapture (Otis et al., 1981; Royle, 2009), spatial capture-recapture (Efford,
2004; Borchers and Efford, 2008; Royle and Young, 2008), and useful
combinations of such methods (Amundson et al., 2014; Miller et al., 2019).
Many extensions of these modeling techniques exist to meet objectives of
different types of monitoring programs, such as to estimate population size of
multiple species at a given point in time (e.g., Chandler et al. 2013; Farr et
al. 2019) and assess trends in species abundance across larger time periods
(Dail and Madsen, 2011; Zhao and Royle, 2019). Recently, the multi-species
occupancy model (Dorazio and Royle, 2005; Gelfand et al., 2005; Dorazio et
al., 2006; Gelfand et al., 2006) was extended to a multi-region framework that
enables modeling of occupancy and species richness across independent spatial
units that comprise a broader network (Sutherland et al., 2016; Wright et al.,
2020). An analogous approach could model abundance across a large network to
provide inference on factors driving spatio-temporal changes in the abundances
of common and rare species.

In 2006, the Northeast Temperate Inventory and Monitoring Network (NETN) of
the National Park Service (NPS) began a volunteer monitoring program to
collect data on bird community abundance and composition across upland forest
habitats in eight parks. We used these data in a network-wide analysis to
characterize spatio-temporal dynamics of forest bird abundance from 2006-2019.
We estimated park-specific trends in bird abundance to evaluate whether
patterns in the bird community change among protected areas, which may signify
where and whether management efforts would be most beneficial. Within each
park, we also estimated trends of functionally similar groups of birds (i.e.,
guilds) and of individual species to determine if management should target the
whole community, specific guilds, or individual species. We also sought to
understand whether bird abundance is more driven by presence of breeding
habitat (i.e., percent of forest within a 1km radius) or quality of the
breeding habitat characterized by the amount of live tree basal area and
forest regeneration, which have been shown to influence bird species diversity
(Flaspohler et al., 2002; Rankin and Perlut, 2015; Zipkin et al., 2010; Vander
Yacht et al., 2016). To do this, we extend the basic removal sampling model
(Farnsworth et al., 2002) and the multi-species removal sampling model
(Chandler et al., 2013) to a multi-region, multi-species framework that
directly accounts for imperfect detection and shares information across
species within each park and across the network of parks. This novel
hierarchical approach enables estimation for rare species and estimation of
parameters at parks where there is a paucity of data.

We predicted high variability in abundance trends across parks and species,
which we hypothesized would reveal which species and locations should be
targeted for management interventions. We also predicted the amount of local
forest cover (percent forest within a 1km radius) would have the largest (and
positive) effect on forest bird abundance (Ladin et al., 2016). We expected
little variation in the direction and magnitude of trends of species within a
guild and larger variation among the different guilds (O’Connell et al.,
2000). For instance, we predicted forest regeneration would have a positive
effect on ground and shrub nesting species, as the amount of regeneration is
negatively correlated with the amount of deer browsing (Russell et al., 2001;
Augustine and DeCalesta, 2003), which has been shown to reduce the occupancy
probability of ground-nesting and understory species (Zipkin et al., 2010). We
expected forest interior obligate and canopy nesting species to show positive
relationships with the amount of basal area (Rankin and Perlut, 2015), while
shrub nesters and forest ground nesters will show either a negative
relationship with basal area or a curvilinear relationship, peaking at
intermediate values (Flaspohler et al., 2002; Rankin and Perlut, 2015; Vander
Yacht et al., 2016).

## Materials and Methods

### Study Site and Sampling Methods

The NETN monitored eight parks across the northeastern United States of
America (USA) from 2006-2019 (Figure 1, Faccio et al. 2015). Table 1 includes
the park acronyms used throughout the remainder of the manuscript, the total
amount of forested land in each park, and a summary of data available for
parks (as sampling intensity varied across years and parks). The number of
point count locations used to monitor forest birds (henceforth called points)
ranged across parks from 5 (WEFA) to 51 (ACAD). The NETN established point
count locations based on four criteria: 1) points were 200-250m apart to avoid
duplicate sampling; 2) points were located at least 50m from forest edges; 3)
points were located in the dominant forest type of each park; and 4) at least
one forest vegetation sampling plot was located within 50m of at least one
point in each dominant forest type. To meet these criteria, the NETN used a
systematic grid sampling frame to select points at all parks except ACAD,
where the NETN used the Generalized Random Tesselation Stratified algorithm
(Stevens, 1997; Stevens and Olsen, 2004). See Faccio et al. (2015) for further
details on point selection and protocol methods.

Volunteers performed 10-minute surveys at points located along permanent
forested, transects, annually during the breeding season (May-July). At each
point, an observer recorded the time of day and species of individual birds
they saw or heard within ten, one-minute increments. Observers recorded each
individual bird only the first time they detected the bird during the point
count. Observers performed point count surveys in the early morning and during
adequate weather conditions, when vocalization is most likely. There was large
variability in the number of points surveyed in any specific year within parks
as a result of volunteer availability (Table 1).

We used the time-interval of detections to estimate species’ abundances with a
model that incorporated detection probabilities and the removal sampling data
(Farnsworth et al., 2002; Nichols et al., 2009). We developed a model that
shares information across parks and species within a hierarchical framework
because data were sparse at some parks (i.e., WEFA, SAGA, SARA) and for rarely
detected species. This approach enables detection corrected estimates of
abundance at lesser sampled parks and for rare and/or elusive species, and
also yields more precise estimates for all species at all parks, analogous to
multi-species occupancy models (e.g., Zipkin et al. 2010).

### Forest Covariates

In addition to bird monitoring stations, the NETN has over 300 forest
vegetation monitoring plots implemented across the same eight parks (Tierney
et al., 2016). Some of these plots are co-located within 50m of the bird point
count survey points, which enables assessment of the relationship between
forest structure and bird abundance. For points that did not have an
associated forest plot within 50m, we used covariate data from the closest
forest vegetation monitoring plot occurring in the same general forest cover
type.

We evaluated the influence of three forest covariates: amount of forest
regeneration, amount of live basal area, and percent forest within a 1km
radius around the survey points. We used the amount of forest regeneration
(stems ha-1) as a metric to quantify advanced tree regeneration in the forest
understory, and computed regeneration as the total observed density of
seedlings (diameter at breast height (DBH) < 1 cm and at least 15 cm tall)
plus the total observed density of saplings (1 cm $\leq$ DBH $\leq$ 10 cm) per
plot. We calculated the amount of living basal area ($\text{m}^{2}/\text{ha}$)
for a given canopy tree in a plot from the diameter at breast height (DBH) and
then summed over all species in the plot. Field crews measured these forest
variables once every four years in each plot, resulting in 3-4 measurements of
basal area and regeneration for a given point over the study period. We
averaged these values so that each point was associated with a single value of
basal area and forest regeneration that represented the average of these
forest characteristics across the 14-year study period to provide a broad
overview of how forest condition influences forest bird abundance. We obtained
the percent forest within a 1km radius around the survey location in 2016
using the National Landcover Database (Homer et al., 2015). We standardized
all covariates to have mean zero and a standard deviation of one to facilitate
comparison of effect sizes and relative importance. Preliminary analyses
suggested linear effects of all three variables and a quadratic effect of
basal area on abundance would yield the best-fitting model.

### Modeling Framework

We developed a hierarchical model (e.g., Berliner 1996; Gelman et al. 2004;
Royle and Dorazio 2006; Royle and Dorazio 2008; Hobbs and Hooten 2015) to
estimate trends in community abundance and individual species abundance across
multiple geographically distinct parks (Sutherland et al. 2016; Wright et al.
2020; see Figure 2 for overview of modeling framework). To share information
across species and effectively model communities of species within each park,
we view species-level parameters as arising from common park-level
hyperparameters, which enables improved precision on species level estimates
(Yamaura et al., 2011, 2012). We extend this concept to another hierarchical
layer to share information across the entire meta-community by drawing park-
level hyperparameters from a common network-level distribution (Sutherland et
al., 2016; Wright et al., 2020). The model is thus a multi-level extension of
a generalized linear model (GLM) with the following hierarchical levels: 1)
detection model to accommodate imperfect detection for each species; 2)
process model of the true latent abundance of each species; 3) park level that
models species-specific parameters from a common park (i.e., community) level
distribution to induce dependence among species in a park and share
information across the community; and 4) network-level that views park-level
parameters as arising from a common network-level distribution to share
information across all parks in the network. While alternative approaches
exist to model species jointly via direct dependency of species detections
and/or abundances (e.g., Ovaskainen et al. 2010; Ovaskainen and Soininen 2011;
Warton et al. 2015), the hierarchical approach we employ is a mechanistic
alternative that enables inference at multiple ecological levels while
accounting for imperfect detection. This approach has been shown to provide
good model fit, ecological insights, and key information to support management
decisions (e.g., Zipkin et al. 2009; Chandler et al. 2013; Linden and Roloff
2013; Zipkin et al. 2020).

Following the removal sampling protocol, our model uses the time period of
first observation to account for imperfect detection of individuals during the
survey by estimating a detection probability that is the product of
availability and detectability (Farnsworth et al., 2002; Nichols et al., 2009;
Kéry and Royle, 2015). We used intervals of two-minutes in length to summarize
detections from the available data resulting in a total of $B=5$ intervals in
which an individual bird could be detected during the survey.

We implemented point count surveys using the removal sampling protocol at
$r=1,\dots,8$ geographically distinct parks for $i=1,\dots,n_{r}$ species at
$j=1\dots,J_{r}$ points for each of $t=1,\dots,T_{r}$ years. We define the
vector $\bm{y}_{r,i,j,t}$ as the number of individuals of species $i$
encountered during year $t$ at point $j$ within park $r$ in each of the $B$
time intervals. Subsequently, we define $\bm{y}^{*}_{r,i,j,t}$ as equivalent
to $\bm{y}_{r,i,j,t}$ with an additional value representing the number of
individuals that were not detected in any time interval. Following the basic
removal sampling model, we view the data $\bm{y}^{*}_{r,i,j,t}$ as arising
from a multinomial distribution with cell probabilities
$\bm{\pi}^{*}_{r,i,j,t}$ conditioned on the latent abundance $N_{r,i,j,t}$. We
then view the latent abundance $N_{r,i,j,t}$ as a Poisson distributed random
variable with mean $\lambda_{r,i,j,t}$. The multinomial observation model thus
allocates the total latent number of birds of species $i$ at point $j$ in park
$r$ during year $t$ ($N_{r,i,j,t}$) across the $B$ removal time periods and
the additional category for individuals not observed. Accordingly our models
for both the sampling and biological processes, respectively, take the
following forms

$\begin{split}\bm{y}^{*}_{r,i,j,t}&\sim\text{Multinomial}(N_{r,i,j,t},\bm{\pi}^{*}_{r,i,j,t})\\\
N_{r,i,j,t}&\sim\text{Poisson}(\lambda_{r,i,j,t}).\end{split}$

Given this formulation, we can use the relationship between the multinomial
and Poisson distributions to directly model $\bm{y}_{r,i,j,t}$ as a series of
conditionally independent Poisson distributions by analytically deriving the
marginal likelihood (Royle, 2004a; Dorazio et al., 2005; Kéry and Royle, 2015;
Yackulic et al., 2020). For each $b=1,\dots,B$ removal period, our marginal
likelihood is defined as

$y_{r,i,j,t,b}\sim\text{Poisson}(\lambda_{r,i,j,t}\cdot\pi_{r,i,j,t,b}),$

where $\bm{\pi}_{r,i,j,t}$ is equivalent to $\bm{\pi}^{*}_{r,i,j,t}$ but does
not include the probability of failing to detect an individual that was truly
present. Using this conditional likelihood, we do not directly obtain
estimates of the latent abundance $N_{r,i,j,t}$, but we can subsequently draw
these values from a Poisson distribution with mean $\lambda_{r,i,j,t}$.

#### Detection Model

Under the removal sampling protocol, the cell probability of being observed in
the $b$th time interval is defined as

$\pi_{r,i,j,t,b}=p_{r,i,j,t}(1-p_{r,i,j,t})^{b-1},$

where $p_{r,i,j,t}$ is the probability of an individual of species $i$ being
detected in at least one time interval during year $t$ at point $j$ within
park $r$. We incorporated covariate and random effects in the detection
probability models as

$\displaystyle\text{logit}(p_{r,i,j,t})=\alpha_{0,r,i}+\alpha_{1,r,i}\cdot\text{DAY}_{r,j,t}+\alpha_{2,r,i}\cdot\text{DAY}^{2}_{r,j,t}+\alpha_{3,r,i}\cdot\text{TIME}_{r,j,t}+\alpha_{4,r,t},$

where $\alpha_{0,r,i}$ is the intercept for species $i$ in park $r$, and
$\alpha_{1,r,i},\alpha_{2,r,i},\alpha_{3,r,i}$ are regression coefficients
representing the linear effect of Julian date (DAYr,j,t), the quadratic effect
of Julian date ($\text{DAY}^{2}_{r,j,t}$), and the time since sunrise
($\text{TIME}_{r,j,t}$), respectively, on the detection probability of species
$i$ at park $r$. We included a random year effect stratified by park,
$\alpha_{4,r,t}$, to account for variability in detection probability across
years and parks that follows a normal distribution with mean 0 and variance
$\sigma^{2}_{4,p}$. The species/park specific effects
$\alpha_{0,r,i},\alpha_{1,r,i},\alpha_{2,r,i}$, and $\alpha_{3,r,i}$ follow
distributions with common park-specific parameters to enable sharing of
information across species within a given park. Specifically, we model
$\alpha_{1,r,i}$ as

$\alpha_{1,r,i}\sim\text{Normal}(\mu_{p,1,r},\sigma^{2}_{p,\text{REG},1}),$

where $\mu_{p,1,r}$ represents the mean linear effect of day on detection
probability across the entire species community in park $r$, and
$\sigma^{2}_{p,\text{REG},1}$ represents the variability in this effect across
the different species comprising the community in park $r$. Models of
$\alpha_{0,r,i},\alpha_{2,r,i}$, and $\alpha_{3,r,i}$ are defined analogously.
The $\mu_{p,1,r}$ coefficients provide a straightforward way to compare how
covariates influence the bird communities across different parks. Further,
$\mu_{p,1,r}$ follow a common distribution to enable sharing of information
across parks, which allows for estimates in sparsely sampled parks (i.e.,
WEFA, SAGA, SARA). To do this, each $\mu_{p,1,r}$ is modeled as

$\mu_{p,1,r}\sim\text{Normal}(\bar{\mu}_{p,1},\sigma^{2}_{p,\text{META},1}),$

where $\bar{\mu}_{p,1}$ is the mean linear effect of day on detection
probability across all species and parks (i.e., the meta-community), and
$\sigma^{2}_{p,\text{META},1}$ represents the variability in the linear effect
of day on detection probability across parks within the meta-community. Models
of $\mu_{p,0,r},\mu_{p,2,r}$, and $\mu_{p,3,r}$ are defined analogously.

#### Abundance Model

As with the detection model, we incorporated random effects and covariates at
multiple levels of the hierarchy to explain variation in $\lambda_{r,i,j,t}$.
We model $\lambda_{r,i,j,t}$ using the following specification

$\begin{split}\text{log}(\lambda_{r,i,j,t})=\;&\beta_{0,r,i}+\beta_{1,r,i}\cdot\text{YEAR}_{r,t}+\beta_{2,r,i}\cdot\text{REGEN}_{r,j}+\beta_{3,r,i}\cdot\text{FOR}_{r,j}+\\\
&\beta_{4,r,i}\cdot\text{BA}_{r,j}+\beta_{5,r,i}\cdot\text{BA}^{2}_{r,j},\end{split}$

where $\beta_{0,r,i}$ is the abundance intercept for species $i$ in park $r$,
and $\beta_{1,r,i}$, $\beta_{2,r,i}$, $\beta_{3,r,i}$, $\beta_{4,r,i}$, and
$\beta_{5,r,i}$ are regression coefficients representing the linear effect of
year ($\text{YEAR}_{r,t}$), regeneration ($\text{REGEN}_{r,j}$), local forest
cover ($\text{FOR}_{r,j}$), basal area ($\text{BA}_{r,j}$), and quadratic
effect of basal area ($\text{BA}^{2}_{r,j}$), respectively, on the abundance
of species $i$ at park $r$. We define all species and park-specific parameters
as random effects that follow a common distribution within each park. Park-
specific parameters subsequently follow a common meta-community distribution.
For $\beta_{1,r,i}$ we have

$\displaystyle\beta_{1,r,i}\sim\text{Normal}(\mu_{\lambda,1,r},\sigma^{2}_{\lambda,\text{REG},1})$
$\displaystyle\mu_{\lambda,1,r}\sim\text{Normal}(\bar{\mu}_{\lambda,1},\sigma^{2}_{\lambda,\text{META},1}),$

where $\mu_{\lambda,1,r}$ is the mean year effect across the entire community
in park $r$, $\sigma^{2}_{\lambda,\text{REG},1}$ is the variability of the
year effect across all species in the meta-community, $\bar{\mu}_{\lambda,1}$
is the mean year effect across the entire meta-community, and
$\sigma^{2}_{\lambda,\text{META},1}$ is the variability of the park year
effects across the entire meta-community. Models of
$\beta_{0,r,i},\beta_{2,r,i},\beta_{3,r,i},\beta_{4,r,i}$, and $\beta_{5,r,i}$
are defined analogously.

### Bird Guilds

To assess the effects of local forest structure on species’ abundances, we
assigned birds to behavioral and physiological response guilds following
O’Connell et al. (2000). Bird guilds are groups of species that respond in
similar ways to environmental changes as a result of similar uses of the
environment (Root, 1967; Verner, 1984; Szaro, 1986). Overall, we categorized
16 guilds in three biotic integrity elements (Functional, Compositional,
Structural) as specialist or generalist guilds depending on each guild’s
relationship to the landscape structure and function. Specialist guilds
consist of species with highly specific habitat requirements, and generalist
guilds consist of species that can utilize a wide range of habitats. We
selected guilds to indicate different aspects of species’ life history traits,
which can lead to a single species being assigned to multiple guilds (Faccio
et al. 2015; Table 3). Based on the recommendations in Pacifici et al. (2014),
we only used bird guilds in post-hoc assessment.

### Model Estimation

We implemented the model using a Bayesian framework (e.g., Berliner 1996;
Clark 2007; Hobbs and Hooten 2015). The Bayesian approach is particularly
useful in our setting because it easily accommodates missing data and
facilitates inference about derived quantities such as covariate effects
across different bird guilds with propagated uncertainty across all
hierarchical levels (e.g., species and parks). We specified vague normal
priors for regression coefficients and vague gamma priors for variance
parameters. We fit our models using Markov chain Monte Carlo (MCMC) in JAGS
(Plummer, 2003) within the `R` statistical environment (R Core Team, 2019)
using the `jagsUI` (Kellner, 2018) package. We ran models for 50,000
iterations with a burn-in period of 45,000 iterations and a thinning rate of
2. We assessed model convergence using visual assessment of trace plots and
the Gelman-Rubin R-hat diagnostic, where convergence was considered to occur
for all values of R < 1.1 (Brooks and Gelman, 1998; Gelman et al., 2004). We
assessed model fit using a Bayesian p-value approach (Gelman et al., 1996;
Hobbs and Hooten, 2015), which indicated a successful model fit (mean = 0.21).
We interpret a parameter as “significant” if the 95% credible interval does
not include zero. We performed all subsequent analysis in `R` using the `coda`
package (Plummer et al., 2006). All code and data are available on Zenodo
(Doser et al., 2021).

## Results

We observed 106 bird species from 2006-2019 across the eight parks, with an
average of 59 species observed at an individual park (Table 2). SARA had the
highest observed species richness with a total of 71 observed species, while
WEFA had the lowest richness, with a total of 45 observed species. Most
species in each park’s community were rare, as the average point-level
observed abundance of a species across all parks was 0.17, with average point-
level abundance highest at MORR (0.22) and lowest in ACAD (0.11).

### Temporal trends in bird abundance

There was wide variation in species abundance trends among the eight parks
with three parks (MABI, ACAD, MORR) showing significant declines in overall
abundance (across all species), two parks (SAGA, MIMA) showing no trend, and
three parks (ROVA, SARA, WEFA) showing significant increases in abundance over
the time frame of the study (Figure 3b). Park-level species richness showed
similar trends to that of abundance (Figure 4). Interestingly, the wide
variation across parks resulted in no temporal trend in total bird abundance
across the entire network (Figure 3a), revealing the benefit of the
hierarchical modeling approach to evaluating park-specific dynamics that might
be masked at the overall network level (i.e., Simpson’s paradox). To assess if
the variability across spatial units was a result of bird community
composition, we computed the Bird Community Index (BCI; O’Connell et al.
2000), a measure of bird communities that provides inference on the ecological
integrity of a site (See Appendix S1 for details of BCI computation). We
performed a correlation analysis between the estimated year trend and the mean
estimated BCI, which revealed a strong negative correlation (median = -0.77).

Individual species abundance trends were highly variable across parks, but
species within a given park tended to show similar trends (Figures 3c,
Appendix S1: Figure S1). Further, species in the same guild showed similar
trends in abundance within a given park. However, the direction and magnitude
of these trends varied widely across parks (Figure 5).

### Forest structure effects

At the network level, we found moderate support for a peak in bird abundance
at intermediate values of basal area (85% and 78% probability of significant
linear and quadratic effects, respectively), high amounts of local forest
cover (92% probability of a positive effect), and to a lesser extent high
amounts of forest regeneration (70% probability of a positive effect). Most
parks followed the network level pattern, although there was variation among
parks, especially in the effects of percent forest and forest regeneration
(Figure 6). Percent forest has a more prominent effect than either basal area
or forest regeneration, as evidenced by comparison of the magnitudes of the
standardized regression coefficients. There was little variability in the
effects of the covariates across different bird guilds (Appendix S1: Table
S1).

## Discussion

With the growing threats of habitat loss, invasive species, and climate
change, large-scale monitoring networks are becoming increasingly important
for monitoring trends in species distributions and abundance across space and
time. Such large-scale programs are often limited by resources at certain
locations, indicating a need for sophisticated modeling approaches to estimate
species trends across the entire region of interest. Our novel hierarchical
model reveals large variations in bird abundance trends across eight protected
forests but not across bird guilds within a park, suggesting that ecological
processes, biological invasions, and management activities that affect local
forest condition appear to have consistent effects on local forest bird
communities. An understanding of how these variables influence diverse
communities of bird species is crucial to informed and lasting management
solutions for both forests and birds.

Trends in bird abundance differed across space (i.e., parks) but not by
species guilds (Figure 5), suggesting that local forest conditions might have
broad and consistent effects on bird communities within parks. The consistent
trends in species across guilds within a park are surprising given the
variable responses bird guilds show towards insect disturbances (Janousek et
al., 2019), silviculture (Thiollay, 1997; Augenfeld et al., 2008), and
elevation (Tenan et al., 2017). However, life-history characteristics,
including diet, foraging strategy, habitat preference, and nesting location,
do not predict the effects of climate change on bird species distributions in
northeastern North America (Zuckerberg et al., 2009; Cohen et al., 2020),
suggesting responses of bird species to stressors such as climate change might
be independent of life-history traits. The lack of variability we found in
trends across guilds is partially attributable to species being assigned to
multiple guilds. However, this cannot explain the variability in trends of
guilds across the different parks (Figure 5). Thus, our results suggest that
diverse communities of forest birds might show similar responses to variations
in local forest condition and structure, as well as other spatially-dependent
stressors, such as climate change (Zuckerberg et al., 2009; Cohen et al.,
2020).

We identified three parks (MABI, ACAD, MORR) with significant negative
declines in bird abundance, two parks (MIMA, SAGA) with no trend, and three
parks (ROVA, SARA, WEFA) with significant, positive trends (Figure 3).
Understanding the causes of variability in bird abundance trends across parks
and why species within parks behave fairly consistently is critical. Our
findings correspond with previous studies showing variability in bird
abundance and distribution trends across large spatial regions in the U.S.
(Rosenberg et al., 2019; Rushing et al., 2020). The strong negative
correlation between the BCI and the estimated year trend (median = -0.77)
suggests bird communities reflective of higher ecological integrity are
showing the fastest declines over time (Appendix S1: Figure S2). This result
potentially suggests communities of birds within these parks might respond
differently over time as a result of differences in local environmental
stressors and interactions with other species (Hutchinson, 1957), a phenomenon
commonly addressed in species distribution models (Pollock et al., 2014;
Wilkinson et al., 2019; Lany et al., 2020).

While overall declines of bird abundance in MABI, ACAD, and MORR are
concerning, it is important to emphasize that declines in each of these guilds
should be viewed differently. Declines in specialist guilds (e.g., interior
forest obligates, canopy nesters) are of highest concern, as these species are
indicative of bird communities with high ecological integrity (O’Connell et
al., 2000; Ladin et al., 2016). Declines in generalist guilds (e.g., nest
predators/brood parasites, forest generalists) are of lesser concern, and such
declines can actually lead to an increase in the ecological integrity of the
bird community as measured by the BCI. Thus, future analyses and management
efforts should focus on declining specialist guilds in these three parks that
require the interior and older forest habitat these parks are designed to
protect.

The amount of local forest cover has the largest effect on bird abundance
across the parks (Table 4, Figure 6). Four parks (MORR, MIMA, ROVA, SARA)
showed significant positive relationships between bird abundance and percent
forest, which is consistent with our hypothesis and previous findings
(Willson, 1974; Ladin et al., 2016). However, ACAD showed a significant
negative relationship and SAGA showed a non-significant negative relationship,
indicating higher abundance of birds at points surrounded by lower amounts of
forest cover. The non-significant results at SAGA are likely a result of the
low amount of variation in local forest cover (range of 49-70%) near survey
points. However, local forest cover at ACAD varied considerably among points
(51% to 97%). The forests in ACAD and the surrounding forest matrix show
distinct characteristics compared to the other seven parks in terms of their
structural stage, density of large (>30 cm DBH) trees, and tree species
diversity (Miller et al., 2016, 2018)—underscoring that local forest cover
alone does not account for all forest characteristics potentially important to
bird breeding and foraging ecology. Further, it is important to note that
while abundance in ACAD declines with forest cover, these sites were
surrounded by relatively high forest cover, and even the points with the
lowest forest cover are still surrounded by substantial forested habitat that
could be enough habitat to attract and maintain many forest dwelling species
(Willson, 1974; Zuckerberg and Porter, 2010). Overall, these results suggest
management should focus on limiting forest fragmentation and maintaining or
increasing the amount of forested cover in the surrounding landscape matrix,
followed by maintenance of forest structure and diversity. However, the
inconsistent effect of local forest cover across parks suggests that forest
breeding birds are likely affected by interactions between local forest
structure, surrounding land use, local community interactions, climate, and
local stand dynamics (e.g., pest outbreaks, disturbances, succession).

Our finding that overall abundance peaked at intermediate levels of live tree
basal area in five park forests is consistent with previous research
(Flaspohler et al., 2002; Rankin and Perlut, 2015; Vander Yacht et al., 2016).
This effect is likely related to larger variety in vertical forest structure
at intermediate levels that provides suitable habitat to a wider variety of
species (MacArthur and MacArthur, 1961; Crosby et al., 2020). More broadly,
the peaking of abundance at medium levels of live basal area could be
reflective of multi-aged stands supporting a variety of species, individual
tree size, and vertical layering that provides habitat to a wide range of
birds.

Forest regeneration showed the weakest effect of the three covariates, with
only a small positive relationship on overall bird abundance and significant
at only one park (ROVA). Multiple parks had highly positive-skewed
distributions of forest regeneration, which could indicate that there was not
enough regeneration to attract birds that utilize such habitat. Given the
legacy of deer overabundance, invasive species, and regional declines in
regeneration that have shaped the understory in these parks (Miller and
McGill, 2019), more attention to examining the effects of regeneration and
invasive plant abundance on forest birds seems warranted. Additionally, forest
management aimed at increasing regeneration abundance and diversity within
these forests should also be evaluated for its benefits to the overall bird
community.

While there was support for effects of tree basal area, forest regeneration,
and local forest cover on bird abundance, these variables alone cannot explain
the consistent trends of species within parks. ACAD, MABI, and MORR on average
had the most forest within a 1km radius of each point. However, these forests
vary considerably in structural stage (Miller et al., 2016), type, and
diversity (Miller et al., 2018), suggesting that declines in bird abundance
might be driven by factors other than forest structure. MABI is the only park
in the network that is subject to logging, which has led to decreases in basal
area, crown closure, and tree density and increases in regeneration over the
study period that could contribute to the strong decreasing trends in bird
abundance. MORR has experienced a decrease in crown closure as a result of
intense storms and the invasive Emerald Ash Borer (_Agrilus planipennis_),
leading to increased gaps in the forests that could contribute to the
decreasing trend in overall bird abundance. Reasons for the declines in ACAD
are more ambiguous, but they could partially be attributed to declines in
early successional habitat throughout the park. Overall, these variations in
trends across spatial units are likely a composite result of numerous effects
such as agricultural intensification and land use change outside of the parks
(Pino et al., 2000), variations in forest age, structure, and size (Miller et
al., 2016), differing levels of anthropogenic mortality, loss of sensitive
areas or resources, migratory behavior, and climate change effects occurring
during both the breeding and nonbreeding seasons (Rushing et al., 2020).

Because the point counts did not contain specific information on the spatial
location of each individual bird, we chose to describe population sizes in
terms of point-level abundances rather than a measure of density, as the true
effective sampling area of the point count is not well-defined due to movement
and temporary immigration/emigration of birds (Kéry and Royle, 2015). This
inability to accurately describe the sampling area is an inherent problem of
any non-spatial method of counting birds (e.g., repeated counts, removal
sampling) and is only remedied via approaches that incorporate spatial
information of each individual, such as distance sampling (Buckland et al.,
1993) and spatial capture-recapture (Efford, 2004; Borchers and Efford, 2008;
Royle and Young, 2008). However, inference on point-level abundances, as used
in the current study, provides information on population dynamics in time and
space, which can yield management-relevant insights (Yamaura et al., 2011;
Chandler et al., 2013).

We developed a novel hierarchical model that enables inference on abundance
for all 106 species observed across the eight parks and provides broad
inference on how forest conditions affect bird species. We extended previous
work on community abundance models (Yamaura et al., 2011, 2012; Chandler et
al., 2013) to efficiently model communities of species across multiple
geographically distinct regions. Our flexible modeling approach is widely
applicable to monitoring programs where inference is desired on abundance of
rare species and abundance of species at sparsely sampled locations. Further,
our modeling approach can be easily extended to other sampling protocols
commonly used to estimate abundance (e.g., distance sampling). Many of the 106
species in our data set are rare (i.e., lots of zeros in data), which leads to
a trade-off between the number of species for which we can estimate trends and
the number of covariates we can include in the analysis. This, in addition to
the temporal misalignment of the forest covariates, limited our ability to use
time-varying covariates to relate variability in forest structure over time to
bird abundance trends. The addition of variables reflecting changes in species
or vertical structure diversity (e.g., Shannon diversity, Gini diversity), as
well as other variables related to climate change could likely provide
additional insight on our results. Future work might also seek to combine
additional data sources (e.g., Breeding Bird Survey, acoustic recordings,
survival/productivity data) with our data using an integrated modeling
approach (e.g., Miller et al. 2019) to help address this trade off, provide
further inference on what is driving these trends in specialist guilds
indicative of high integrity bird communities, and provide a more mechanistic
approach at assessing population declines (Zipkin and Saunders, 2018; Saunders
et al., 2019).

## Acknowledgements

JWD, ASW, EFZ, and AOF conceived the modeling approach; ASW and KMM managed
and organized the data collection; JWD analyzed the data and led writing of
the manuscript. All authors contributed critically to the drafts and gave
final approval for publication. We declare that there are no conflicts of
interest. We are grateful to Camilla Seirup and the field crews that
contributed to the collection of NETN forest data in this study, to the
volunteers, Steve Faccio, and Ed Sharron for their contributions to the NETN
bird monitoring program, and to Adam Kozlowski and Brian Mitchell for their
contributions to both programs. We thank John Marzluff, Marc Kéry, and an
anonymous reviewer for insightful comments that greatly improved the
manuscript. This work was supported by the National Park Service’s Inventory
and Monitoring Division and National Science Foundation grants DMS-1916395,
EF-1253225, and DBI-1954406.

## Open Research

Data and code (Doser et al., 2021) are available on Zenodo at:
https://doi.org/10.5281/zenodo.4701477.

## References

  * Amundson et al. (2014) Amundson, C. L., Royle, J. A., and Handel, C. M. (2014).  A hierarchical model combining distance sampling and time removal to estimate detection probability during avian point counts.  The Auk, 131(4):476–494.
  * Augenfeld et al. (2008) Augenfeld, K. H., Franklin, S. B., and Snyder, D. H. (2008).  Breeding bird communities of upland hardwood forest 12 years after shelterwood logging.  Forest Ecology and Management, 255(3-4):1271–1282.
  * Augustine and DeCalesta (2003) Augustine, D. J. and DeCalesta, D. (2003).  Defining deer overabundance and threats to forest communities: From individual plants to landscape structure.  Ecoscience, 10(4):472–486.
  * Berliner (1996) Berliner, M. L. (1996).  Hierarchical bayesian time series models.  Maximum entropy and Bayesian methods, 79.
  * Bibby (1999) Bibby, C. J. (1999).  Making the most of birds as environmental indicators.  Ostrich, 70(1):81–88.
  * Borchers and Efford (2008) Borchers, D. L. and Efford, M. G. (2008).  Spatially explicit maximum likelihood methods for capture–recapture studies.  Biometrics, 64(2):377–385.
  * Brooks and Gelman (1998) Brooks, S. P. and Gelman, A. (1998).  General methods for monitoring convergence of iterative simulations.  Journal of Computational and Graphical Statistics, 7(4):434–455.
  * Bruns (1960) Bruns, H. (1960).  The Economic Importance of Birds in Forests.  Bird Study, 7(4):193–208.
  * Buckland et al. (1993) Buckland, S. T., Anderson, D. R., Burnham, K., and Laake, J. (1993).  Distance sampling: estimating abundance of biological populations.  Oxford University Press.
  * Canterbury et al. (2000) Canterbury, G. E., Martin, T. E., Petit, D. R., Petit, L. J., and Bradford, D. F. (2000).  Bird communities and habitat as ecological indicators of forest condition in regional monitoring.  Conservation Biology, 14(2):544–558.
  * Chandler et al. (2013) Chandler, R. B., King, D. I., Raudales, R., Trubey, R., Chandler, C., and Arce Chávez, V. J. (2013).  A small-scale land-sparing approach to conserving biological diversity in tropical agricultural landscapes.  Conservation Biology, 27(4):785–795.
  * Clark (2007) Clark, J. S. (2007).  Models for Ecological Data.  Princeton University Press.
  * Cohen et al. (2020) Cohen, J. M., Fink, D., and Zuckerberg, B. (2020).  Avian responses to extreme weather across functional traits and temporal scales.  Global Change Biology, (March):1–11.
  * Crosby et al. (2020) Crosby, A. D., Porter, W. F., Roloff, G. J., Walters, M. B., and Donovan, M. L. (2020).  Combining conservation value with conservation filters to guide forest management for avian biodiversity.  Forest Ecology and Management, 466(January).
  * Dail and Madsen (2011) Dail, D. and Madsen, L. (2011).  Models for Estimating Abundance from Repeated Counts of an Open Metapopulation.  Biometrics, 67(2):577–587.
  * Dorazio et al. (2005) Dorazio, R. M., Jelks, H. L., and Jordan, F. (2005).  Improving Removal-Based Estimates of Abundance by Sampling a Population of Spatially Distinct Subpopulations.  Biometrics, 61(4):1093–1101.
  * Dorazio and Royle (2005) Dorazio, R. M. and Royle, J. A. (2005).  Estimating size and composition of biological communities by modeling the occurrence of species.  Journal of the American Statistical Association, 100(470):389–398.
  * Dorazio et al. (2006) Dorazio, R. M., Royle, J. A., Söderström, B., and Glimskär, A. (2006).  Estimating species richness and accumulation by modeling species occurrence and detectability.  Ecology, 87(4):842–854.
  * Doser et al. (2021) Doser, J. W., Weed, A. S., Zipkin, E. F., Miller, K. M., and Finley, A. O. (2021).  Code and data for Trends in bird abundance differ among protected forests but not bird guilds Zenodo, https://doi.org/10.5281/zenodo.4701477.
  * Efford (2004) Efford, M. (2004).  Density estimation in live-trapping studies.  Oikos, 106(3):598–610.
  * Faccio et al. (2015) Faccio, S., Mitchell, B. R., and Pooler, P. S. (2015).  Northeast Temperate Network breeding landbird monitoring protocol: 2015 revision.  Technical Report Natural Resource Report NPS/NETN/NRR-2015/942, National Park Service, Fort Collins, Colorado.
  * Farnsworth et al. (2002) Farnsworth, G. L., Pollock, K. H., Nichols, J. D., Simons, T. R., Hines, J. E., and Sauer, J. R. (2002).  A Removal Model for Estimating Detection Probabilities from Point-Count Surveys.  The Auk, 119(2):414–425.
  * Farr et al. (2019) Farr, M. T., Green, D. S., Holekamp, K. E., Roloff, G. J., and Zipkin, E. F. (2019).  Multispecies hierarchical modeling reveals variable responses of African carnivores to management alternatives.  Ecological Applications, 29(2):1–11.
  * Flaspohler et al. (2002) Flaspohler, D. J., Fisher-Huckins, C. J., Bub, B. R., and Dusen, P. J. V. (2002).  Temporal patterns in aquatic and avian communities following selective logging in the upper Great Lakes region.  Forest Science, 48(2):339–349.
  * Gelfand et al. (2005) Gelfand, A. E., Schmidt, A. M., Wu, S., Silander Jr, J. A., Latimer, A., and Rebelo, A. G. (2005).  Modelling species diversity through species level hierarchical modelling.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 54(1):1–20.
  * Gelfand et al. (2006) Gelfand, A. E., Silander, J. A., Wu, S., Latimer, A., Lewis, P. O., Rebelo, A. G., Holder, M., et al. (2006).  Explaining species distribution patterns through hierarchical modeling.  Bayesian Analysis, 1(1):41–92.
  * Gelman et al. (2004) Gelman, A., Carlin, J. B., Stern, H. S., and Rubin, D. B. (2004).  Bayesian Data Analysis.  Chapman and Hall/CRC, 2nd ed. edition.
  * Gelman et al. (1996) Gelman, A., Meng, X. L., and Stern, H. (1996).  Posterior predictive assessment of model fitness via realized discrepancies.  Statistica Sinica, 6(4):733–807.
  * Gregory et al. (2003) Gregory, R. D., Noble, D., Field, R., Marchant, J., Raven, M., and Gibbons, D. W. (2003).  Using birds as indicators of biodiversity.  Ornis Hungarica, 12(13):11–24.
  * Guisan and Thuiller (2005) Guisan, A. and Thuiller, W. (2005).  Predicting species distribution: Offering more than simple habitat models.  Ecology Letters, 8(9):993–1009.
  * Hobbs and Hooten (2015) Hobbs, N. T. and Hooten, M. B. (2015).  Bayesian models: a statistical primer for ecologists.  Princeton University Press.
  * Homer et al. (2015) Homer, C., Dewitz, J., Yang, L., Jin, S., Danielson, P., Xian, G., Coulston, J., Herold, N., Wickham, J., and Megown, K. (2015).  Completion of the 2011 national land cover database for the conterminous united states–representing a decade of land cover change information.  Photogrammetric Engineering & Remote Sensing, 81(5):345–354.
  * Hutchinson (1957) Hutchinson, G. E. (1957).  Concluding remarks. Cold Spring Harbor Symposia on Quantitative Biology.  Quantitative Biology, (22):415–427.
  * Inger et al. (2015) Inger, R., Gregory, R., Duffy, J. P., Stott, I., Voříšek, P., and Gaston, K. J. (2015).  Common European birds are declining rapidly while less abundant species’ numbers are rising.  Ecology Letters, 18(1):28–36.
  * Janousek et al. (2019) Janousek, W. M., Hicke, J. A., Meddens, A. J., and Dreitz, V. J. (2019).  The effects of mountain pine beetle outbreaks on avian communities in lodgepole pine forests across the greater Rocky Mountain region.  Forest Ecology and Management, 444(April):374–381.
  * Karr and Roth (1971) Karr, J. R. and Roth, R. R. (1971).  Vegetation Structure and Avian Diversity in Several New World Areas.  The American Naturalist, 105(945):423–435.
  * Kellner (2018) Kellner, K. (2018).  jagsUI: A Wrapper Around rjags to Streamline JAGS Analyses.  R package version 1.5.0.
  * Kéry and Royle (2015) Kéry, M. and Royle, J. A. (2015).  Applied Hierarchical Modeling in Ecology: Analysis of distribution, abundance and species richness in R and BUGS: Volume 1: Prelude and Static Models.  Academic Press.
  * Ladin et al. (2016) Ladin, Z. S., Higgins, C. D., Schmit, J. P., Sanders, G., Johnson, M. J., Weed, A. S., Marshall, M. R., Campbell, J. P., Comiskey, J. A., and Shriver, W. G. (2016).  Using regional bird community dynamics to evaluate ecological integrity within national parks.  Ecosphere, 7(9).
  * Lany et al. (2020) Lany, N. K., Zarnetske, P. L., Finley, A. O., and McCullough, D. G. (2020).  Complementary strengths of spatially-explicit and multi-species distribution models.  Ecography, 43(3):456–466.
  * Linden and Roloff (2013) Linden, D. W. and Roloff, G. J. (2013).  Retained structures and bird communities in clearcut forests of the Pacific Northwest, USA.  Forest Ecology and Management, 310:1045–1056.
  * MacArthur et al. (1966) MacArthur, R., Recher, H., and Cody, M. (1966).  On the Relation between Habitat Selection and Species Diversity.  The American Naturalist, 100(913):319–332.
  * MacArthur (1958) MacArthur, R. H. (1958).  Population ecology of some warblers of northeastern coniferous forests.  Ecology, 39(4):599–619.
  * MacArthur (1964) MacArthur, R. H. (1964).  Environmental Factors Affecting Bird Species Diversity.  The American Naturalist, 98(903):387–397.
  * MacArthur and MacArthur (1961) MacArthur, R. H. and MacArthur, J. W. (1961).  On Bird Species Diversity.  Ecology, 42(3):594–598.
  * Miller et al. (2019) Miller, D. A., Pacifici, K., Sanderlin, J. S., and Reich, B. J. (2019).  The recent past and promising future for data integration methods to estimate species’ distributions.  Methods in Ecology and Evolution, 10(1):22–37.
  * Miller et al. (2016) Miller, K. M., Dieffenbach, F. W., Campbell, J. P., Cass, W. B., Comiskey, J. A., Matthews, E. R., McGill, B. J., Mitchell, B. R., Perles, S. J., Sanders, S., Schmit, J. P., Smith, S., and Weed, A. S. (2016).  National parks in the eastern United States harbor important older forest structure compared with matrix forests.  Ecosphere, 7(7):1–20.
  * Miller and McGill (2019) Miller, K. M. and McGill, B. J. (2019).  Compounding human stressors cause major regeneration debt in over half of eastern US forests.  Journal of Applied Ecology, 56(6):1355–1366.
  * Miller et al. (2018) Miller, K. M., McGill, B. J., Mitchell, B. R., Comiskey, J., Dieffenbach, F. W., Matthews, E. R., Perles, S. J., Schmit, J. P., and Weed, A. S. (2018).  Eastern national parks protect greater tree species diversity than unprotected matrix forests.  Forest Ecology and Management, 414:74–84.
  * Nichols et al. (2009) Nichols, J. D., Thomas, L., and Conn, P. B. (2009).  Inferences about landbird abundance from count data: Recent advances and future directions. In Environmental and Ecological Statistics, volume 3.  Springer.
  * O’Connell et al. (2000) O’Connell, T. J., Jackson, L. E., and Brooks, R. P. (2000).  Bird Guilds as Indicators of Ecological Condition in the Central Appalachians.  Ecological Applications, 10(6):1706–1721.
  * Otis et al. (1981) Otis, D. L., Burnham, K. P., White, G. C., and Anderson, D. R. (1981).  Statistical Inference From Capture Data on Closed Animal Populations.  Wildlife Monographs, 62:3–135.
  * Ovaskainen et al. (2010) Ovaskainen, O., Hottola, J., and Siitonen, J. (2010).  Modeling species co-occurrence by multivariate logistic regression generates new hypotheses on fungal interactions.  Ecology, 91(9):2514–2521.
  * Ovaskainen and Soininen (2011) Ovaskainen, O. and Soininen, J. (2011).  Making more out of sparse data: hierarchical modeling of species communities.  Ecology, 92(2):289–295.
  * Pacifici et al. (2014) Pacifici, K., Zipkin, E. F., Collazo, J. A., Irizarry, J. I., and Dewan, A. (2014).  Guidelines for a priori grouping of species in hierarchical community models.  Ecology and Evolution, 4(7):877–888.
  * Pino et al. (2000) Pino, J., Rodà, F., Ribas, J., and Pons, X. (2000).  Landscape structure and bird species richness: Implications for conservation in rural areas between natural parks.  Landscape and Urban Planning, 49(1-2):35–48.
  * Plummer (2003) Plummer, M. (2003).  Jags: A program for analysis of bayesian graphical models using gibbs sampling.
  * Plummer et al. (2006) Plummer, M., Best, N., Cowles, K., and Vines, K. (2006).  Coda: Convergence diagnosis and output analysis for mcmc.  R News, 6(1):7–11.
  * Pollock et al. (2014) Pollock, L. J., Tingley, R., Morris, W. K., Golding, N., O’Hara, R. B., Parris, K. M., Vesk, P. A., and Mccarthy, M. A. (2014).  Understanding co-occurrence by modelling species simultaneously with a Joint Species Distribution Model (JSDM).  Methods in Ecology and Evolution, 5(5):397–406.
  * R Core Team (2019) R Core Team (2019).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Rankin and Perlut (2015) Rankin, D. T. and Perlut, N. G. (2015).  The effects of Forest Stand Improvement Practices on occupancy and abundance of breeding songbirds.  Forest Ecology and Management, 335:99–107.
  * Root (1967) Root, R. B. (1967).  The niche exploitation pattern of the blue-gray gnatcatcher.  Ecological monographs, 37(4):317–350.
  * Rosenberg et al. (2019) Rosenberg, K. V., Dokter, A. M., Blancher, P. J., Sauer, J. R., Smith, A. C., Smith, P. A., Stanton, J. C., Panjabi, A., Helft, L., Parr, M., et al. (2019).  Decline of the north american avifauna.  Science, 366(6461):120–124.
  * Royle (2004a) Royle, J. A. (2004a).  Generalized estimators of avian abundance from count survey data.  Animal Biodiversity and Conservation, 27(1):375–386.
  * Royle (2004b) Royle, J. A. (2004b).  N-Mixture Models for Estimating Population Size from Spatially Replicated Counts.  Biometrics, 60(1):108–115.
  * Royle (2009) Royle, J. A. (2009).  Analysis of Capture-Recapture Models with Individual Covariates Using Data Augmentation.  Biometrics, 65(1):267–274.
  * Royle and Dorazio (2006) Royle, J. A. and Dorazio, R. M. (2006).  Hierarchical models of animal abundance and occurrence.  Journal of Agricultural, Biological, and Environmental Statistics, 11(3):249–263.
  * Royle and Dorazio (2008) Royle, J. A. and Dorazio, R. M. (2008).  Hierarchical modeling and inference in ecology: the analysis of data from populations, metapopulations and communities.  Elsevier.
  * Royle and Young (2008) Royle, J. A. and Young, K. V. (2008).  A hierarchical model for spatial capture–recapture data.  Ecology, 89(8):2281–2289.
  * Rushing et al. (2020) Rushing, C. S., Rohrbaugh, R. W., Fiss, C. J., Rosenberry, C. S., Rodewald, A. D., and Larkin, J. L. (2020).  Long-term variation in white-tailed deer abundance shapes landscape-scale population dynamics of forest-breeding birds.  Forest Ecology and Management, 456:117629.
  * Russell et al. (2001) Russell, F. L., Zippin, D. B., and Fowler, N. L. (2001).  Effects of White-Tailed Deer ( Odocoileus Virginianus ) on Plants , Plant Populations and Communities : A Review.  The American Midland Naturalist, 146(1):1–26.
  * Saunders et al. (2019) Saunders, S. P., Farr, M. T., Wright, A. D., Bahlai, C. A., Ribeiro, J. W., Rossman, S., Sussman, A. L., Arnold, T. W., and Zipkin, E. F. (2019).  Disentangling data discrepancies with integrated population models.  Ecology, 100(6):1–14.
  * Stanton et al. (2018) Stanton, R. L., Morrissey, C. A., and Clark, R. G. (2018).  Analysis of trends and agricultural drivers of farmland bird declines in North America: A review.  Agriculture, Ecosystems and Environment, 254(December 2017):244–254.
  * Stem et al. (2005) Stem, C., Margoluis, R., Salafsky, N., and Brown, M. (2005).  Monitoring and Evaluation in Conservation: a Review of Trends and Approaches.  Conservation Biology, 19(2):295–309.
  * Stevens (1997) Stevens, D. L. (1997).  Variable density grid-based sampling designs for continuous spatial populations.  Environmetrics, 8(3):167–195.
  * Stevens and Olsen (2004) Stevens, D. L. and Olsen, A. R. (2004).  Spatially balanced sampling of natural resources.  Journal of the American Statistical Association, 99(465):262–278.
  * Sutherland et al. (2016) Sutherland, C., Brambilla, M., Pedrini, P., and Tenan, S. (2016).  A multiregion community model for inference about geographic variation in species richness.  Methods in Ecology and Evolution, 7(7):783–791.
  * Szaro (1986) Szaro, R. C. (1986).  Guild management: an evaluation of avian guilds as a predictive tool.  Environmental Management, 10(5):681–688.
  * Tenan et al. (2017) Tenan, S., Brambilla, M., Pedrini, P., and Sutherland, C. (2017).  Quantifying spatial variation in the size and structure of ecologically stratified communities.  Methods in Ecology and Evolution, 8(8):976–984.
  * Thiollay (1997) Thiollay, J. M. (1997).  Disturbance, selective logging and bird diversity: A neotropical forest study.  Biodiversity and Conservation, 6(8):1155–1173.
  * Tierney et al. (2016) Tierney, G., Mitchell, B., Miller, K., Comiskey, J., A., K., and D., F.-L. (2016).  Northeast Temperate Network long-term forest monitoring protocol: 2016 revision.  Technical Report Natural Resource Report NPS/NETN/NRR-2016/1184, National Park Service, Fort Collins, Colorado.
  * Vander Yacht et al. (2016) Vander Yacht, A. L., Keyser, P. D., Buehler, D. A., Harper, C. A., Buckley, D. S., and Applegate, R. D. (2016).  Avian occupancy response to oak woodland and savanna restoration.  Journal of Wildlife Management, 80(6):1091–1105.
  * Verner (1984) Verner, J. (1984).  The guild concept applied to management of bird populations.  Environmental Management, 8(1):1–13.
  * Warton et al. (2015) Warton, D. I., Blanchet, F. G., O’Hara, R. B., Ovaskainen, O., Taskinen, S., Walker, S. C., and Hui, F. K. (2015).  So many variables: joint modeling in community ecology.  Trends in Ecology & Evolution, 30(12):766–779.
  * Whelan et al. (2010) Whelan, C., Wenny, D., and Marquis, R. (2010).  Policy implications of ecosystem services provided by birds.  Synesis: A Journal of Science, Technology, Ethics, and Policy, 1(1):T11—-T20.
  * Whelan et al. (2008) Whelan, C. J., Wenny, D. G., and Marquis, R. J. (2008).  Ecosystem services provided by birds.  Annals of the New York Academy of Sciences, 1134(1):25–60.
  * Whelan et al. (2015) Whelan, C. J., Şekercioğlu, Ç. H., and Wenny, D. G. (2015).  Why birds matter: from economic ornithology to ecosystem services.  Journal of Ornithology, 156:227–238.
  * Wilkinson et al. (2019) Wilkinson, D. P., Golding, N., Guillera-Arroita, G., Tingley, R., and McCarthy, M. A. (2019).  A comparison of joint species distribution models for presence–absence data.  Methods in Ecology and Evolution, 10(2):198–211.
  * Willson (1974) Willson, M. F. (1974).  Avian Community Organization and Habitat Structure.  Ecology, 55(5):1017–1029.
  * Wright et al. (2020) Wright, A. D., Grant, E. H. C., and Zipkin, E. F. (2020).  A hierarchical analysis of habitat area, connectivity, and quality on amphibian diversity across spatial scales.  Landscape Ecology.
  * Yackulic et al. (2020) Yackulic, C. B., Dodrill, M., Dzul, M., Sanderlin, J. S., and Reid, J. A. (2020).  A need for speed in bayesian population models: a practical guide to marginalizing and recovering discrete latent states.  Ecological Applications, 30(5):e02112.
  * Yamaura et al. (2011) Yamaura, Y., Andrew Royle, J., Kuboi, K., Tada, T., Ikeno, S., and Makino, S. (2011).  Modelling community dynamics based on species-level abundance models from detection/nondetection data.  Journal of Applied Ecology, 48(1):67–75.
  * Yamaura et al. (2012) Yamaura, Y., Royle, J. A., Shimada, N., Asanuma, S., Sato, T., Taki, H., and Makino, S. (2012).  Biodiversity of man-made open habitats in an underused country: a class of multispecies abundance models for count data.  Biodiversity and Conservation, 21(6):1365–1380.
  * Zhao and Royle (2019) Zhao, Q. and Royle, J. A. (2019).  Dynamic N-mixture models with temporal variability in detection probability.  Ecological Modelling, 393(November 2018):20–24.
  * Zipkin et al. (2010) Zipkin, E. F., Andrew Royle, J., Dawson, D. K., and Bates, S. (2010).  Multi-species occurrence models to evaluate the effects of conservation and management actions.  Biological Conservation, 143(2):479–484.
  * Zipkin et al. (2009) Zipkin, E. F., Dewan, A., and Andrew Royle, J. (2009).  Impacts of forest fragmentation on species richness: A hierarchical approach to community modelling.  Journal of Applied Ecology, 46(4):815–822.
  * Zipkin et al. (2020) Zipkin, E. F., DiRenzo, G. V., Ray, J. M., Rossman, S., and Lips, K. R. (2020).  Tropical snake diversity collapses after widespread amphibian loss.  Science, 367(6479):814–816.
  * Zipkin and Saunders (2018) Zipkin, E. F. and Saunders, S. P. (2018).  Synthesizing multiple data types for biological conservation using integrated population models.  Biological Conservation, 217(March 2017):240–250.
  * Zuckerberg and Porter (2010) Zuckerberg, B. and Porter, W. F. (2010).  Thresholds in the long-term responses of breeding birds to forest cover and fragmentation.  Biological Conservation, 143(4):952–962.
  * Zuckerberg et al. (2009) Zuckerberg, B., Woods, A. M., and Porter, W. F. (2009).  Poleward shifts in breeding bird distributions in New York State.  Global Change Biology, 15(8):1866–1883.

## Tables

Table 1: General information on NETN parks included in the landbird monitoring program. Park Name | Code | Forested Area | Years Sampled | Points/Year
---|---|---|---|---
|  | (Ha) | (Year Established) | Mean (sd)
Acadia | ACAD | 8178 | 13 (2007) | 35.8 (9.6)
Saratoga | SARA | 687 | 12 (2007) | 13.3 (6.5)
Morristown | MORR | 626 | 14 (2006) | 20.6 (6.4)
Roosevelt-Vanderbilt | ROVA | 338 | 14 (2007) | 24.4 (6.1)
Minute Man | MIMA | 234 | 14 (2006) | 21.5 (3.0)
Marsh-Billings-Rockefeller | MABI | 196 | 14 (2006) | 23.6 (2.3)
Saint-Gaudens | SAGA | 48 | 12 (2006) | 7.9 (0.3)
Weir Farm | WEFA | 18 | 11 (2006) | 5 (0)

Table 2: Summary statistics of the observed bird community at each park in the NETN. Species abundance is the observed number of individuals of a given species at a given point. Park | Observed Species | Average Species Abundance | Species Abundance Range
---|---|---|---
|  | (birds/point) | (birds/point)
ACAD | 69 | 0.11 | [0.0021, 1.43]
MABI | 58 | 0.20 | [0.0030, 2.95]
MIMA | 60 | 0.22 | [0.0033, 1.46]
MORR | 57 | 0.22 | [0.0035, 1.31]
ROVA | 63 | 0.15 | [0.0029, 0.83]
SAGA | 49 | 0.21 | [0.011, 1.79]
SARA | 71 | 0.20 | [0.0063, 1.44]
WEFA | 45 | 0.18 | [0.018, 1.16]
NETN | 106 | 0.17 | [0.0021, 2.95]

Table 3: Information on bird guilds used in hierarchical model. Adopted from O’Connell et al. (2000). Biotic Integrity Element | Response Guild | Type | Number of Species
---|---|---|---
Functional | Omnivore | Generalist | 34
Functional | Bark prober | Specialist | 10
Functional | Ground gleaner | Specialist | 7
Functional | Upper canopy forager | Specialist | 11
Functional | Lower canopy forager | Specialist | 20
Compositional | Exotic | Generalist | 4
Compositional | Resident | Generalist | 29
Compositional | Single-brooded | Specialist | 65
Compositional | Nest predator/brood parasite | Generalist | 7
Compositional | Temperate migrant | Generalist | 26
Structural | Canopy nester | Specialist | 31
Structural | Shrub nester | Generalist | 20
Structural | Forest-ground nester | Specialist | 14
Structural | Interior forest obligate | Specialist | 29
Structural | Forest generalist | Generalist | 25
Structural | Open-ground nester | Specialist | 9

Table 4: Park-specific posterior medians of covariate effects on abundance and detection. Boldface indicates significance. Park |  |  | Abundance |  |  |  | Detection |
---|---|---|---|---|---|---|---|---
| BA |  BA2 | %Forest | Regen | Year | Day |  Day2 | Time
MABI | -0.16 | 0.075 | 0.050 | 0.093 | -0.18 | 0.0013 | -0.087 | -0.16
ACAD | 0.042 | -0.15 | -0.25 | -0.046 | -0.18 | 0.091 | -0.02 | -0.12
MORR | -0.28 | -0.24 | 0.54 | -0.077 | -0.11 | 0.44 | -0.12 | -0.032
SAGA | 0.009 | -0.12 | -0.15 | 0.11 | -0.041 | 0.093 | -0.012 | 0.14
MIMA | -0.16 | -0.069 | 0.45 | 0.025 | 0.019 | 0.11 | -0.041 | -0.22
ROVA | -0.006 | 0.091 | 0.91 | 0.18 | 0.12 | 0.16 | -0.077 | 0.31
SARA | 0.037 | -0.071 | 0.29 | -0.046 | 0.14 | -0.078 | -0.024 | -0.071
WEFA | -0.094 | -0.02 | 0.12 | 0.23 | 0.24 | -0.12 | -0.0045 | 0.035
NETN | -0.074 | -0.064 | 0.26 | 0.047 | 0.0014 | 0.088 | -0.047 | -0.012

## Figure Legends

Figure 1: Location of NETN parks participating in the breeding landbird
monitoring program in the northeastern United States. See Table 1 for
description of park codes.

Figure 2: Overview of hierarchical model components. Observations
$\bm{y}_{r,i,j,t}$ are obtained for species $i$ in park $r$ at point $j$ in
year $t$ using traditional point count surveys. Species level detection
($\bm{\alpha}_{r,i}$) and abundance ($\bm{\beta}_{r,i}$) parameters follow
common park-level distributions with park level coefficients
($\bm{\mu}_{p,r}$, $\bm{\mu}_{\lambda,r}$), which subsequently follow common
network level distributions with network level coefficients
$\bar{\bm{\mu}}_{p}$ and $\bar{\bm{\mu}}_{\lambda}$, respectively.

Figure 3: Trends in bird abundance from 2006-2019. Panel (a) shows the median
linear trend of year across the entire network with the 95% credible interval
in parentheses. The trend estimates represent the average change in abundance
on the log scale of a species per point every four years. Panel (b) shows the
park level trends. Red highlight indicates a significant negative year trend,
white highlight indicates no significant trend, and blue highlight indicates
positive significant trend. Panel (c) shows the number of species with median
linear trend of year estimates that are negative (red) and positive (blue)
within each park. The number of species with significant trends is shown in
boldface in parentheses.

Figure 4: Trends in estimated park-level species richness from 2006-2019 at
each of the eight parks. Points are the mean species richness across all
points sampled in each given park, gray regions denote the 95% credible
intervals. Inset text is the median (95% credible interval) linear trend
estimate for a post-hoc analysis between park-level species richness and year,
representing the estimated change per year in species richness per park.

Figure 5: Mean linear year effect at each park averaged across all species
defined within 16 bird guilds. * denotes significance. Gray color indicates no
species in that guild were detected.

Figure 6: Relationship between mean expected number of birds of a single
species per point and basal area (a), percent forest in a 1km radius (b), and
forest regeneration (c). Abundance estimates in each graph are estimated at
the median levels of all other covariates within a given park.

## Figures

Figure 1:

Figure 2:

Figure 3:

Figure 4:

Figure 5:

Figure 6: