Integrated distance sampling models for simple point counts

Point counts (PCs) are widely used in biodiversity surveys, but despite numerous advantages, simple PCs suffer from several problems: detectability, and therefore abundance, is unknown; systematic spatiotemporal variation in detectability produces biased inferences, and unknown survey area prevents formal density estimation and scaling-up to the landscape level. We introduce integrated distance sampling (IDS) models that combine distance sampling (DS) with simple PC or detection/nondetection (DND) data and capitalize on the strengths and mitigate the weaknesses of each data type. Key to IDS models is the view of simple PC and DND data as aggregations of latent DS surveys that observe the same underlying density process. This enables estimation of separate detection functions, along with distinct covariate effects, for all data types. Additional information from repeat or time-removal surveys, or variable survey duration, enables separate estimation of the availability and perceptibility components of detectability. IDS models reconcile spatial and temporal mismatches among data sets and solve the above-mentioned problems of simple PC and DND data. To fit IDS models, we provide JAGS code and the new IDS() function in the R package unmarked . Extant citizen-science data generally lack adjustments for detection biases, but IDS models address this shortcoming, thus greatly extending the utility and reach of these data. In addition, they enable formal density estimation in hybrid designs, which efficiently combine distance sampling with distance-free, point-based PC or DND surveys. We believe that IDS models have considerable scope in ecology, management, and monitoring.


