Sharing detection heterogeneity information among species in community models of occupancy and abundance can strengthen inference

Abstract The estimation of abundance and distribution and factors governing patterns in these parameters is central to the field of ecology. The continued development of hierarchical models that best utilize available information to inform these processes is a key goal of quantitative ecologists. However, much remains to be learned about simultaneously modeling true abundance, presence, and trajectories of ecological communities. Simultaneous modeling of the population dynamics of multiple species provides an interesting mechanism to examine patterns in community processes and, as we emphasize herein, to improve species‐specific estimates by leveraging detection information among species. Here, we demonstrate a simple but effective approach to share information about observation parameters among species in hierarchical community abundance and occupancy models, where we use shared random effects among species to account for spatiotemporal heterogeneity in detection probability. We demonstrate the efficacy of our modeling approach using simulated abundance data, where we recover well our simulated parameters using N‐mixture models. Our approach substantially increases precision in estimates of abundance compared with models that do not share detection information among species. We then expand this model and apply it to repeated detection/non‐detection data collected on six species of tits (Paridae) breeding at 119 1 km2 sampling sites across a P. montanus hybrid zone in northern Switzerland (2004–2020). We find strong impacts of forest cover and elevation on population persistence and colonization in all species. We also demonstrate evidence for interspecific competition on population persistence and colonization probabilities, where the presence of marsh tits reduces population persistence and colonization probability of sympatric willow tits, potentially decreasing gene flow among willow tit subspecies. While conceptually simple, our results have important implications for the future modeling of population abundance, colonization, persistence, and trajectories in community frameworks. We suggest potential extensions of our modeling in this paper and discuss how leveraging data from multiple species can improve model performance and sharpen ecological inference.

1. The estimation of abundance and distribution and factors governing patterns in these parameters is central to the field of ecology. The continued development of hierarchical models that best utilize available information to inform these processes is a key goal of quantitative ecologists. However, much remains to be learned about simultaneously modeling true abundance, presence, and trajectories of ecological communities.
2. Simultaneous modeling of the population dynamics of multiple species provides an interesting mechanism to examine patterns in community processes and, as we emphasize herein, to improve species-specific estimates by leveraging detection information among species. Here, we demonstrate a simple but effective approach to share information about observation parameters among species in hierarchical community abundance and occupancy models, where we use shared random effects among species to account for spatiotemporal heterogeneity in detection probability.
3. We demonstrate the efficacy of our modeling approach using simulated abundance data, where we recover well our simulated parameters using N-mixture models. Our approach substantially increases precision in estimates of abundance compared with models that do not share detection information among species. We then expand this model and apply it to repeated detection/nondetection data collected on six species of tits (Paridae) breeding at 119 1 km 2 sampling sites across a P. montanus hybrid zone in northern Switzerland (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020). We find strong impacts of forest cover and elevation on population persistence and colonization in all species. We also demonstrate evidence for interspecific competition on population persistence and colonization probabilities, where the presence of marsh tits reduces population persistence and colonization probability of sympatric willow tits, potentially decreasing gene flow among willow tit subspecies.