Introduction
Point count methods are among the most widely used and longest-standing protocols in wildlife surveys worldwide (Rosenstock et al. 2002;Darras et al. 2021).Simple point counts (PC) are short surveys in which observers count all individuals of some set of species (single species to entire communities) detected without distance constraints (i.e., unlimited distance), or within a predefined distance from the observer.Point count methods are logistically uncomplicated and thus have been adopted in the design of countless surveys worldwide, e.g., the North American Breeding Bird Survey/BBS (Sauer et al. 2017), and many European national BBS or bird atlas schemes (Balmer et al. 2013).In addition, a wide range of what in essence are point count methods, although varying from highly standardized to essentially design-free, is at the core of rapidly growing citizen-science projects such as eBird (Sullivan et al. 2009).
Despite the prevalence of simple point counts, their simplicity is not without drawbacks; e.g., it is not possible to estimate true abundance or occupancy if visits to points are unreplicated (but see Lele et al. 2012 andSólymos et al. 2012).In addition, PC data are non-spatial in the sense that the area from which the detected animals are drawn is unknown.This prevents spatial extrapolation for rigorous estimation of regional population sizes.Similarly, integrated analyses of data from different schemes is hampered due to commonly occurring spatial mismatches.Finally, variable survey duration is commonplace and creates temporal mismatches in the data (Pacifici et al. 2019).This further complicates joint analyses from multiple survey schemes that use point count methods.
In planned surveys, additional information is often collected during PC surveys that permits estimation of detection probability (Nichols et al. 2009).Such extra information includes replicated counts (Royle 2004), double-observer surveys (Nichols et al. 2000), removal counts (Wyatt 2002;Dorazio et al. 2005), distance information (Marquez et al. 2007;Buckland et al. 2015), or locational information from recognizable individuals which enables spatial capture-recapture (SCR) models to be fit (Borchers and Efford 2008;Royle et al. 2014).These survey protocols permit estimation of abundance, and thus assessment of status and trends free from any bias produced by imperfect detection and by unmodelled temporal or spatial patterns in detectability (Kéry & Royle 2016, process permits estimation of separate detection functions, along with different parameters linking these functions to covariates, for all three data types.Estimation of separate detection functions, when needed, can accommodate any systematic differences between survey schemes.Combining PC and DND data with the more information-rich DS data not only enables estimation of detection probability, but also makes the resulting abundance estimates area-explicit: effective survey areas for PC and DND surveys become estimable, and population density is estimated with improved precision.Thus, IDS models can reconcile all discrepancies, including spatial and temporal mismatches, among these widespread data types. In this article, we begin by formally describing IDS models.Then, we use simulation to first demonstrate identifiability of our model when separate detection functions are estimated for each data type, including separate parameters for detection function covariates.Next, we explore the effects of adding variable amounts of the more information-rich but more "expensive" DS surveys to a larger sample of the less information-rich but "cheaper" PC data.Following that, we demonstrate the identifiability of a model that includes availability in addition to perceptibility, provided that survey duration is variable; variable survey duration is one of the types of extra information which enable availability in addition to perceptibility to be estimated (e.g., Solymos et al. 2013).Finally, we showcase IDS models with real data, using the Oregon 2020 Bird Survey Project (Robinson et al. 2020) as a case study.As part of the case study, we demonstrate the ability of IDS models to allow for different levels of heterogeneity in the detection functions estimated for different portions of the data (Oedekoven et al. 2015).Such accommodation of survey-specific features of the observation process will be particularly important when reconciling data from heterogeneous survey protocols in a single integrated model.
We have implemented a range of IDS models in the new fitting function IDS() in the R package unmarked (Fiske & Chandler 2011), to permit model fitting by maximum likelihood, and we provide BUGS code for Bayesian inference using JAGS (Plummer 2003).We believe that IDS models have a large scope of application for exploiting PC and DND data in a more rigorous and synthetic manner, and to obtain less biased and larger-scale inferences about abundance and density, particularly for large amounts of citizen-science data.

Integrated Distance Sampling (IDS) models
We develop joint likelihoods, i.e., integrated models (Besbeas et al. 2002;Miller et al. 2019;Kéry & Royle 2021: chapter 10;Schaub and Kéry 2022), for the following data types, which we assume to observe the same density and availability processes.We note that i indexes different sites across data types.Also, our current models assume closure and the absence of any temporal replicates at a site, though both assumptions may be relaxed in future developments.For inference about density, first, for the DS data we adopt a hierarchical distance sampling model (Royle et al. 2004)   This is a key advantage of DS methods, which associates i N with a well-defined area.It makes abundance estimates in a DS protocol area-explicit, in contrast to most abundance estimation protocols other than SCR (Borchers & Efford 2008;Royle et al. 2014).For songbirds, the availability probability  will be mainly a function of singing rates (Solymos et al. 2013), which cannot be estimated from distance data alone.Hence, conventional distance sampling (CDS) requires either the assumption of perfect detection at distance 0 or else acceptance that inferences will be restricted to the available part of a population only (Buckland et al. 2015).However, availability becomes estimable in a DS model if certain extra information is collected, e.g., from multiple observers (Borchers et al. 1998), replicated surveys (Chandler et al. 2011), time-removal (Farnsworth et al. 2002), or alternatively from variable survey duration in an IDS model, as we will show.
Second, for the PC data we adopt a variant of the binomial N-mixture model (Royle 2004).
Simple PCs are neither area-explicit, nor can detection probability be estimated without temporal replication (though see Lele et al. 2012 andSólymos et al. 2012) To summarize our IDS models, for regular DS data we specify likelihood ds L (Royle et al. 2004), for PC data pc L (Royle 2004), and for DND data dnd L (Royle & Nichols 2003).Importantly, for both PC and DND data, we assume a latent DS observation process protocol and estimate p by integration of a detection function with parameters that become estimable in an IDS model.Under independence assumptions, i.e., when at most a negligible portion of sites appears in more than one data set, we define the following joint likelihoods for three variants of an IDS model: (which we call model IDS1) and

Tests and demonstrations of integrated distance sampling models with simulated and real data Simulation 1: Identifiability of separate observation process parameters in IDS1 and IDS2
To demonstrate identifiability of the new models, we analyzed simulated data sets and estimated parameters for separate detection functions in an IDS model with either DS + PC data or DS + DND data, i.e., in the IDS1 and IDS2 case.We used function simHDS() in the R package AHMbook to simulate two data sets with DS data from 250 sites, and PC or DND data from another 1000 sites.To obtain PC data, we first generated DS data, and then discarded all distance information, just retaining one count per site each, and to produce DND data we additionally quantized the resulting counts.Mean density was kept constant at 1, following our assumption of a shared density process.The scale parameter  in the half-normal detection function was set at 100 m for the DS data and was varied randomly between 10 and 130 m for the PC and DND data sets.
Thus, the key criterion for identifiability of the new models was how well estimates of  matched their true values in the data simulation.In the submodel for the DS data sets, we chose a truncation distance of 200 m.In this simulation we aimed to establish the identifiability of the new models in their simplest form only.That is, we implicitly assumed availability to be 1 and did not use any covariates in either density or detection.We used JAGS (Plummer 2003) to fit IDS1 or IDS2 to 1000 data sets each.
In Simulation 1B (Appendix S2: Section S1) we continued our investigation of parameter identifiability and of estimator performance in model IDS1.We varied all of the following four settings independently according to a response-surface design: average density, detection function scale for both DS and PC data, and DS truncation distance.We again used JAGS to fit all models.

Simulation 2: Identifiability with distinct covariate effects in the observation model
We conducted two sets of simulations to answer the following related questions: (i) Does the IDS model allow DS and PC detection to have different covariate relationships?(ii) Are relationships still identifiable if the same covariates are related to both detection and density?We answered these questions by simulating data sets with DS and PC data from 200 and 1000 sites, respectively.
Density was governed by an intercept of 1 on the natural scale and an effect of 1 of one covariate ("habitat").The half-normal detection function  had an intercept of 100 and 150 m on the natural scale for DS and PC data, respectively.In a first analysis we used simHDS() to simulate 1000 data sets with these specifications, and where the half-normal detection function  , on the log-scale, was affected by another covariate "wind" by independently drawing two random numbers from a Uniform(-0.5, 0.5) distribution.In the second analysis, we used a modified version of function simHDS() to simulate another 1000 data sets with the same specifications as above, except that now we generated log-scale effects of the same covariate as for density (i.e., "habitat") by independently drawing two U(-0.5, 0.5) random numbers for the DS and PC data sets as their coefficients.We used the new IDS() function in R package unmarked to fit the data-generating model.We discarded numerical failures, which we conservatively identified by standard errors that were either NA or >5, or by MLEs that were >10 times their true values.

Simulation 3: How many DS sites are required to obtain adequate estimates of density?
We simulated 1000 data sets with PC data from 200 sites, to which we added DS data from 1-100 sites in six mixing ratios.Density was governed by an intercept of 1 on the natural scale, with one habitat covariate with coefficient 1.The detection function  was 70 m in the PC and 100 m in the DS data, and we chose a truncation distance of 200 m in the latter.We generated a total of 6000 data sets (1000 for each level of the mixing ratio factor) and fit the IDS1 model using function IDS(), discarding numerical failures as in Simulation 2.

Simulation 4: How well can availability be estimated in an IDS model?
We simulated 1000 data sets that resembled our case study below: each had DS data from 3000 sites, and PC data from either 1000, 3000, or 6000 sites.DS survey duration was kept constant at 5 min, but was varied between 3 and 30 min in PC surveys, with a strong right skew, as found in the case study data.Density was goverened by an intercept of 1 on the natural scale and with a habitat covariate with coefficient 1, detection function  was 70 m in the PC and 100 m in the DS data, with a truncation distance of 200 m in the latter.Singing rates varied between 0.1 and 2, corresponding to a probability of 0.1-0.86 to sing at least once over a 5 min interval.We fit IDS1 using the IDS() function and discarded numerical failures as in Simulation 2.

Case study: American Robins in the Oregon 2020 Bird Survey Project
We used the IDS3 model to estimate population density of American Robin (Turdus migratorius) in Benton and Polk counties, Oregon.The 3680 km 2 area in Western Oregon is bounded on the east by the Willamette River and its floodplain, while the western portions include the Coast Range mountains.Silviculture of coniferous forest is the dominant land use in the mountains.Nearly every square kilometer contains a narrow, lightly travelled road for timber harvest, which allowed access for bird surveys.The eastern floodplain sections contain a mix of agricultural uses, mostly festucoid grass seed fields and orchards, and suburban development.
DS surveys were conducted every 0.8 km along accessible roads throughout the study area, and every 200-m in an off-road grid placed over the William L. Finley National Wildlife Refuge, producing a total of 2,912 sites sampled and 2020 American Robins detected (Robinson et al. 2020).DS surveys were conducted during the breeding season (April 30-July 11) from 2011 to 2013 by WDR.Each survey followed the Oregon 2020 protocol (Robinson et al. 2020), which used 5-minute stationary counts initiated between 30 min before sunrise and noon on days with no or little rain.All birds detected by sight or sound were recorded with an estimated distance from the observer (verified with a range finder when possible) following standard distance sampling protocols (Buckland et al. 2015).
We combined DS surveys with opportunistically-gathered citizen-science PC data from the eBird database (Sullivan et al. 2009), using checklists from 2011-2017 in Benton and Polk counties.
After stringent filtering (see Appendix S2: Section S2) and geographic sampling, 1060 PCs with 819 detections of American Robins were included.We filtered data to include only complete checklists using stationary protocols and personal locations, conducted during the breeding season.
We further filtered data to include only checklists with durations between 3-30 minutes that were conducted between sunrise and seven hours after sunrise.Finally, we applied geographic sampling to reduce the effects of popularly birded sites by overlaying a 200m grid over the study area and randomly selecting only a single checklist from each grid cell.To illustrate the combination of DS data with both simple PC and DND data, we randomly subset half of the PC data and reduced it to DND.This resulted in 521 simple PCs with a total of 394 individuals counted, plus 538 DNDs, among them 228 with at least one detection of the species.See Appendix S2: Section S2 for eBird query details.
For DS data, we selected a truncation distance b ds of 300 meters.We binned the distance data into 50 m distance classes.For the analysis of PC and DND data, an upper distance limit b pc and b dnd of 500 meters was adopted, assuming that observers do not detect individuals further away than 0.5 km (the 0.99 quantile in the Oregon 2020 database (Robinson et al. 2020) was 0.4 km).For all data types, we assumed identical parameters for annual density and availability.We modelled density  with a random intercept for year, and with quadratic terms for elevation and percentage of canopy cover in a 315 m radius around the observer location.This radius was selected as it was previously found to be the most predictive of abundance for this species of the radii considered (Hallman and Robinson 2020).For availability, we adopted the model of Solymos et al. ( 2013) linking availability probability with activity rate ϕ, and used quadratic terms for day of the year and minutes since dawn on the log activity rate.
The observation process in the designed DS surveys in the Oregon 2020 project may differ from point count surveys recorded in eBird, even after stringent filtering.Therefore, we allowed for different detection functions for the DS and the PC/DND portions in our analysis by fitting separate intercepts in the half-normal detection scale  .Moreover, to accommodate possibly different levels of detection heterogeneity among sites, we specified site-specific random effects in  with a potentially different variance for the DS and PC/DND portions of the data (Oedekoven et al. 2015).
In addition, we modeled  using the percentage of urban area and percentage of canopy cover, both in a 165 m radius around the observer location; these slope parameters were shared between DS and PC/DND data.We used a smaller radius for canopy cover in the detection function as the distance that an observer can detect is impacted more heavily by nearby environmental conditions.
Elevation was obtained from the Oregon Spatial Data Library (Oregon Spatial Data Library 2017), urban land cover was obtained from the U.S. Geological Survey's National Gap Analysis Project (USGS 2011), and canopy cover was obtained from Landscape Ecology, Modeling, Mapping and Analysis's gradient nearest neighbor structure maps (LEMMA 2014).
We processed data in R (R Core Team 2019) and fitted the model in JAGS, using the R package jagsUI (Kellner 2016).For all parameters, we chose vague priors (see BUGS model on ZENODO for specifics).We assessed the goodness of fit for our model for each data portion separately using posterior predictive checks (Conn et al. 2018) with a Freeman-Tukey discrepancy measure between the observed and the expected counts for the DS and PC data (Kéry & Royle 2016).For the binary DND we used the observed number of sites with detections as a test statistic and checked whether the actual value for the DND data was within a 95% credible interval (CRI) of this statistic in data sets replicated under the model.This approach suggested an adequate fit of the model overall: Bayesian p-values for the DS part of the model revealed slight underdispersion of the data, while the PC and DND parts of the data indicated a fitting model (Appendix S2: Table S6).
We obtained posterior predictive distributions of abundance to project density onto a raster map with elevation and canopy cover for each of the 3874 1-km 2 grid cells in Benton and Polk counties to obtain an abundance-based species distribution map of the Robin in our study area.We also fit a simpler variant of the model using the IDS() function in unmarked (Fiske & Chandler 2011) to illustrate both the Bayesian and the maximum likelihood approach.Code and data to replicate the case study can be found on Zenodo.

Simulation 1: Identifiability of separate observation process parameters in IDS1 and IDS2
In an IDS model, separate detection functions can clearly be estimated under both IDS1 (combining DS and PC data) and IDS2 (combining DS and DND data); see Fig. 1 and Appendix S2: Table S1.
There were no signs of bias in either model: % relative bias was <<1% for all sigma's and <2% for the abundance estimates at sites with N>0.Credible interval (CRI) coverage was close to the nominal level of 95% for all estimators.Not surprisingly, precision was slightly lower in model IDS2 than in IDS1 (see middle of Fig. 1).In addition, simulation 1B confirmed the good frequentist operating characteristics of the estimators in IDS models under an even wider range of conditions (Appendix S2: Section S1, Appendix S2: Table S2).

Simulation 2: Identifiability with distinct covariate effects in the observation model
In the first set of simulations, where two distinct covariates affected density and the detection function, but where the effects on the latter were distinct for the DS and PC portions of the data, we discarded 23 sets of estimates as invalid.The remaining 977 sets of estimates indicated that this model was identifiable and showed little or no signs of bias (Fig. 2, left; Appendix S2: Table S3 left).In the second set of simulations, where a single covariate independently affected density and the two detection functions, we discarded 66 invalid sets of estimates.The remainder again showed this model to be identifiable (Fig. 2, right; Appendix S2: Table S3 right).

Simulation 3: How many DS sites are required to obtain adequate estimates of density?
The IDS model showed excellent performance with as few as 20 DS sites (Fig. 3), with relative bias <1% for all estimators and CI coverage at or near nominal levels (Appendix S2: Table S4).
However, the number of numerical failures increased greatly when decreasing numbers of DS data were added in the integrated model; from only 2 out of 1000 when 100 DS sites added, to 85 with 20 DS sites, and to 490 out of 1000 when 1 DS site was added.

Simulation 4: How well can availability be estimated in an IDS model?
Sampling distributions of density estimators were all concentrated around the true value.There were long right tails, but these became more symmetrical with larger sample sizes.Singing rate ( ) estimators were precise up to values of about 0.8, 1.3 and 1.4, respectively, for 1000, 3000 and 6000 PC sites, but became very imprecise for greater values of the singing rate.Presumably, this was because overall availability reached an asymptote close to 1 when singing rates were very high, making precise estimation of  hard (Fig. 4).Overall, there was a slight positive bias in both density and singing rates, but it declined from 14 and 10% with 1000 PC sites to 3000 and 2% with 6000 sites, while CI coverage was always at nominal levels (Appendix S2: Table S5).Relative bias of the detection function  for both data types was 0.76% or less throughout.

Case study: American Robins in Oregon
Over all surveys considered, mean survey date was June 7, and time since dawn ranged 17-519 min (mean 229).At mean date and time since dawn, availability within a one minute survey was estimated at 0.27 (95% CRI 0.15-0.46;Appendix S2: Table S6).Estimated availability peaked soon after dawn, decreased during the next five hours, then increased again, and tended to increase slightly throughout the season.Density was estimated to be highest on plots with a canopy cover of ~40% and to decrease with elevation.Median estimates for density varied between 2 and 47 individuals per km 2 .Maxima were found in the foothills where open woodlands transition from the floodplain agricultural zones into the denser forests at higher elevations, while minima were found in the most intensively harvested woodlands (Fig. 5).Over the entire study area, we estimated a population size of 74,344 American Robins (95% CRI 48,039-117,944).Interestingly, on average the estimated detection function was not different between the DS and PC/DND portions of the data (see parameter 'mean.sigma' in Appendix S2: Table S6).However, there was greater heterogeneity in the detection function  for surveys on eBird than for regular DS surveys conducted within the Oregon 2020 project (see parameter 'sd.eps').

Discussion
We discovered how simple point count (PC) or detection/nondetection (DND) data can be formally integrated in a model together with distance sampling (DS) data, to estimate separate parameters of an underlying latent DS observation process.This allows estimation of full detection probability parameters for all three data types.Moreover, integrating DS data makes abundance estimates from PC and DND data area-explicit.Thereby, IDS models achieve a formal spatial calibration of PC and point-indexed DND data, as well as a reconciliation of spatial mismatches between all three data types.IDS models thus solve two major problems that plague simple point count surveys producing PC or DND data: that detection probability and effective survey areas are both unknown.The main assumption of our IDS model is a shared density process: that either density is identical among all sample locations, or that differences can be explained by covariate regressions with coefficients that are identical for all data types in the integrated model.These assumptions should be reasonable when all data types are collected in the same general area, and they may also hold when data sets are from disjoint regions.However, as in perhaps all cases where different data sets are combined in a single analysis, this is a judgment call on the part of the analyst.
We believe IDS models have a large scope of application and can facilitate use of the large amounts of currently available PC data, such as the North American BBS (Sauer et al. 2017) in more formal analyses of abundance that account for imperfect detection They may also be applied for carefully quality-controlled eBird data (Sullivan et al. 2009), as illustrated in our case study.We have shown that only a relatively small amount of regular DS data are required to supplement simple PC data when used in an IDS model.Our findings agree with related work with other types of integrated models that demonstrate the benefits of combining even small amounts of data with a higher information content with less informative, but cheaper data (Dorazio 2014, Zipkin et al. 2017, Doser et al. 2021).This suggests that the scope of inference of point count surveys may be substantially increased by adding even a relatively small number of sites where the additional distance information is collected.
In our case study we found that perceptibility was not different on average between the DS data contributed by the Oregon 2020 project and the PC data obtained on eBird: the intercepts of the detection function  were hardly different between these two portions of the data used in the IDS model.However, allowing for random variation of the detection function  between surveys (i.e., sites) and for possibly different magnitudes of that variation between the two types of data revealed greater heterogeneity among surveys in the eBird data base than among surveys in the Oregon 2020 project.Arguably, this makes intuitive sense, since it appears likely that consistency between surveys was higher in the latter than in eBird.In addition, our case study emphasizes how careful modeling of patterns in the detection function of an IDS model can help to make data from different protocols more 'comparable', or 'alike', by explicitly allowing for their differences in terms of the observation processes that produce them.This is a great strength of IDS models and of parametric statistical inference in general.
Many survey data typically have large variation in duration (Solymos et al. 2013) and thus there is also a need for temporal mismatch among datasets to be addressed.We conceive of this as an availability process (Kendall et al. 1997), where over time some activity such as singing puts individuals at an increasing risk of being detected.Hence, survey duration is naturally informative about availability.However, this part of our model presents more challenges.Population closure is required and hence surveys should probably not be very long in duration.More importantly, this part of our IDS models has the form of a single-visit occupancy or N-mixture model (Lele et al. 2012): estimability hinges upon a continuous, "private" covariate that affects detection, and in our case, availability.Such models are identifiable (Dorazio 2012), but they rely strongly on parametric assumptions and may lack robustness to violations of those assumptions (Knape & Korner-Nievergelt 2015).Our simulations and the case study both showed availability to be identifiable in an IDS model.However, our study species was chosen to be fairly common.In contrast, in rarer species, there may well be challenges when attempting to estimate that parameter.When information to estimate availability is too sparse, estimates may tend towards the boundary of 1; however, this is what any model ignoring availability assumes implicitly anyway.
Therefore, IDS models that estimate availability must be developed and applied with care.
Any extra information about availability should be incorporated in the model, such as data from multiple observers (as in mark-recapture distance sampling; Borchers et al. 1998), replicated surveys (Chandler et al. 2011), time-of-detection and time-removal data (Farnsworth et al. 2005;Alldredge et al. 2007;Solymos et al. 2013, Amundson et al. 2014).Alternatively, availability parameters may be estimated from altogether different data types, such as recordings of individual singing behavior, or perhaps even taken from the literature.We note that Solymos et al. ( 2013) had good success with the integration of time-removal and distance sampling data, but in a simpler model that did not involve estimation of a detection function for the time-removal data.
We can envision at least four major extensions to the IDS models described in this paper.
First is the accommodation of survey sites included in the dataset which were sampled using multiple protocols.This induces a dependence that must be addressed in the construction of the joint likelihood.Second, IDS models could be developed for other survey geometries, such as line transects or search-encounter designs (Royle et al. 2014, Mizel et al. 2018).Third, allowing for open populations and demographic processes (Kéry & Royle, 2021: chapters 1 and 2) will be an important extension that may open up avenues for truly large-scale demography models (see also Appendix S1).Fourth, additional data types may be incorporated in the model, such as opportunistic data conceptualized as point patterns (Farr et al. 2021), time-to-detection data (Strebel et al. 2021), aggregated counts (Schmidt et al. 2022), and data from autonomous recording units (ARUs) (Doser et al. 2021).For instance, IDS models may be beneficial for ARU data by allowing estimation of the "listening range" of these devices under widely varying conditions, while additionally exploiting the information on singing rate contributed by the ARU data.
In summary, we believe that IDS models have a substantial scope to improve analyses of widely available simple PC and DND data obtained in citizen-science schemes, as well as the increasing amount of ARU data.IDS models may serve as a keystone of the formal, model-based unification of various data types both from designed and less-designed to even design-free surveys, to great mutual benefit.We find it fascinating to see how DS and simple PC or DND data both contribute two essential pieces of information towards the full IDS model: DS data contain most information about the detection function, while the heterogeneity in survey duration commonly found in simple PC/DND data enables estimation of the availability process.This neatly illustrates the fact that the future of biodiversity monitoring arguably lies in a combination of both designed surveys and carefully chosen citizen-science schemes.Sample size in both simulations is 1000 data sets.See also Appendix S2: Figure S1 and Tables S1-S2.S3.S4.

Figure captions
Fig. 4: Sampling distributions of estimators of density (lambda,  ) and of activity/singing rate (phi,  ) in an IDS model with availability (  ) fit to data from 3000 DS sites, plus 1000, 3000, or 6000 PC sites added (Simulation 4, n = 930, 866, and 997 valid analyses).Red denotes truth or absence of estimation error, dashed blue shows mean of valid estimates.See also Appendix S2: Table S5.
Section S1: Overview of the three types of data As currently described in the article and implemented in the function IDS() in the R package unmarked, three closed-population data types can be jointly analyzed: distance sampling (DS) data, an obligate ingredient of every IDS model, and then one or both of two further kinds of distance-free data: simple point count (PC) data or detection/nondetection (DND) data.All three are assumed to be collected by a stationary observer, i.e., during point-based surveys.However, the extension to data collected along line transects, i.e., by a moving observer, is conceptually straightforward.
Figure S1 represents such point-based survey data in a hypothetical landscape inhabited by 100 birds (red circles) preferring yellow over blue.Nine stationary observers (black triangles) survey birds out to some maximum distance (black circles) and collect either distance sampling data (left), distance-free point counts (middle), or distance-free detection/nondetection data (right).To avoid clutter in these maps, we here assume perfect detection, but in reality some individuals would typically be missed in distance sampling and point count surveys, and likewise for presences in detection/nondetection surveys.We also note that we depict all three survey types as being conducted from the identical nine locations.In practice, however, different locations (and also maximum survey distances) will typically be chosen, either by design or opportunistically (see also Figure S2).

Section S2: Overview of IDS models
An IDS model is a probabilistic representation of the processes that we imagine underlie these three data types: (1) the completely shared state process, (2) a partially shared observation process, and (3) what can be called an aggregation/summary process; the latter is completely different for each data type.
(1) Shared latent state process: a point process summarized by density per unit area Figure S2 shows the shared state process, red points are again bird locations.Triangles are the locations of a human observer doing "point-transect" distance sampling surveys (left), distance-free point counts (middle), or distance-free detection/nondetection surveys (right).Survey locations are here shown to be different among survey types, as they will most often be in practice.Circles indicate the maximum distance out to which birds are surveyed, and this truncation distance is here depicted to be different among survey types, as will also typically be the case in the real world.The parameters linking density to the environmental covariates are assumed to be identical for all data types combined in an IDS analysis.Therefore, we call this a completely shared latent state process, both structurally (i.e., in terms of the density model) as well as in terms of the actual parameter values linking density with the landscape.(2) Partially shared observation process: a distance sampling detection model For all data types, the fundamental observation process in an IDS model is a distance sampling detection function (Figure S3).This is a novel and key assumption of our models and a distinguishing feature of our models from other types of integrated models that contain distance sampling as one of their ingredients, such as Solymos et al.In the multinomial observation model for the binned distance sampling data, the cell probabilities in  are given by an integration of the detection function over the chosen distance breaks that define the distance bins.In the binomial/Bernoulli observation model for the point counts and the detection/nondetection data, ̅ is obtained as the integral of the detection function from 0 out to either a truncation distance (if there is one), or else to a sufficiently large distance that is chosen so that no individual can be detected further than that.

Section S3: Further comments and outlook on future model extensions (1) Shared state process
In our IDS models the state process is assumed to be identical among all data sets combined in a single analysis.That means that the parameters relating density and the environmental covariates must be identical for all data sets.In practice this will mean that all surveys are conducted in the same region, as depicted in Figure S2.Sometimes, however, we may want to combine data types in an IDS model that do not show an adequate geographical interspersion, but rather may come from separate regions.Whether the integration of these data makes sense or not is a judgment call exactly as in any other study where we combine in a single analysis data from multiple places or from different times. (

2) Observation process
The key idea of an IDS model is to treat distance sampling as the fundamental observation process for all three data types.That is, we posit a distance-related, monotone decline of detection probability with increasing observer-bird distance and we describe it with a detection function such as the half-normal (though other detection functions could be chosen instead; Buckland et al. 2015; Chapter 8 in Kéry & Royle, 2016).The parameters of this detection function can be directly estimated from distance sampling data alone, but not from point count or detection/nondetection data alone.However, as part of an IDS model, i.e., in combination with distance sampling data and the assumption of a shared state process, the parameters of this ("latent") detection function become estimable also for point count and detection/nondetection data, as we demonstrate in Simulations 1-3 in our paper.
(3) Availability process This conceptual outline describes how a basic IDS model "work".That is, an IDS model that is based conceptually on what is called conventional distance sampling (CDS; Buckland et al. 2015), where either availability probability is equal to 1 (e.g., every individual present in the surveyed area around each point is likely to sing at least once during a survey) or else where all data types combined in an IDS model stem from surveys with identical duration; in this case availability can be assumed to be identical for all.In many surveys of songbirds, the availability process will be governed mainly by the singing rate of an individual.Sometimes, we will have extra-information which enables us to estimate the parameters of this availability process in addition to those of the other processes described above; we show this in Simulation 5 of the paper.Such extra information may come in the form of varying survey duration, removal counts (Solymos et al. 2013), or external estimates of singing rates from outside of a study.
(4) Temporal replicates (within "season" or "among seasons") and open populations In our paper we develop basic IDS models, with or without availability, for a static population.Incorporation of temporal replicates is the subject of future work (see also the Discussion of our paper).In principle, such extensions can be achieved by integrating the appropriate more complex likelihoods in the same manner as done in this paper.For example, Chandler et al. (2011) for repeated counts (or DS data) collected over time subject to temporary emigration, or the fully dynamic model of Dail and Madsen (2010).To the extent that these models have been described previously, there are no technical challenges to this generalized IDS model, although the code to implement them may be challenging even for specific cases let alone for general application.
Section S3: Frequentist operating characteristics of the IDS estimators in Simulations 1-4 Table S1: Simulation 1: Absolute bias, % relative bias, and 95% CRI coverage of the estimators (posterior means) of the half-normal scale parameter (sigma,  ) in the distance-sampling (DS) and the point count (PC) parts of the data, and of the site-level abundance (N) at the point-count sites under IDS models 1 and 2.An asterisk means that in this statistic all sites with true zero abundance were discarded due to division by zero.Sample size is 1000 sets of estimates from converged runs of the MCMC algorithm for each.Compare also with Figure 1  Estimates based on a run of JAGS with 4 chains of 120000 iterations, adaptation = 10000 iterations, burn-in = 60000 iterations, and thin rate = 60, yielding 4000 total samples from the joint posterior.

( 1 )
Distance-sampling (DS) data , ds ij y , possibly with truncation distance ds i b and survey duration ds i t , where j indexes J distance classes, and where ,. national BBS or bird Atlas schemes.(3) Detection/nondetection (DND) data dnd i y , indicating the observed presence or absence of a species during a point-location survey of duration dnd i t out to an optional truncation distance dnd i b , as they are similarly produced frequently by biological surveys.
availability ( ) and perceptibility ( ds p ) are the two components of detection probability(Marsh and Sinclair 1989;Nichols et al. 2009).Perceptibility will primarily be a function of distance and is estimated from distance data by integrating out to distance i b a suitable detection function such as a half-normal with scale parameter  .In turn, that truncation distance i b defines the survey area, which for a point count survey (assuming perfect detection) is 2 ii Ab  = .
combination of data.These likelihoods can be maximized numerically to obtain MLEs, or we can place priors on their parameters and use MCMC methods to obtain Bayesian posterior inferences.See Appendix S1 for a conceptual outline of IDS models and of how they conceptualize PC and DND data as the outcome of a "latent" distance sampling observation process.

Fig. 3 :
Fig. 3: Sampling distributions of estimators of density (intercept and slope of a continuous

Figure S1 :
Figure S1: Schematic of three types of point-based surveys in a hypothetical landscape, showing how the same underlying process, combined with different survey methods, gives rise to the three data types that we combine in an IDS model.

Figure S2 :
Figure S2: Another schematic of three types of point-based surveys in the same hypothetical landscape, this time emphasizing two typical features of real data sets fused in an IDS model: different survey locations and detection radii/truncation distances.
(2013),Farr et al. (2021),and Schmidt  et al. (2022).But although this structural part of the model is shared by all data types, different detection function parameters can be estimated for each data type, as we show in Simulation 1 of our paper.The truncation distance is shown by the vertical dotted line in FigureS3.For all three data types, the average per-individual detection probability in the area associated with the distance Distance models for observed data The typical observation models in an IDS, and those in the new IDS() function in R package unmarked, are a multinomial for the vector of binned distance counts y for the distance sampling data(Royle et al. 2004), a binomial (as in a binomial mixture model; Royle 2004) for point counts, and a Bernoulli distribution (as in the model of Royle & Nichols 2003) for detection/nondetection data.Distance sampling   ~ Multinomial(  , ) Point count   ~ Binomial(  , ̅ ) Detection/nondetection   ~ Bernoulli(1 − (1 − ̅ )   )

Table S2 :
Simulation 1B: Absolute bias, % relative bias, and 95% CRI coverage of the estimators (posterior means) of the half-normal scale parameter (sigma,  ) in the distance-sampling (DS) and the point count (PC) parts of the data, of the site-level abundance (N) at the point-count sites, and of mean density under the IDS model 1 (which combines DS + PC data).Asterisk: in this statistic all sites with zero true abundance were lost due to division by zero.Sample size is 1000 sets of estimates from converged runs of the MCMC algorithm.Compare also with FigureS1above.

Table S3 :
Simulation 2: Absolute bias, % relative bias, and 95% CRI coverage of the estimators (MLEs) in IDS1 combining distance-sampling (DS) and point count (PC) data and where the covariate effects in the detection functions of the DS and the PC data differ (Simulation 2a).Left: independent covariates are simulated in the density and the detection function parts of the model; right: a single covariate affects density as well as the two detection functions (Simulation 2b).Sample size is 977 (left) and 934 (right) valid model fits; see main text for more explanation.Compare also with Figure2in main paper.

Table S4 :
Simulation 3: Absolute bias, % relative bias, and 95% CRI coverage of the estimators (MLEs) of the shared intercept and the slope of density (  ), and of the scale parameter (sigma,  ) in the distance-sampling (DS) and the point count (PC) parts of the data under the IDS model 1 (which combines DS + PC data), when between 1 and 100 DS sites are added in the analysis of 200 PC sites.Original sample size is 1000 data sets at each added value of DS sites.Compare also with Figure3in main paper.