| INTRODUC TI ON
An understanding of changes in population abundance, persistence, and colonization and the demographic processes that drive these changes is a key focus of ecology and many related fields such as conservation biology and wildlife management. Further, scaling up local demographic processes and traits to understand range-wide population changes is a fundamental question in population and community ecology (Sutherland et al., 2013). Thus, the accurate estimation of abundance and its demographic drivers is critical. The increasing accessibility of complex ecological models (e.g., Kéry & Royle, 2021;Schaub & Kéry, 2021) and "big data," and the introduction of novel quantitative approaches (e.g., Besbeas et al., 2002;Dail & Madsen, 2011;Royle, 2004;Zhao et al., 2017Zhao et al., , 2019 have allowed ecologists to begin to estimate these parameters more widely and efficiently. The last two decades have seen a tremendous increase in the use and development of quantitative models to estimate the distribution, abundance, and demographic rates of populations (e.g., Hobbs & Hooten, 2015;Hooten & Hefley, 2019;Kéry & Royle, 2016, 2021Royle et al., 2013;Royle & Dorazio, 2008;Schaub & Kéry, 2021;Williams et al., 2002), and it seems evident that the next two decades will see a tremendous increase in the application of these techniques to communities (e.g., Montaño-Centellas et al., 2021) or groups of multiple species.
Dynamic N-mixture and occupancy models allow for the formal modeling of demographic processes with counts or detection/ non-detection data of unmarked individuals (Dail & Madsen, 2011;MacKenzie et al., 2003;Royle & Kéry, 2007;Sollmann et al., 2015;Zipkin et al., 2014). However, a critical consideration in the effective use of N-mixture models is unmodeled heterogeneity in detection probability Duarte et al., 2018;Link et al., 2018) or ecological processes (Bellier et al., 2016). Unmodeled heterogeneity can lead to poor model fit, parameter estimates far from ecological truth, and consequently flawed inference and conservation action Duarte et al., 2018;Link et al., 2018).
Thus, effectively modeling patterns in detection probability and ecological processes is important for the utility of these model types (Bellier et al., 2016;Kéry & Royle, 2021).
Multispecies or community models (Dorazio & Royle, 2005;Gelfand et al., 2005Gelfand et al., , 2006MacKenzie et al., 2005) allow for sharing information about demographic and abundance parameters among species. Recent work on N-mixture models has examined the potential for estimation of visit-specific heterogeneity in detection probability (Dorazio et al., 2013;Martin et al., 2011), as well as hierarchical modeling of variation in mean detection probability among species (Gomez et al., 2018). Further, previous research has suggested approaches to use information from multiple species to inform the observation process (MacKenzie et al., 2005;Nichols et al., 2000). In this study, we build on this research to estimate speciesand visit-specific detection probabilities in a Bayesian hierarchical framework. Our central idea is simple: unknown and unmodeled factors that affect the detection probability of one species during surveys will often affect the detection probability of other similar species in the same way. This is true because components that affect detection such as weather, habitat, date, time, and observer are the same during each multi-species survey. In fact, we note that simultaneous observations of multiple species are not statistically independent observations. Thus, we use random effects to share information about heterogeneity in nuisance parameters (i.e., detection probability) among species and clearly demonstrate in a simulation study that this approach can lead to more accurate and precise estimates of model parameters when model assumptions are met.
To further demonstrate the efficacy of our approach beyond static N-mixture models, we use our shared parameterization for the detection process with a dynamic occupancy model (MacKenzie et al., 2003;Royle & Kéry, 2007) applied to detection/non-detection data collected on willow (Poecile montanus montanus, P. m. salicarius, and P. m. rheananus), marsh (P. palustris), crested (Lophophanes cristatus), coal (Periparus ater), blue (Cyanistes caeruleus), and great (Parus major) tits (Paridae) breeding across a P. montanus hybrid zone and elevational gradient in northern Switzerland (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020). These six species can be loosely sorted into low and high elevation guilds in Switzerland, where marsh, blue, and great tits typically peak in abundance near 500 m in deciduous forests, and willow, coal, and crested tits typically peak in abundance near 1500 m in coniferous forests (Knaus et al., 2018). Further, tits have often been used as model organisms to examine inter-and intra-specific competition (Alatalo et al., 1986;Dhondt, 1977) during the breeding (Dhondt, 2010;Minot & Perrins, 1986) and non-breeding (Dhondt & Eyckerman, 1980) seasons. Although the majority of studies have focused on competition within the previously described elevation guilds (e.g., Wittwer et al., 2015), previous research has demonstrated that in the absence of marsh tit, willow tits occupy lower elevation deciduous habitats . As global climates warm, montane vegetation communities shift upslope 4. While conceptually simple, our results have important implications for the future modeling of population abundance, colonization, persistence, and trajectories in community frameworks. We suggest potential extensions of our modeling in this paper and discuss how leveraging data from multiple species can improve model performance and sharpen ecological inference.

K E Y W O R D S
abundance, community, demography, detection, hierarchical model, multilevel model, Nmixture model, occupancy model, Paridae, species interactions (Beckage et al., 2008), albeit as a function of a variety of complex processes (Alexander et al., 2018;Scherrer et al., 2020). In response, avian communities often shift upslope as well (Freeman et al., 2018;Gasner et al., 2010;Paxton et al., 2016). Changes in the distribution and abundance of species can lead to novel competitive landscapes and communities (Alexander et al., 2015), as well as opportunities to examine the effects of interspecific competition on species distribution (Alexander et al., 2016;Legault et al., 2020). Thus, as marsh tits follow deciduous vegetation upslope, we might expect displacement of congeneric willow tits. In this study, we describe a simple multispecies approach to improve precision of estimates of abundance and occupancy. We then demonstrate the utility of this approach by examining interspecific competition in Swiss tits, where we demonstrate that range expansion of marsh tits may lead to range contractions in willow tits, potentially affecting hybridization rates and gene flow among distinct willow tit (P. montanus) subspecies.

| Simulation study: estimating abundance and shared detection
We simulated data to test the efficacy of sharing observation process information among species (s) to effectively model detection parameters under a static multi-species N-mixture model. We To simulate heterogeneity in detection probability among surveys, we first simulated a shared species-specific mean detection probability (p = 0.5). We simulated heterogeneity in detection ( i,j ) during each visit ( j) to each site (i ) that was also shared among species, where i,j are added on the logit scale and may represent unmodeled variation in observer skill and effort, local weather, date, time, and a myriad of other factors that affect simultaneously the detection of all species present. Crucially, this part of the model creates a covariation in the detection probability among the modeled species in a manner that we believe is frequently very realistic for animal surveys. The number of individuals of each species observed at each site during each visit (y s,i,j ) was simulated as a binomial trial given the number of individuals present (N s,i ) and detection probability (p s,i,j ), We analyzed the multi-species abundance data with two separate models. First, we used a classic approach to account for heterogeneity in detection (Martin et al., 2011), but treated the data for each species separately (i.e., without any shared parameters), Second, we jointly analyzed the data with shared heterogeneity in detection probability among species by simplifying the estimation of heterogeneity and retaining the structure for abundance from the previous approach, We simulated 100 datasets using R 4.0.3 (R Core Team, 2018), saved the medians, 95% credible intervals, and standard deviations of posterior distributions, and report mean squared error, mean signed difference and parameter calibration, or coverage (Little, 2006;Williams & Hooten, 2016) for each parameter estimated in the simulation.

| Case study: interspecific competition in Swiss tits
To further demonstrate the use of our shared-detection parameterization, we apply it to a community occupancy dataset to examine interspecific competition between the willow tit (P. montanus) and a congener, the marsh tit (P. palustris). We use detection/non-detection data for willow tit, marsh tit, crested tit, coal tit, blue tit, and great tit collected at 119 1 km 2 long-term monitoring sites (i ) below 1500 m in the Jura, Central Plateau, and pre-Alps ecoregions of Switzerland.
Each site was surveyed three times ( j) during each breeding season (mid-April to late June; 2004-2020; t) as part of the long-term Swiss Breeding Bird Survey (Monitoring Häufige Brutvögel; Schmid et al., 2004).
We modeled the observation process for each species as a function of a species-specific intercept and visit-specific random-effects (2) (3) To estimate heterogeneity in detection probability among species, sites, and surveys, we first modeled species-specific variation in mean detection probability ( s , Gomez et al., 2018) as a function of mean detection probability of all species (p) and random variation in detection probability among species ( ).
Thus, in this example, we can think of variation in speciesspecific mean detection probability ( ) as representing latent heterogeneity among species due to differences in visibility, behavior, or song distinctiveness and volume (Gomez et al., 2018), and in the case of occupancy models, abundance (Royle & Nichols, 2003).
We modeled the presence of each species at each 1 km 2 survey site (Z) as Bernoulli trials. During the first occasion (z s,i,1 ), we estimated presence as a function of initial occupancy probability ( s,i ).
During subsequent occasions, we estimated occupancy (z s,i,t ) as a function of the occupancy status at the previous sampling occasion (z s,i,t−1 ), and species-, site-, and time-specific persistence ( ) and colonization ( ) probabilities, We modeled the initial occupancy probability of each species at each site ( s,i ) as a function of a species-specific intercept ( 0,s ) and linear species-specific relationships with percent forest cover (f i ) and elevation (e i ), We used a similar approach to model persistence ( ) and colonization ( ) probability for crested, coal, blue, and great tits, Given the potential for interspecific competitive effects of marsh and willow tits, we included an effect of the presence of congeners at each site (z c,i,t ) in addition to the effects of percent forest cover and elevation, In other words, we included an intercept adjustment that accounted for the effect of marsh tit presence on willow tit persistence and colonization probability, and vice versa. We transformed covariates to standard normal deviates (i.e., forest and elevation covariates were z-standardized). We used vague priors for slope and intercept parameters for the models of initial presence, persistence, and colonization, : (0, 100).
We sampled using three MCMC chains for 500,000 iterations with an adaptive phase of 1000 iterations for each model for the simulated data. We discarded the first 250,000 iterations as a burn-in and retained every fifth saved iteration. For the tit data, we sampled using three MCMC chains for 1m iterations and discarded the first 500,000 as a burn-in. We report posterior distribution medians, 95% credible intervals, and the proportion of the posterior distribution on the same side of zero as the mean ( ) for parameters estimated in the case study. This value ( ) represents our confidence that a parameter is either greater or less than 0.

| Simulation study: estimating abundance and shared detection
Only 44 of the 100 single-species models of simulated data converged ( � R < 1.1) with the specified MCMC settings. All of the joint analyses converged. Parameter estimates from both the joint and separate analyses were centered around the true values used to generate the data ( Figure 1, Table 1). However, joint analyses substantially increased precision in estimates of abundance, mean detection probability, and detection heterogeneity among visits (Figure 1), where both mean signed difference and mean squared error were reduced compared with single-species models (Table 1).
persistence probability, and colonization probability of marsh, blue, and great tit were negatively affected by elevation (Table 2, Figure 3). The effects of forest coverage were mixed, where only blue tits consistently demonstrated negative effects of forest cover on demographic parameters (Figure 3). Initial presence probability was positively affected by forest cover for willow, crested, and coal tits and negatively affected by forest cover for blue tits (Table 2). Forest cover positively affected persistence probability for willow, crested, coal, and marsh tits and negatively affected persistence probability for blue and great tits (Table 1).
Colonization probability was positively affected by forest cover for willow, crested, coal, and great tits and negatively affected by forest cover for blue tits (Table 2).
In sites occupied by both marsh and willow tits, we found weak evidence for a decline in persistence probability of willow tits ( 3 = − 0.418, = 0.820; Figure 4), and no effect on persistence probability of marsh tits ( 3 = 0.242, = 0.734). In contrast, we found strong evidence that marsh tit occupancy was associated with reduced colonization by willow tits ( 3 = − 0.831, = 0.980 ; Figure 4) and willow tit occupancy with reduced colonization probability by marsh tits ( 3 = − 1.292, = 0.999).

| DISCUSS ION
In many statistical community models, researchers focus on understanding interspecific effects on occupancy, abundance, and demography. In this paper, we demonstrate the importance and utility of sharing information among species to inform the observation process, thereby strengthening our inferences about the demographic processes. First, we use this approach with simulated count data where we perfectly meet model assumptions. Our results demonstrate that sharing information about heterogeneity in detection probability among species can be an effective tool for improving the precision and accuracy of detection estimates. We then apply this approach to more complex, dynamic occupancy data on tits breeding at lower elevations (< 1500 m) in northern Switzerland.
We use data from all six tit species breeding in northern Switzerland to inform detection probability and habitat relationships, while constraining inference regarding interspecific competition to the two most closely related species, marsh and willow tit.
In northern Switzerland, marsh tits occupy lower elevation sites, while three subspecies of willow tit occupy higher elevation habitats (Knaus et al., 2018). Previous research in island systems has F I G U R E 1 Boxplots of medians (top) and standard deviations (bottom) of posterior distributions for estimates of species-specific mean abundance on the log-scale ( ; left), mean detection on the logit scale ( ; middle), and heterogeneity in detection probability among visits ( ; right) for the simulated count data. Blue boxes and points are results from models with shared heterogeneity in detection among species, while red boxes and points are from models that analyze each species data separately. Note that joint estimates are shared among species for and Species TA B L E 1 Mean signed difference (MSD) and mean squared error (MSE) for the medians of posterior distributions given the true values used to simulate the data, means of the standard deviations of the posterior distributions ( ) from the simulations, and posterior distribution coverage for N-mixture models with shared detection heterogeneity parameters (Joint) and models that treat each species individually (Separate) Note: is the initial abundance for each simulated species on the log-scale, is the intercept for detection probability of each species on the logit scale, and is the amount of detection heterogeneity on the logit scale.
TA B L E 2 Medians ( ), 95% credible intervals, and proportion of the posterior distribution on the same side of zero as the median ( ) for posterior distributions of the effects of percent forest cover ( 1 ) and elevation ( 2 ) on initial presence probability ( ), persistence probability ( ), and colonization probability ( ) of six species of tit breeding in northern Switzerland (2004-2020)  . In response to climate warming, montane vegetation (Beckage et al., 2008) and avian communities (Freeman et al., 2018;Gasner et al., 2010;Paxton et al., 2016) move upslope. Our results demonstrate that the presence of either species reduces the probability of colonization by the other (Table 2), and we found weak evidence ( 3 = − 0.430, = 0.825) that the presence of marsh tits reduces the persistence probability of willow tit populations ( Figure 4). However, over time, even a weak tendency toward displacement of willow tit by marsh tit may lead to range reduction and reduced interbreeding among willow tit subspecies (Knaus et al., 2018).
The primary tenet of the approach we describe is that detection probabilities of all surveyed species positively covary to some degree, and ideally strongly so. We expect that such positive covariation is common when species are sampled using the same protocol and do not have diametric behavior. For instance, the tit populations in this study are most often detected acoustically and singing activity commonly depends on weather, date, time of day, or habitat. If species are included that differ radically in behavior (e.g., tits and owls), we would not expect positive covariation of detection parameters.
Readers should note that major violations of model assumptions will lead to incorrect parameter estimates in all hierarchical models, and these issues can be particularly problematic when using the various classes of N-mixture models ( Elevation (m) Link et al., 2018). Thus, we suggest that more complex and flexible random-effects structures may be employed in the future as well.
For instance, researchers might employ covariance approaches, estimating correlations among detection probabilities as a function of shared traits (see Nichols et al., 2000) or even phylogeny (e.g., Abadi et al., 2014). Researchers might also use informative priors to inform covariance in processes among species, or explicitly estimate correlations (Riecke et al., 2019). The approach described herein also allows for the simple inclusion of survey-specific covariates, such as survey date, local weather covariates, or nesting visit-specific random effects within observer-specific random effects. While fully shared observation parameters may seem to be too generalized, we emphasize that this has been a hidden assumption of every existing community modeling framework that does not explicitly account for the observation process. Finally, given the strong linkage between abundance and detection of species presence (Royle & Nichols, 2003), we urge caution when using these model types with detection/non-detection data for groups of species with extreme variation in abundance.
The estimation of abundance, distribution, and patterns in ecological communities is of central importance to a wide number of ecological fields, and continued advances in the estimation of these ecological processes will be critical to the progression of ecological knowledge (Hooten & Hefley, 2019;Hooten et al., 2017;Kéry & Royle, 2016  2018; Link et al., 2018), these models can estimate abundance and demographic parameters when detection heterogeneity is adequately modeled (Ficetola et al., 2018;Kéry, 2018;Royle, 2004;Veech et al., 2016). Numerous studies have recently demonstrated methods to integrate other data types with N-mixture models to improve performance. For example, telemetry (Ketz et al., 2018;Popescu et al., 2017;Schmidt et al., 2015), band-recovery (Zhao et al., 2019), and genetic (Furnas et al., 2018) data have all been used to help inform demographic parameters of single species in integrated modeling frameworks. In this paper, we use the statistical dependence of community-level count and detection/nondetection data to inform heterogeneity in detection probability.
Further, we suggest that linking the observation process among species may even allow for the inclusion of other data types (e.g., capture-recapture) from a single species to inform the observation of communities. While we urge continued caution and further simulation-based validation of new techniques, we suggest that the thoughtful sharing of information among species within community-level models (Gomez et al., 2018;MacKenzie et al., 2005;Nichols et al., 2000, this paper) has strong potential to improve ecological inference moving forward.

ACK N OWLED G M ENTS
We thank the hundreds of volunteers and employees at the Swiss Ornithological Institute for development and maintenance of the Swiss Breeding Bird Survey dataset (Monitoring Häufige Brutvögel; Schmid et al., 2004). We thank the associate editor and two anonymous reviewers for thoughtful and constructively critical contributions to the manuscript.

CO N FLI C T O F I NTE R E S T
There are no conflict of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
R script for simulating and analyzing the data and the data used in this manuscript are archived at the Dryad Digital Repository (https:// doi.org/10.5061/dryad.573n5 tb8s).