Hélène Labach
email: hlabach@miraceti.org

Caroline Azzinari

Maxime Barbier

Cathy Cesarini

Boris Daniel

Léa David

Frank Dhermain

Nathalie Di-Méglio

Benjamin Guichard

Julie Jourdan

Valentin Lauret

Nicolas Robert

Marine Roul

Nicolas Tomasi

Olivier Gimenez

Distribution and abundance of common bottlenose dolphin (Tursiops truncatus) over the French Mediterranean continental shelf

Keywords: Abundance, bottlenose dolphin, capture-recapture, distribution, French Mediterranean Sea, photo-identification, Tursiops truncatus

The common bottlenose dolphin (Tursiops truncatus) subpopulation in the Mediterranean is listed as vulnerable by the International Union for Conservation of Nature. This species is strictly protected in France and the designation of Special Areas of Conservation (SAC) is required under the EU Habitats Directive (92/43/EEC). However, little information is available about the structure and dynamics of bottlenose dolphins in French Mediterranean waters. We collected photo-identification data over the whole French Mediterranean continental shelf year round between 2013 and 2015. We sighted 151 groups of bottlenose dolphins allowing the individual photo-identification of 766 animals. The encounter rate distribution showed the presence of bottlenose dolphins over the whole continental shelf year round. We estimated for the first time, using capture-recapture methods, the size of this bottlenose dolphin population at 2,350 individuals (95% credible interval 1,827-3,135). Our results were used in support of the designation of a new dedicated SAC in the Gulf of Lion and provide a baseline for the bottlenose dolphin monitoring in the French Mediterranean waters in the context of the Marine Strategy Framework Directive.

Introduction

The common bottlenose dolphin (Tursiops truncatus, Montagu, 1821; hereafter bottlenose dolphin) is considered as a regular species in the Mediterranean Sea [START_REF] Reeves | The status and distribution of Cetaceans in the Black sea and Mediterranean sea[END_REF]. It has been observed along most of the Mediterranean coast [START_REF] Bearzi | Ecology and conservation of common bottlenose dolphins Tursiops truncatus in the Mediterranean Sea[END_REF], most often over the continental shelf [START_REF] Gannier | Summer distribution and relative abundance of delphinids in the Mediterranean Sea[END_REF][START_REF] Gnone | Distribution, abundance, and movements of the bottlenose dolphin (Tursiops truncatus) in the Pelagos Sanctuary MPA (north-west Mediterranean Sea)[END_REF][START_REF] Di Sciara | Cetaceans in the central mediterranean sea: Distribution and sighting frequencies[END_REF], even though groups have also been observed offshore [START_REF] Laran | Seasonal distribution and abundance of cetaceans within French waters-Part I: The North-Western Mediterranean, including the Pelagos sanctuary[END_REF]. Both resident populations and transient individuals have been reported [START_REF] Gnone | Distribution, abundance, and movements of the bottlenose dolphin (Tursiops truncatus) in the Pelagos Sanctuary MPA (north-west Mediterranean Sea)[END_REF]. The Mediterranean bottlenose dolphin subpopulation is genetically differentiated from populations inhabiting the contiguous eastern North Atlantic and the Black Sea and is structured into a Western and an Eastern population, corresponding to habitat boundaries [START_REF] Natoli | Habitat structure and the dispersal of male and female bottlenose dolphins (Tursiops truncatus)[END_REF].

The Mediterranean bottlenose dolphin subpopulation is considered as vulnerable on the IUCN (International Union for Conservation of Nature) Red List. It is listed in Annex II of the Washington Convention on International Trade in Endangered Species, in Appendix II of the Bern Convention for the Conservation of European Wildlife and Natural Habitats, in and compared with reference values through standardized long-term monitoring [START_REF] Cairns | A Proposed Framework for Developing Indicators of Ecosystem Health for the Great Lakes Region[END_REF][START_REF] Dale | Challenges in the development and use of ecological indicators[END_REF].

In France, the monitoring program set up for the implementation of the European Marine Strategy Framework Directive (2008/56/EC; MSFD) recommends specific monitoring of resident coastal populations of marine mammal species using photo-identification, including bottlenose dolphins. Bottlenose dolphins can be identified individually from their natural markings, and the resulting longitudinal sighting histories of individuals can be analyzed with capture-recapture (CR) models [START_REF] Hammond | Mark-Recapture[END_REF][START_REF] Hammond | Individual Recognition of Cetaceans : Use of Photo-Identification and Other Techniques to Estimate Population Parameters[END_REF][START_REF] Rosel | Photo-identification Capture-Mark-Recapture Techniques for Estimating Abundance of bay, Sound and Estuary Populations of Bottlenose Dolphins along the U.S. East Coast and Gulf of Mexico: A Workshop Report[END_REF] to estimate population abundance and survival. Photo-identification data have been used widely to monitor populations of bottlenose dolphins (e.g. [START_REF] Defran | Occurrence, distribution, site fidelity, and school size of bottlenose dolphins (Tursiops truncatus) off San Diego, California[END_REF][START_REF] Gnone | Distribution, abundance, and movements of the bottlenose dolphin (Tursiops truncatus) in the Pelagos Sanctuary MPA (north-west Mediterranean Sea)[END_REF][START_REF] Karczmarski | Matrix photo-identification technique applied in studies of free-ranging bottlenose and humpback dolphins[END_REF][START_REF] Louis | Social structure and abundance of coastal bottlenose dolphins , Tursiops truncatus, in the Normano-Breton Gulf, English Channel[END_REF][START_REF] Shane | Ecology, behavior and social organization of the bottlenose dolphin: a review[END_REF].

In French Mediterranean waters, several studies on bottlenose dolphins have been conducted since the 1990s, and many of them were based on photo-identification [START_REF] Bompar | Estimation of the Corsican population of bottlenose dolphins Tursiops truncatus: is there a real conflict with fishermen[END_REF][START_REF] Dhermain | First evidence of the movement of a bottlenose dolphin Tursiops truncatus between Corsica and Hyères archipelago, south-eastern France[END_REF][START_REF] Labach | Suivi de la population de grands dauphins Tursiops truncatus le long des côtes provençales (Méditerranée nord-occidentale)[END_REF][START_REF] Labach | Suivi des Grands Dauphins (Tursiops truncatus) et Dauphins de Risso (Grampus griseus) sur le secteur des îles d' Hyères en 2009[END_REF][START_REF] Ripoll | First summer population estimate of bottlenose dolphins along the north-western coasts of the occidental Mediterranean basin[END_REF]. The knowledge of the population structure, ecology, and dynamics remains poor and unequal, in part because these studies were limited to small areas over short periods.

In this study, we conducted the first large-scale bottlenose dolphin photo-identification survey in French Mediterranean waters. Standardized photo-identification data were collected throughout the French Mediterranean continental shelf year round over two years through a standard protocol by a network of organizations. The objectives of our study were to evaluate the distribution of bottlenose dolphins and to provide the first population abundance estimate over the French continental shelf.

Methods

Study area

The French Mediterranean waters contain high diversity and a richness of habitats and seabed.

The Gulf of Lion, from the Spanish border to Marseille, is a vast continental shelf limited to the north by a sandy and lagoon coastline and to the south by a broad slope cut by numerous canyons. The Corso-Liguro-Provençal basin (Riviera and west coast of Corsica) presents a rocky coastline prolonged by a narrow continental shelf quickly giving way to an abrupt slope, cut by deep canyons. To the east of Corsica, the reliefs are shallower with a wider continental shelf. The Corso-Liguro-Provençal basin and the Gulf of Lion are highly productive areas, attracting a high diversity of species [START_REF] Ribera Dalcaì | On the trophic regimes of the Mediterranean Sea: a satellite analysis[END_REF].

The study area covers the continental shelf of the French Mediterranean waters between the coast and the 500 m isobath, bounded by the Spanish border to the west, the Italian border to the east, and includes the whole Corsican coastline (Fig. 1). The overall study area covers 24,481 km 2 and was divided into three regions based on geographic and topographic characteristics: Gulf of Lion (14,731 km 2 ), Riviera (2,866 km 2 ), and Corsica (6,884 km 2 ). To ensure homogeneous sampling over the whole study area, each region was divided into subregions of similar area, covering on average 2,500 (+/-500) km 2 (4 in Gulf of Lion, 2 in Riviera and 3 in Corsica) and assigned to 5 local organizations associated with marine mammal monitoring (BREACH, CARI, EcoOcéan Institut, GECEM, and Parc naturel régional de Corse) according to their location and usual study areas, with each organization covering one or two subregions.

Data collection

The objective of the survey was to conduct four days of boat-based survey effort per season within each subregion. We defined the survey effort as the length (in km) of track actively traveled, prospecting the area with naked eyes by three observers. Surveys were conducted only in favorable weather conditions (Beaufort scale 0-3 and good visibility). Seasons were defined as spring (March 22-June 21), summer (June 22-September 21), autumn (September 22-December 21), and winter (December 22-March 21). This survey design was conducted twice between 2013 and 2015 using small sailing and motor boats. Routes were planned on a variable basis, according to weather conditions, to optimize the study area sampling coverage.

All partners applied a standard common protocol using a digital application for the data collection specifically designed with Cybertracker (https://www.cybertracker.org/), systematically recording survey tracks with a GPS receiver. When we encountered a group of bottlenose dolphins, we recorded the position of first contact, group size, and composition along with the group's main activity. We defined a group as all dolphins seen with naked eye during the sighting. The estimated group size is the estimated number of individuals observed or photo-identified when the latter figure is greater than the estimated one.

Whenever possible, we took photographs with a digital reflex camera of both sides of the dorsal fins of all individuals of the group, regardless of their markings.

Photo-identification

We identified individuals using natural marks: scars, nicks, and scratches on their dorsal fins (Würsig & Jefferson, 1990;Würsig & Würsig, 1977). The best photographs of each side of each individual at every encounter were selected and scored according to their quality (Q1 for good, Q2 for medium, and Q3 for bad) and the distinctiveness of animals (D1 for wellmarked, D2 for moderately marked and D3 for poorly marked) [START_REF] Berrow | Abundance estimate of bottlenose dolphins (Tursiops truncatus) in the Lower River Shannon candidate special area of conservation, Ireland[END_REF][START_REF] Ingram | The ecology and conservation of bottlenose dolphins in the Shannon Estuary, Ireland[END_REF] to generate catalogs of identified dolphins. Catalog and associated data were incorporated into a common database which was uploaded onto the international web database INTERCET (http://www.intercet.it/). Each partner compared its own catalog with all the others compiled within the project, and a global matching was conducted by one additional scientist, ensuring that matching was conducted by at least three people, minimizing bias due to false positive and negative matching. The matching process led to three regional catalogs and one global one. A matrix of sighting histories was associated with each catalog gathering all the sighting data and information of all the dolphins identified. For the analyses, we used only medium and good quality photographs (quality scores = 1 or 2) of moderately and well-marked individuals (distinctiveness score = 1 or 2).

The proportion of moderately and well-marked individuals was evaluated as the number of individuals identified during the project, considering each side (right or left side) separately, that scored D1 or D2 (removing all pictures scored D3) among individuals identified by at least one photograph scored Q1 or Q2 (removing all pictures scored Q3) (see also the Abundance estimation section).

Distribution

We calculated the encounter rate (ER) as the number of sightings per km of effort traveled in each region and within each 5'x5' cell of the Marsden grid WGS 84. All maps and spatial analyses were done in R 3.5.0 (R Core [START_REF] Team | R: A language and environment for statistical computing[END_REF].

Abundance estimation

To estimate the abundance of bottlenose dolphins occurring within the study area, we fitted CR models to the photo-identification data [START_REF] Hammond | Individual Recognition of Cetaceans : Use of Photo-Identification and Other Techniques to Estimate Population Parameters[END_REF]. We defined a capture when an individual was identified using photo-identification, and a recapture as the resighting of an individual already identified during the project. Because some individuals were sighted both in Gulf of Lion and in the Riviera during the study period, we performed three separate analyses, corresponding to the sightings made in the Gulf of Lion, in the Riviera, and along the continental coast (Gulf of Lion plus Riviera). We did not pursue CR analyses with the Corsican sightings because of insufficient recaptures (Table1). We used Bayesian closed population models [START_REF] Mcclintock | Multimark: An R package for analysis of capture-recapture data consisting of multiple "noninvasive" marks[END_REF] to estimate abundance while accounting for a capture probability less than one. We considered the eight seasons as our capture occasions. The main assumptions underlying closed population CR models are 1) the population is demographically closed (i.e., natality and mortality events do not occur) during the study period; 2) all individuals are correctly identified at each capture occasion and 3) the marks are considered permanent.

To fit CR models, we used the package multimark [START_REF] Mcclintock | Multimark: An R package for analysis of capture-recapture data consisting of multiple "noninvasive" marks[END_REF] in R (R CoreTeam, 2018), which implements Markov Chain Monte Carlo (MCMC) simulations. We performed an analysis for each data set, that is the Gulf of Lion, the Riviera, and the continental coast (Gulf of Lion plus Riviera). For each analysis, we fitted eight models, including a model with constant detection probability (M0), a model with time variation in the detection probability (Mt), a model with heterogeneity in the detection probability (Mh) under the form of an individual random effect, a model with behavioral response in the detection probability (Mb), and combinations of these effects (Mbt, Mht, Mbh, Mbht). To determine the model best supported by the data, we calculated posterior model probabilities [START_REF] Barker | Bayesian multimodel inference by RJMCMC: A Gibbs sampling approach[END_REF].

More precisely, we relied on an extension of the standard Bayes theorem where the posterior distribution of all parameters is now defined over both the parameter and model space. In addition to posterior summaries for parameters (abundance and detection probabilities), we also obtained the posterior probability for each model obtained as the proportion of the time the MCMC simulations spend in each model. Because we used only moderately-and wellmarked individuals (assumed to be adults) in the CR analyses, the abundance, including poorly marked individuals (juveniles and neonates), was obtained by correcting the CR abundance by the proportion of moderately and well-marked individuals (Williams, Dawson a Slooten, 1993) estimated as : Δ = nQD12/nQ12 where Δ is the proportion of moderately and well-marked individuals, nQ12 is the number of individuals identified after removing bad quality photos (Q3) and nQD12 is the number of individuals identified after removing bad quality photos (Q3) and poorly marked individuals (D3) evaluated on the global data set (all captures included) for each data set. We estimated this proportion for left-and right-side photographs then calculated the average of these two proportions for each data set. The standard error was obtained with the formula 𝑝̂ (1 -𝑝) / n for a proportion estimate 𝑝̂ where n is the sample size.

In practice, the Bayesian approach using MCMC made it easy to propagate the uncertainty in the proportion of moderately and well-marked individuals. We divided each MCMC value drawn in the posterior distribution of abundance by a random draw from a normal distribution with mean Δ and standard deviation the standard error, therefore providing the posterior distribution of the corrected CR abundance. We reported posterior means and 95% credible intervals for abundance and detection probabilities.

Results

Survey effort

We traveled 21,464 km in survey effort. The distribution of the effort between the 3 regions was heterogeneous with a high coverage of Riviera but low coverage of Corsica and the offshore areas of Gulf of Lion. Summer was the best prospected season, autumn and winter being less prospected in the three regions (Fig. 2).

Distribution

We sighted bottlenose dolphins across the study area in all seasons (Fig. 3). Global ER was higher in Corsica (0.012) than in Gulf of Lion (0.007) and in Riviera (0.003) (Table 2). In Riviera, ER appeared higher in spring (Fig. 3), while in Gulf of Lion and Corsica, ER appeared higher in summer (Fig. 3).

Sightings and photo-identification

We sighted 151 groups of bottlenose dolphins during the project. Group size was highly variable in the three regions, mean group size was similar in Riviera (15.7 SD 10.3) and Gulf of Lion (16.6 SD 13.2) and lower in Corsica (5.3 SD 4.5) (Table 2).

After photos scoring and sorting, 766 different moderately and well-marked dolphins were identified on the basis of good and medium quality photos (Table 2), of which 30% were observed more than once during the project. The percentage of individuals recaptured was higher in Gulf of Lion (29%) than in Riviera (19%) and in Corsica (18%). Six individuals were sighted in both the Gulf of Lion and Riviera, while no recaptures were made between continental and Corsican coast.

Abundance estimates

We excluded 15% of the 1,705 photographs from the analyses because of their low quality (Q3). The proportion of moderately-and well-marked individuals was 0.73 in Riviera (SE = 0.02), 0.84 (SE = 0.01) in Gulf of Lion and 0.86 (SE = 0.01) in the whole continental coast (Riviera and Gulf of Lion). Many dolphins (68% in continental coast) were seen only once.

The maximum number of captures was 6 for two dolphins (Table 1).

The model best supported by the Gulf of Lion and whole continental coast data included heterogeneity and temporal variation in the detection probability (posterior probability = 0.70 and 0.78 respectively), while the model best supported by the Riviera data was the one considering heterogeneity, a behavioral response and temporal variation in the detection probability (posterior probability = 1). Detection probabilities varied between 0.02 (0.01-0.03) and 0.16 (0.11-0.22) for Gulf of Lion, between 0.07 (0.01-0.19) and 0.44 (0.19-0.68) for Riviera, and between 0.02 (0.02-0.03) and 0.15 (0.11-0.19) for the whole continental coast.

Abundance (corrected by the proportion of moderately-and well-marked individuals) of population was estimated at 223 (152-385) individuals in Riviera,2,231 (1,175) in Gulf of Lion and 2,350 (1,135) along the whole continental coast.

Discussion

Our study provides the first large-scale dedicated photo-identification survey for the bottlenose dolphin in the French Mediterranean waters. We demonstrate the power of a collaborative and coordinated survey to study a mobile species at population's scale. Our results show that the whole continental shelf is frequented by bottlenose dolphins, including the entire Gulf of Lion, all year round and provides the first abundance estimate of bottlenose dolphins frequenting the French continental shelf of Riviera and Gulf of Lion.

The effort of 21,464 km covered 87% of the study area. We found heterogeneity in this effort, mainly between Corsica and the continental coast, which we explain by a later start of the survey in Corsica and difficult survey conditions during the study. The results obtained in Gulf of Lion and Riviera show that the survey effort is sufficient to provide consistent estimates of abundance and distribution.

The global encounter rate (0.007) was higher than the encounter rates obtained with the program "Surveillance Aérienne de la Mégafaune Marine" (SAMM) (0.0041 in winter and 0.0028 in summer) [START_REF] Laran | Seasonal distribution and abundance of cetaceans within French waters-Part I: The North-Western Mediterranean, including the Pelagos sanctuary[END_REF]. This difference could be explained by the fact that this comprehensive aerial survey of marine megafauna conducted by the French Biodiversity Agency in 2011 and 2012 over the whole French Exclusive Economic Zone (EEZ), encompasses continental shelf, slope, and oceanic waters [START_REF] Laran | Seasonal distribution and abundance of cetaceans within French waters-Part I: The North-Western Mediterranean, including the Pelagos sanctuary[END_REF]. Another likely reason might be that aerial surveys consistently have lower encounter rates than boat-based surveys because of the order of magnitude difference in survey speed. The ER in Riviera (0.003) and in Corsica (0.012) were also higher than the maximum ER obtained by [START_REF] Gnone | Distribution, abundance, and movements of the bottlenose dolphin (Tursiops truncatus) in the Pelagos Sanctuary MPA (north-west Mediterranean Sea)[END_REF]Gnone et al. ( ) between 1994Gnone et al. ( and 2007 in in Provence (ER = 0.0006) and in Corsica (ER = 0.0086), which suggest an increase in dolphin abundance in these two regions.

The distribution of ER showed that bottlenose dolphins were present over the entire French Mediterranean continental shelf all year round. The higher ER in summer in the Gulf of Lion and Corsica was consistent with the results of the SAMM survey, which despite showing higher ER in winter than summer in the global EEZ, also showed contrasting seasonal distributions, with encounters concentrating in coastal areas of the Gulf of Lion and Corsica in the summer [START_REF] Laran | Seasonal distribution and abundance of cetaceans within French waters-Part I: The North-Western Mediterranean, including the Pelagos sanctuary[END_REF]. These results suggest a seasonal migration of bottlenose dolphins between offshore waters in winter to coastal waters in summer, especially in Gulf of Lion and Corsica. The sighting of dolphins both in Riviera and Gulf of Lion also points towards some eastward and westward movements. No movement between the continental areas and Corsica was observed during the project. In previous studies [START_REF] Gnone | Distribution, abundance, and movements of the bottlenose dolphin (Tursiops truncatus) in the Pelagos Sanctuary MPA (north-west Mediterranean Sea)[END_REF], 5 individuals were identified both in Corsica and the continental coast, highlighting that some dolphins perform long distance travels. The high percentage of dolphins captured only once (70%) during the project can be explained by the short period and the large study area of the project decreasing the recapture probability. It may also highlight that an important proportion of bottlenose dolphins sighted over the French Mediterranean continental shelf are transient animals coming from remote areas, as suggested by the seasonal differences in the ER and the movements identified, whereas other animals are resident, as suggested by sightings of some individuals all year round for more than 20 years (authors unpublished data). Pursuing photoidentification at this scale over the long term will allow further exploration of and characterizing residency patterns. In their study, [START_REF] Gnone | Distribution, abundance, and movements of the bottlenose dolphin (Tursiops truncatus) in the Pelagos Sanctuary MPA (north-west Mediterranean Sea)[END_REF] identified two subpopulations inside the Pelagos Sanctuary coinciding with the national boundaries between French and Italian territories. The identification of distinct units and the characterization of connections between them along the French continental coast is the object of ongoing work using population genetic and social structure analyses based on photo-identification and biopsy data collected during the present study. Sharing photo-identification catalogs and associated metadata through the INTERCET platform will make possible the characterization of bottlenose dolphin movements and social structure at wider Mediterranean scale. The higher percentage of poorly marked individuals in Riviera (27%) suggests a higher percentage of immature dolphins in this region than in Gulf of Lion (16%).

The robust estimation of abundance relies on the validation of CR model assumptions. The two-year sampling period and the fact that newborns were observed in the study area suggest that assumption 1 of the capture-recapture model is likely to have been violated. We expect little underestimation bias in the abundance estimates because of birth occurring during the study, as we considered adults only, and they have high survival probability. Mortality occurring during the study could also lead to abundance overestimation, nevertheless we expect little bias considering the relative short period for this long-lived species and the high survival probability. Assumptions 2 and 3 are ensured by the fact that only moderately and well-marked individuals with medium and good-quality photographs were included in the analysis. Also, if the marks evolve, the short sampling period would allow recognition of the animals.

The population abundance along the continental coast was higher than the estimates of the only previous census dedicated to bottlenose dolphins in the same area, which estimated based on observed count (not corrected by imperfect detection) the number of bottlenose dolphins between 200 and 209 in the Gulf of Lion and 16 in Provence [START_REF] Ripoll | First summer population estimate of bottlenose dolphins along the north-western coasts of the occidental Mediterranean basin[END_REF]. These figures are not inconsistent with our abundance estimates which accounted for imperfect detection by correcting the observed counts by the estimated detection probability.

Our abundance estimates are coherent with the results obtained from the program SAMM with the distance sampling methodology, which estimated the absolute abundance of bottlenose dolphins in the Gulf of Lion at 63 (95% CI 17-241) in winter and at 1,331 (95% CI 466-3,805) in summer and over the continental slope inside the French EEZ waters and including some Italian and Spanish waters at the eastern and western borders at 1,795 (95% CI 190) in winter and 10 (95% CI 3-30) in summer [START_REF] Laran | Seasonal distribution and abundance of cetaceans within French waters-Part I: The North-Western Mediterranean, including the Pelagos sanctuary[END_REF].

Further efforts should be implemented to complete the survey planned for this study in Corsica in order to provide consistent estimates for this region also.

Implications for conservation

Our study provides an operational framework as well as a baseline for the implementation of a long-term large-scale monitoring of bottlenose dolphin population in the French Mediterranean waters. We therefore recommend pursuing the monitoring initiated in this study for the long-term, taking into consideration the evaluation of the efficacy of the survey design, on which we're currently working, to allow the identification of trends in the population as required for the surveillance program of the MSFD.

We shared the data on the international webGIS platform INTERCET (http://www.intercet.it/), which will allow enlarging the study of this species beyond French boundaries to the basin and Mediterranean scale.

The results of our study, together with those from the SAMM survey [START_REF] Laran | Seasonal distribution and abundance of cetaceans within French waters-Part I: The North-Western Mediterranean, including the Pelagos sanctuary[END_REF], led to an update of the Mediterranean bottlenose conservation status in the national IUCN Red List, which was changed from "vulnerable" in 2009 to "near threatened" in 2017 because of the improved knowledge. They also contributed to the update of the Mediterranean subpopulation status initiated in 2020. Our demonstration of the presence of bottlenose dolphins in the entire Gulf of Lion led France to submit the designation of a dedicated offshore SAC encompassing the whole Gulf of Lion continental shelf beyond the territorial waters and to the recognition of this area as an important marine mammal area (IMMA) for bottlenose dolphins (https://www.marinemammalhabitat.org/imma-eatlas/). Our results will also contribute to updating the ACCOBAMS bottlenose dolphin conservation plan. of moderately-and well-marked individuals identified on the basis of good and medium quality photos, and number of recaptured individuals (sighted more than once) in each region. The bathymetry is also displayed on the map.

Tables

Figure legends Figure 1 :

 legends1 Figure legends

Figure 2 :

 2 Figure 2: Seasonal distribution of survey effort (number of kilometers actively traveled per

Figure 3 :

 3 Figure 3: Seasonal distribution of bottlenose dolphins over French Mediterranean waters

  Fig.1

Table 1 .

 1 Distribution of individuals per number of captures. Number of moderately-and well-marked individuals identified (on the basis of good and medium quality photos) 1, 2, 3, etc. times in each data set. Continental coast refers to Riviera plus Gulf of Lion

		1	2	3	4	5	6	Total
	Corsica	78	16	0	1	0	0	95
	Riviera	79	9	5	3	1	0	97
	Gulf of Lion	411	100	51	15	1	2	580
	Continental coast 458	123	61	21	6	2	671

Table 2 :

 2 Sightings and photo-identification of bottlenose dolphins

		Sightings Encounter rate	Mean Group size	Identified	Recaptured
				(SD)	individuals	individuals
	Corsica	41	0.012	5.3 (4.5)	95	17 (18%)
	Riviera	18	0.003	15.7 (10.3)	97	18 (19%)
	Gulf of Lion 92	0.007	16.6 (13.2)	580	169 (29%)
	Global	151	0.007	13.6 (12.5)	766	230 (30%)
	Number of sightings, encounter rates, mean group size and standard deviation (SD), number

Acknowledgements

We acknowledge MAVA Foundation, the French Biodiversity Agency and the Pelagos Sanctuary, who financially supported the GDEGeM project and also thank Fondation de France for funding through the INTERACT project. We thank all the partners for their involvement in the data collection: BREACH, CARI Corse, EcoOcéan Institut, GECEM, and Parc naturel régional de Corse, as well as all the people who participated to the survey. We also acknowledge Guido Gnone and Michela Bellingeri for their support in the use of the INTERCET platform.

 (1990). Methods of photo-identification for small cetaceans.

Report of the


Daniel Turek

Claudia Wehrhahn
email: cwehrhah@ucsc.edu

Olivier Gimenez
email: olivier.gimenez@cefe.cnrs.fr

Bayesian non-parametric detection heterogeneity in ecological models

Keywords: Bayesian non-parametrics, Capture-recapture, Detection heterogeneity, Markov chain Monte Carlo, nimble, nimbleEcology, Occupancy models, Statistical ecology

Detection heterogeneity is inherent to ecological data, arising from factors such as varied terrain or weather conditions, inconsistent sampling effort, or heterogeneity of individuals themselves. Incorporating additional covariates into a statistical model is one approach for addressing heterogeneity, but is no guarantee that any set of measurable covariates will adequately address the heterogeneity, and the presence of unmodelled heterogeneity has been shown to produce biases in the resulting inferences. Other approaches for addressing heterogeneity include the use of random effects, or finite mixtures of homogeneous subgroups. Here, we present a non-parametric approach for modelling detection heterogeneity for use in a Bayesian hierarchical framework. We employ a Dirichlet process mixture which allows a flexible number of population subgroups without the need to pre-specify this number of subgroups as in a finite mixture. We describe this non-parametric approach, then consider its use for modelling detection heterogeneity in two common ecological motifs: capture-recapture and occupancy modelling. For each, we consider a homogeneous model, finite mixture models, and the non-parametric approach. We compare these approaches using two simulation studies, and observe the non-parametric approach as the most reliable method for addressing varying degrees of heterogeneity. We also present two real-data examples, and compare the inferences resulting from each modelling approach. Analyses are carried out using the nimble package for R, which provides facilities for Bayesian non-parametric models.

Introduction

Inferring species distribution and demography are two key questions in ecology [START_REF] Begon | Ecology: From Individuals to Ecosystems[END_REF]. However, addressing these questions is challenging when studying animals and plants in their natural environmental because of the inherent imperfection in the detection process of species or individuals. When ignored, this issue of imperfect detection can lead to biased estimates of species distribution or survival [START_REF] Gimenez | The risk of flawed inference in evolutionary studies when detectability is less than one[END_REF][START_REF] Guillera-Arroita | Modelling of species distributions, range dynamics and communities under imperfect detection: advances, challenges and opportunities[END_REF]. To cope with imperfect detection, occupancy models [START_REF] Mackenzie | Occupancy Estimation and Modeling: Inferring Patterns and Dynamics of Species Occurrence[END_REF] and capture-recapture models [START_REF] Mccrea | Analysis of Capture-Recapture Data[END_REF] were developed to provide unbiased estimates of species range and individual survival, with numerous applications in all fields of ecology.

As with any statistical model, inferences made from occupancy and capture-recapture analyses rely on assumptions which must be satisfied at least to a reasonable degree.

Common to occupancy and capture-recapture models is the assumption of homogeneity of the detection process, which asserts that there is no unmodelled heterogeneity in species and individual detection probabilities. If ignored, heterogeneity in detection leads to flawed inference in occupancy and capture-recapture models [START_REF] Royle | Site occupancy models with heterogeneous detection probabilities[END_REF][START_REF] Gimenez | Individual heterogeneity and capture-recapture models: what, why and how?[END_REF]. Detection heterogeneity can be due to heterogeneous sampling effort, variation in animal abundance or behavior, site characteristics or even on account of varying observer skills. Ideally, covariates could be measured and incorporated in ecological models to account for detection heterogeneity. However, unexplained variation may still remain, or measuring the relevant covariates may simply be impossible.

When unmodelled heterogeneity exists, it can be accommodated using finite mixtures in which discrete latent variables are used to assign sites or individuals to mixture components (i.e., uniform groups) each characterized by group-specific parameters [START_REF] Royle | Site occupancy models with heterogeneous detection probabilities[END_REF][START_REF] Pledger | Open capture-recapture models with heterogeneity: II. Jolly-Seber model[END_REF][START_REF] Pledger | Open capture-recapture models with heterogeneity: I. Cormack-Jolly-Seber model[END_REF][START_REF] Louvrier | Accounting for misidentification and heterogeneity in occupancy studies using hidden Markov models[END_REF]. In simulation studies, finite mixtures were successful in decreasing bias in occupancy and survival probability estimates that was introduced by heterogeneity in the detection process [START_REF] Pledger | The performance of mixture models in heterogeneous closed population capture-recapture[END_REF][START_REF] Louvrier | Accounting for misidentification and heterogeneity in occupancy studies using hidden Markov models[END_REF]. However, from a practical perspective, the issue remains of selecting the number of mixture components for real-life data analyses, which is not straightforward [START_REF] Cubaynes | Assessing individual heterogeneity using model selection criteria: how many mixture components in capture-recapture models?[END_REF][START_REF] Pohle | Selecting the Number of States in Hidden Markov Models: Pragmatic Solutions Illustrated Using Animal Movement[END_REF].

Here, we propose a Bayesian non-parametric (BNP) approach to modeling heterogeneity in occupancy and capturerecapture models. BNP models provide a flexible approach that relaxes typical standard modeling assumptions, such as the choice of a specific parametric kernel in density estimation, or here the choice of a fixed number of groups in finite mixtures. Previous uses of a BNP approach in ecological models include the modelling of wildlife migration patterns [START_REF] Matechou | Modelling individual migration patterns using a Bayesian nonparametric approach for capture-recapture data[END_REF]Diana et al., na), the estimation of population size [START_REF] Manrique-Vallier | Bayesian population size estimation using dirichlet process mixtures[END_REF][START_REF] Dorazio | Modeling unobserved sources of heterogeneity in animal abundance using a dirichlet process prior[END_REF] and that of the probability of remaining in or vacating a given area [START_REF] Ford | Modelling latent individual heterogeneity in mark-recapture data with dirichlet process priors[END_REF].

Perhaps the most widely used BNP model is the Dirichlet process mixture (DPM) model [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF][START_REF] Ferguson | Prior distribution on the spaces of probability measures[END_REF][START_REF] Lo | On a class of Bayesian nonparametric estimates I: Density estimates[END_REF][START_REF] Escobar | Estimating normal means with a Dirichlet process prior[END_REF][START_REF] Escobar | Bayesian density estimation and inference using mixtures[END_REF], which is a mixture model with infinitely many components. DPM models are a suitable fit for addressing the inherent heterogeneity present in ecological models. We consider DPM models of Bernoulli distributions, with the distribution of the Bernoulli detection probability parameters arising from a Dirichlet process. A DPM model can be represented in different, yet equivalent manners. Two of them are the the Chinese Restaurant Process (CRP; [START_REF] Blackwell | Ferguson distributions via Pólya urn schemes[END_REF][START_REF] Pitman | Exchangeable and partially exchangeable random partitions[END_REF][START_REF] Pitman | Some developments of the Blackwell-MacQueen urn scheme[END_REF] and the stick-breaking (SB; [START_REF] Sethuraman | A constructive definition of Dirichlet prior[END_REF] representations.

The flexibility of BNP models is usually translated into a hierarchical model, which relies on Markov chain Monte Carlo algorithms to sample from the resulting posterior distribution. The implementation of these algorithms is usually computationally complex and demanding in standard Bayesian software such as WinBUGS or JAGS (e.g. [START_REF] Ohlssen | Flexible random-effects models using bayesian semi-parametric models: applications to institutional comparisons[END_REF] and often requires writing specific code to implement the BNP model and sampling algorithms [START_REF] Ford | Modelling latent individual heterogeneity in mark-recapture data with dirichlet process priors[END_REF]. Recently, the nimble R package introduced non-parametric functionality to address these difficulties, by supporting the use of non-parametric DPM models. Specialized functions, distributions, and samplers are provided for both the CRP and SB representations. [START_REF] Wehrhahn | Bayesian nonparametric mixture models using nimble[END_REF] contains additional details and examples of using these BNP modelling approaches in nimble.

Here, we focus on the CRP representation to address detection heterogeneity in capture-recapture and occupancy models. The article is organized as follows. In the next section, we present homogeneous capture-recapture and occupancy models, finite mixtures and non-parametric models formulated in a Bayesian framework as hierarchical models. In Section 3, we give the MCMC sampling scheme used to fit non-parametric models and introduce the nimble and nimbleEcology R packages. Section 4 gives the results of two simulation studies, using capture-recapture and occupancy data, which validate the ability of our approach to capture detection heterogeneity in occupancy and survival probabilities. Section 5 illustrates our method using data from a study of gray wolf (canis lupus), in trying to estimate occupancy and survival while accounting for detection heterogeneity. The final section provides general conclusions and discusses the potential of our approach.

Models

We consider three different approaches to handling heterogeneity. The first approach uses a homogeneous model, which disregards any heterogeneity and considers all individuals and sites to be identical. The second approach uses finite mixtures of population subgroups, where individuals or sites within each subgroup are identical, distinct population subgroups may differ in one or more characteristics, and the number of subgroups is pre-specified. The homogeneous model is a special case of the finite mixture model, where the population contains only a single subgroup.

The third approach uses a non-parametric representation for modeling individual and species heterogeneity, where the number of population subgroups, the assignment of individuals and sites into subgroups, and the characteristics of each subgroup are determined by the data at the time of model fitting. The non-parametric approach does not require pre-specifying a fixed number of population subgroups, but rather, the number of population subgroups itself is a model parameter. We now describe the three model formulations in detail, and provide specifications of ecological capturerecapture and occupancy models using each.

Homogeneous Models

In the homogeneous model, we assume that all N observed individuals of the population are identical in their characteristics. There is no heterogeneity between individuals or sites, and thus for any parameter θ of the population, all individuals or sites share a common value of θ.

Homogeneous Capture-Recapture Model

We consider a basic ecological capture-recapture model, for binary-valued observation data of N individuals or sites occurring over T observational time periods. We condition on the first observation of each individual occurring in time period t = 1, although it is straightforward to relax this assumption to allow first-detections to occur in other time periods. We parameterize the model in terms of survival probability between time periods φ, and probability of detection conditional on being alive p.

In subsequent models we introduce heterogeneity in p, but in the homogeneous model all individuals are characterized by the constant probability of detection p. We use a state-space formulation of the model, where binary states xi,t give the alive/dead status of individual i at time t, and using binary observation data yi,t for site i at time t. The homogeneous capture-recapture model is written:

φ ∼ Uniform(0, 1) p ∼ Uniform(0, 1) i = 1, . . . , N : xi,1 = 1 xi,t ∼ Bernoulli(φ • xi,t-1) for t = 2, . . . , T yi,t ∼ Bernoulli(p • xi,t) for t = 1, . . . , T

Homogeneous Occupancy Model

We consider a homogeneous static occupancy model for a total of N sites, each observed at T distinct sampling occasions. We parameterize the model in terms of the constant site occupancy probability ψ, and probability of detection conditional on site occupancy p. In subsequent models we introduce heterogeneity in p, but in the homogeneous model all sites share the probability of detection p. We use a statespace formulation of the model, where binary states zi give the true occupancy status of site i, and using binary observation data yi,t for the observation of site i at sampling occasion t. The homogeneous occupancy model is written:

ψ ∼ Uniform(0, 1) p ∼ Uniform(0, 1) i = 1, . . . , N : zi ∼ Bernoulli(ψ) yi,t ∼ Bernoulli(p • zi) for t = 1, . . . , T

Finite Mixture Models

In a finite mixture model specified as having K 2 distinct population subgroups, each individual or site is considered to be a member of any of the K subgroups with equal probability 1/K. We introduce a discrete indicator variable for each individual or site, gi for individual or site i, where gi denotes the "group" of individual or site i. We use independent discrete uniform prior distributions over the set {1, 2, . . . , K} for each gi, and gi = k indicates that individual or site i is a member of population subgroup k.

Furthermore, each of the K distinct subgroups may differ in one or more demographic characteristics. Again considering the demographic parameter θ, we introduce K model parameters θ1, θ2, . . . , θK , which are given independent and identical prior distributions. Then, all individuals or sites belonging to population subgroup k display θ = θ k . To avoid issues of "label exchanging" between groups, and hence a lack of model identifiability, we also impose the constraint that θ k θ k+1 , or that the ordered set of parameters {θ1, θ2, . . . , θK } is nondecreasing.

Finite Mixture Capture-Recapture Model

We generalize the homogeneous capture-recapture model given in section 2.1 to a K-group finite mixture model, to allow heterogeneity in the probability of detection. Specifically, we introduce K new model parameters p1, p2, . . . , pK , where p k represents the probability of detection for individuals in subgroup k. Each p k is given an independent Uniform(0, 1) distribution, and we impose the constraint on these parameters that p k p k for k < k to maintain identifiability.

Since gi gives the subgroup number which contains individual i, we use probability of detection pg i for individual i. Putting this together, the K-group finite mixture capturerecapture model is written as:

φ ∼ Uniform(0, 1) p k ∼ Uniform(0, 1) for k = 1, . . . , K p k p k for k < k i = 1, . . . , N : gi ∼ DiscreteUniform({1, . . . , K}) xi,1 = 1 xi,t ∼ Bernoulli(φ • xi,t-1) for t = 2, . . . , T yi,t ∼ Bernoulli(pg i • xi,t) for t = 1, . . . , T

Finite Mixture Occupancy Model

We similarly generalize the homogeneous occupancy model from section 2.1 to a K-group finite mixture model using parameters p1, p2, . . . , pK , where p k represents the probability of detection for sites in subgroup k. We use the same independent Uniform(0, 1) prior distribution for each p k , and impose the same constraint to ensure model identifiability. Thus, gi indicates the subgroup which contains site i, which therefore has probability of detection pg i . The K-group finite mixture occupancy model is written as:

ψ ∼ Uniform(0, 1) p k ∼ Uniform(0, 1) for k = 1, . . . , K p k p k for k < k i = 1, . . . , N : gi ∼ DiscreteUniform({1, . . . , K}) zi ∼ Bernoulli(ψ) yi,t ∼ Bernoulli(pg i • zi) for t = 1, . . . , T

Non-Parametric Models

Using a non-parametric approach, the pre-specification of a fixed number of subgroups is no longer required. Furthermore, our previous assumption that individual or site assignments to subgroups must follow a discrete uniform distribution is also relaxed. In a non-parametric model the number of subgroups in the population is considered unknown and is inferred from the data. Theoretically, there could be an infinite number of population subgroups, although in practice there will never exceed N subgroups. This means there cannot exist more population subgroups than the total number of observations. But so as long as there are fewer than N subgroups, then individuals or sites can probabilistically move out of an existing group and into to a newly created group, with its own distinct probability of detection.

As with the finite mixture model, the subgroup assignment structure is encoded in indicator variables gi, where gi = k indicates that individual i belongs to population subgroup k. Consider the following conditional distribution for gi, where g 1:(i-1) = (g1, . . . , gi-1): gi

| g 1:(i-1) , α ∼ 1 i-1+α i-1 j=1 δg j + α i-1+α δgnew , (1)
where g new is an integer not in g1:i, α > 0 is the concentration parameter, and δx is a discrete measure concentrated at x, which translates to a point mass of probability located at x. The discrete distribution (1) for the group assignment of individual or site i (conditional on the group assignments of individuals or sites 1, . . . , i -1) implies that each successive individual or site is assigned to an existing subgroup with probability proportional to the size of each subgroup, and is assigned to a new subgroup with probability proportional to α. The product of the successive conditional distributions given by (1) gives rise to the joint distribution of g1:N , which is the Chinese restaurant process (CRP) prior distribution with concentration parameter α [START_REF] Blackwell | Ferguson distributions via Pólya urn schemes[END_REF][START_REF] Pitman | Exchangeable and partially exchangeable random partitions[END_REF][START_REF] Pitman | Some developments of the Blackwell-MacQueen urn scheme[END_REF]. See [START_REF] Li | A tutorial on dirichlet process mixture modeling[END_REF] for more details and interpretations of the CRP prior distribution.

The strictly positive concentration parameter α of the CRP distribution influences the number of subgroups, through its control over the probability that individuals or sites are assigned into new subgroups. The larger the value of α, the more likely new subgroups are to be created. Figure 1 illustrates the effect of α on the number of subgroups created from the CRP prior distribution, when N = 60. For α = 0.1 generally only one or two subgroups are created, and infrequently three or more. When α = 0.5, the number of subgroups created generally falls between one and six. For α = 1, we seldom observe only one subgroup, and instead generally have between two and ten groups. In our non-parametric models, rather than fixing α to a specific value, we will use a hyperprior distribution for α to allow the degree of heterogeneity within the dataset itself to dictate the plausible range for α.

The possibility of the CRP prior using as many population subgroups as the total number of observations requires the inclusion of N distinct demographic parameters θ1, . . . , θN in our hierarchical non-parametric models. These parameters are assumed to be independent and identically distributed, just as in the finite mixture models.

Non-Parametric Capture-Recapture Model

In the non-parametric capture-recapture model, detection heterogeneity is flexibly addressed using a CRP prior distribution to allow for an unknown number of population subgroups. The number of population subgroups will thereby be inferred from the data, and is not fixed to a pre-specified value as in the finite mixture models described in section 2.2. We assign a CRP(α) prior distribution to g1:N , the vector of subgroup indicator variables. We use a Gamma(1, 1) hyperprior distribution for the CRP concentration parameter α.

We introduce probabilities of detection p1, p2, . . . , pN , each independent and identically following a Uniform(0, 1) prior distribution. In general, fewer than N subgroups are actually used by the CRP prior distribution, and thus only a subset of the probabilities of detection are "active" in terms of their influence on the model likelihood calculation. As in the finite mixture models, gi = k indicates that individual i is a member of subgroup k, and therefore has probability of detection pg i . The full non-parametric capture-recapture model is written as: φ ∼ Uniform(0, 1)

p k ∼ Uniform(0, 1) for k = 1, . . . , N α ∼ Gamma(1, 1) g1:N ∼ CRP(α) i = 1, . . . , N : xi,1 = 1 xi,t ∼ Bernoulli(φ • xi,t-1) for t = 2, . . . , T yi,t ∼ Bernoulli(pg i • xi,t) for t = 1, . . . , T

Non-Parametric Occupancy Model

Similarly, we generalize the finite mixture occupancy model to use a CRP(α) prior distribution for individual group assignments g1:N , and assign a Gamma(1, 1) hyperprior distribution for α. We include the maximum possible necessary number of distinct probabilities of detection p1, p2, . . . , pN , each with independent Uniform(0, 1) prior distributions. Here, gi = k indicates that site i is a member of subgroup k, and therefore has probability of detection pg i . The full non-parametric occupancy model is written as:

ψ ∼ Uniform(0, 1) p k ∼ Uniform(0, 1) for k = 1, . . . , N α ∼ Gamma(1, 1) g1:N ∼ CRP(α) i = 1, . . . , N : zi ∼ Bernoulli(ψ) yi,t ∼ Bernoulli(pg i • zi) for t = 1, . . . , T

Model Fitting via Markov chain Monte Carlo

Mathematical descriptions of our approaches to modeling heterogeneity within a population were given in section 2. These models are formulated in a Bayesian framework as hierarchical models. Prior distributions are specified for model parameters, and the conditional distributions for discrete latent states (e.g., dead/alive status) and the data likelihoods are given. The general tool for fitting hierarchical models to data is Markov chain Monte Carlo (MCMC; [START_REF] Brooks | Handbook of Markov Chain Monte Carlo[END_REF], a stochastic sampling algorithm for generating posterior samples for model parameters, conditional on the data. Demographic inferences are subsequently performed using the empirical posterior distributions of model parameters.

Markov chain Monte Carlo is a powerful tool for fitting hierarchical models to data, but the inherent mathematical complexity of MCMC sampling often requires the use of software. We make use of the recently developed nimble R package (NIMBLE Development Team, 2019), which offers new degrees of freedom for algorithmic development and customization of MCMC sampling [START_REF] De Valpine | Programming with models: writing statistical algorithms for general model structures with NIMBLE[END_REF]. The ability to write custom distributions for use in hierarchical models, and the flexible nature of nimble's MCMC engine has provided noteworthy gains in fitting ecological models [START_REF] Turek | Efficient markov chain monte carlo sampling for hierarchical hidden markov models[END_REF], and more generally in the study of MCMC algorithms [START_REF] Turek | Automated parameter blocking for efficient markov chain monte carlo sampling[END_REF].

We make use of nimble's MCMC engine for fitting these hierarchical models, after expressing them in the BUGS language [START_REF] Lunn | The BUGS project: Evolution, critique and future directions[END_REF]. In addition, we make use of two aspects of nimble's flexibility. For modeling of individual heterogeneity, we make use of the BNP distributions and corresponding sampling algorithms, which are a recent addition to the nimble package. Second, we enhance performance of the specific ecological models by using custom likelihood distributions provided by the nimbleEcology package [START_REF] Goldstein | nimbleEcology: Distributions for ecological models in nimble[END_REF] to remove latent states from the model structures.

Non-Parametric Distributions and MCMC Sampling

When a hierarchical model is formulated using the CRP distribution, the dCRP distribution (available in the nimble package) assigns the joint prior distribution arising from (1) to the labeling vector, g. Correspondingly, a specialized sampler is assigned by nimble's MCMC engine. Because the likelihood function of the ecological models presented in section 2 are not conjugate for the prior distribution of the cluster parameters, p k , the non-conjugate sampling algorithm described in Algorithm 8 of [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF] is assigned to g.

As described in section 2, under a non-parametric approach, the number of population subgroups is not fixed. In terms of the MCMC sampling scheme this means that the number of subgroups, and therefore the number of cluster parameters which are active, can vary with every MCMC iteration. As nimble does not support dynamic length allocation, the number of cluster parameters defined in the model must be fixed. A safe option would be to consider N cluster parameters, however this is highly inefficient both in terms of computation and storage, especially for large values of N . To reduce this inefficiency nimble allows the specification of N < N cluster parameters. If upon any MCMC iteration more than N groups are created, then a warning is issued.

Additionally, to reduce the computational burden of the nonparametric sampling, only the active cluster parameters, p k , are updated.

As discussed in section 2, the concentration parameter α has important implications in the clustering structure of the model. Therefore, efficient sampling of α is an important matter. Although its posterior distribution does not belong to any known class of distributions, when a gamma prior distribution is considered for α, a computationally efficient sampling scheme [START_REF] Escobar | Bayesian density estimation and inference using mixtures[END_REF], section 6) is assigned by nimble's MCMC engine.

Likelihood Distributions Using nimbleEcology

To reduce computation time, we make use of the nimbleEcology R package [START_REF] Goldstein | nimbleEcology: Distributions for ecological models in nimble[END_REF] for specifying the ecological hierarchical models. The nimbleEcology package provides likelihood distributions specific for a variety of common ecological models, which are implemented as custom distributions using the nimble package. The likelihood distributions provided in nimbleEcology include those for capture-recapture models, occupancy and dynamic occupancy models, and more generally for discrete hidden Markov models (HMMs) as appear in multi-state or multi-event capturerecapture (e.g., [START_REF] Turek | Efficient markov chain monte carlo sampling for hierarchical hidden markov models[END_REF].

For each type of ecological model, the likelihood distribution provided by nimbleEcology marginalizes over discrete latent states to directly calculate the unconditional likelihood of observed data. This allows removal of discrete latent state variables -xi,t alive/dead indicator variables in capturerecapture, and zi occupancy indicator variables in occupancy modelling -from the hierarchical model. This reduces model size and the necessary model computations, and increases the speed of MCMC mixing of top-level model parameters (φ in capture-recapture, and ψ in occupancy models) to generate stronger posterior inferences in less computational time.

Using the distributions provided in nimbleEcology to remove latent states does not alter the posterior results generated from each model; it only serves to increase the speed of generating inferences. The only noteworthy difference is that posterior inference for the discrete latent states cannot be performed, since samples for these latent states are never generated. So, for example, we could not perform inference for the alive/dead status xi,t of specific individuals in the capture-recapture setting.

Simulations

We undertake two simulation studies to assess performance of the various approaches to modeling individual and site heterogeneity. The first study is in the context of capturerecapture models, and the second study in that of occupancy models.

Both simulations consider the effect of varying degrees of heterogeneity in detection probability between two population subgroups on the accuracy of inferences. In each, we fit the homogeneous model, 2-group and 3-group finite mixture models, and the non-parametric model as were presented in section 2. R code for all simulations, including the nimble specifications for each model using the likelihood distributions provided in the nimbleEcology package, are provided as supplemental material.

Next we describe the details of each simulation, and then present results.

Capture-Recapture Simulation

In the capture-recapture simulation, we consider two population subgroups each with 800 individuals, for a total population size of 1,600. We use eight observational periods, and condition on all individuals being sighted on the first observational period. Survival probability is fixed at φ = 0.7 for simulating data, and we focus on the ability of various modeling approaches to estimate survival.

One population subgroup is fixed as having individual probability of detection p0 = 0.8, where detection is conditional on being alive. Individuals in the other subgroup have a fixed detection probability p, which varies between simulations. We consider values of p between 0.1 and 0.8 in increments of 0.1, where the terminal case p = 0.8 coincides with the detection probability of the first subgroup, and therefore represents a homogeneous population.

Occupancy Simulation

In the occupancy simulation, we consider two subgroups each with 2,000 sites, for a total of 4,000 sites. We use six independent observations of each site. The proportion of occupied sites is fixed at ψ = 0.7 for simulating data, and we focus on the ability of various modeling approaches to estimate the true occupancy proportion.

One subgroup is fixed as having probability of detection p0 = 0.8 on each independent site visit, where detection is conditional on a site being occupied. The other subgroup has a fixed detection probability p, which varies between simulations. We consider values of p between 0.1 and 0.8 in increments of 0.1, but with a finer resolution on the interval between 0.1 and 0.4. Again, the terminal case p = 0.8 coincides with the detection probability of the first subgroup, and therefore represents a homogeneous population.

Simulation Results

Posterior inferences were performed for individual survival probability φ in the capture-recapture simulation and for site occupancy proportion ψ in the occupancy simulation. We use the posterior median and equal-tailed 95% Bayesian credible intervals for inferences under each model. Fitted models include a homogeneous model (Hom), 2-group and 3group finite mixture models (FM 2 and FM 3, respectively), and a non-parametric model (NP). Each model was fit to simulated data using MCMC, as described in section 3.

Results for the capture-recapture simulation appear in Fig- ure 2. The homogeneous model consistently under-estimates φ, more severely for larger discrepancies in detection probability between the two groups, with the discrepancy diminishing and disappearing as the detection probabilities of the two groups converge. For the lowest value of p = 0.1, and hence a large difference in detection probability between groups, the 3-group mixture model has a slight tendency to inflate estimates of φ, while the opposite is true of the 2-group mixture model; however this difference is minor and quickly disappears as the group difference deceases. Otherwise, the mixture models and non-parametric models are all similarly successful in generating accurate inferences of φ. We also note that all models exhibit a slight negative bias in their estimates of φ, which diminishes as p0 and p jointly approach one. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.60 Results for the occupancy model simulation appear in Figure 3. Once again, for low values of p the homogeneous model underestimates ψ, more severely for larger discrepancies between the two groups. We also see a regime of p values between (approximately) 0.167 and 0.267, in which the 3group mixture model vastly over-estimates ψ, with posterior median estimates in excess of 0.99. The precise location of this regime varied somewhat depending on simulation parameters (specifically the number of individuals in the population, and number of observation periods) but its existence persisted in all simulations. That is, there exists the potential for highly inaccurate inferences when the number of groups in a finite mixture model is chosen either too high, or too low.

The 2-group mixture model (which is also the datagenerating model) generated reliable inferences for all values of p. The non-parametric model admitted greater uncertainty for the upper limit of the 95% credible interval, especially for lower values of p, but the posterior median estimates generated using the non-parametric model were consistently near to the true parameter value ψ = 0.7.

Examples

We present two real-data examples, one in the context of capture-recapture and the other in that of occupancy modeling. For each example, we fit a homogeneous model for detection probability, finite mixture models containing between two q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.5 and ten population subgroups, and a non-parametric model.

In addition to parameter inferences under each model, we also present the WAIC [START_REF] Watanabe | Asymptotic equivalence of bayes cross validation and widely applicable information criterion in singular learning theory[END_REF][START_REF] Gelman | Understanding predictive information criteria for bayesian models[END_REF] of each fitted model. The WAIC value is a measure of the goodness-of-fit of a hierarchical model, calculated using chains of posterior MCMC samples. WAIC is measured on the scale of deviance, and therefore lower values of WAIC indicate a more parsimonious fit to the data. Finally, for each example we also present the posterior distribution for the number of population subgroups in the non-parametric model. This suggests at the degree of heterogeneity present in the data as inferred using the non-parametric model, the only model which does not pre-define the number of population subgroups. The datasets used for each example are available on GitHub at: https://github.com/danielturek/bnp-examples-data.

Capture-Recapture Example

We consider wolf (canis lupus) capture-recapture data collected in France between 1995 and 2003, as studied in [START_REF] Cubaynes | Importance of accounting for detection heterogeneity when estimating abundance: the case of French wolves[END_REF]. The original data contains binary detection data for a total of 87 wolves, over T = 8 observation periods. Since we condition on the first sighting of each individual, we excluded the twenty individuals who were first observed on the final observation period. This leaves a total of N = 67 unique individuals in the dataset. Models for detection heterogeneity described in section 2 were fit to this data using MCMC, and inference was performed for individual survival probability φ under each model.

Posterior median and 95% credible intervals, as well as the WAIC value of each fitted model are presented in Table 1. The homogeneous model produces the largest WAIC value among those models considered (indicating the poorest fit to the We note that the non-parametric model produces the third lowest WAIC value, being 1.5 higher than the lowest two WAIC values (equal from the 4-group and 5-group mixture models). Further, inferences for φ generated under the nonparametric model are nearly identical to those of the 4-group and 5-group mixture models. Given our uncertainty in the structure and degree of heterogeneity, the non-parametric model provides defensible inferences and goodness-of-fit. Figure 4 displays the posterior distribution for the number of population subgroups, as inferred by the non-parametric model. We observe a right-skewed distribution placing the bulk of the posterior mass roughly between three and ten subgroups, with a posterior median of seven, the mode appearing at six subgroups, and a 90% Bayesian credible interval for the number of subgroups as (3,13). This distribution suggests that our consideration of finite mixture models containing between two and ten subgroups was a reasonable choice, even though this was an uninformed and somewhat arbitrary selection.

Occupancy Example

For our occupancy example, we consider a second wolf (canis lupus) dataset [START_REF] Louvrier | Accounting for misidentification and heterogeneity in occupancy studies using hidden Markov models[END_REF] collected in France in 2013. Opportunistic observational data such as tracks, scat, and prey remains were collected from N = 3,211 grid cells, each being a 10km × 10km square. Each site was surveyed on a total of T = 4 independent observation occasions. The categorical data, in total, consisted of 250 "unambiguous detections", 54 "ambiguous detections", and 12,540 "nondetections". We convert this to binary data, wherein both "unambiguous" and "ambiguous" detections are considered to be positive detections. Models for detection heterogeneity described in section 2 were fit to this data using MCMC, and inference was performed for site occupancy proportion ψ under each model. Posterior median and 95% credible intervals, as well as the WAIC value of each fitted model are presented in Table 2. The homogeneous model again produces the largest WAIC value (indicating the poorest fit to the data), and the lowest estimates of ψ. The 2-group mixture model yields the secondhighest WAIC value and also lower estimates of ψ than the remaining models, and the 3-group mixture model yields the third highest WAIC value. This suggests a nontrivial degree of heterogeneity in detection probability between sites.

The remaining mixture models (K 4 groups) give WAIC values between 2168.7 and 2175.9, and exhibit small variations in the inferences for ψ, in particular in the upper limit of the 95% credible interval. In contrast, the non-parametric model yields the lowest WAIC value among all models (2166.7) indicating the best fit to the data. The posterior median estimate from the non-parametric model is similar to that of the K 4 mixture models, but the credible interval is wider. In particular, the non-parametric model suggests a higher upper-bound for the 95% credible interval for ψ than any other model considered. This result may be reasonable, since WAIC suggests the non-parametric model provides the best fit to the data, and we are completely uncertain as to the degrees of detection heterogeneity which exists in this opportunistic sampling dataset.

We again present the posterior distribution for the number of population subgroups as inferred by the non-parametric model in Figure 5. This distribution is again right-skewed, this time spanning a wider range for plausible numbers of subgroups. This makes intuitive sense, as our occupancy example data considers a much larger number of observations (3,211 sites, rather than a total of 67 individuals in the capturerecapture example). The posterior distribution has a median of eight, once again the mode appears at six subgroups, and has a 90% Bayesian credible interval of (3, 19).

Discussion

Here, we have employed a Bayesian non-parametric approach to modeling heterogeneity in the detection process of ecological models. Our approach used a Chinese restaurant process (CRP) prior distribution to model group memberships, which has the benefit of not requiring specification of the number of population sub-groups a priori. Using the CRP prior, both the distribution of individuals or sites among distinct subgroups and the number of subgroups are inferred from the data. This strategy obviates any process of model selection used to choose a "best value" for the number of groups K when using a finite mixture model.

The tendency of the CRP distribution towards using fewer (and hence more diverse) subgroups, or towards using a larger number of tightly specified subgroups is governed by a concentration parameter (α). The effect of different values of α is seen in Figure 1, although we avoid the need to choose a particular value of α by specifying an uninformative Gamma prior distribution. This approach of using a CRP prior distribution to model heterogeneity is used in diverse areas of research including topic modeling [START_REF] Blei | The nested chinese restaurant process and bayesian nonparametric inference of topic hierarchies[END_REF], genomics [START_REF] Qin | Clustering microarray gene expression data using weighted chinese restaurant process[END_REF], and evolutionary clustering [START_REF] Ahmed | Dynamic non-parametric mixture models and the recurrent chinese restaurant process: with applications to evolutionary clustering[END_REF], among others.

The examples and simulations used herein focus on the most basic forms of capture-recapture and occupancy models. Specifically, we have considered the Cormack-Jolly-Seber capture-recapture model, and a single-season static occupancy model. Although our examples focused on the basic forms of each model, the same technique are readily applied to more complex variations of each model. For example, we could apply the same non-parametric CRP prior distribution in multi-state capture-recapture, or in dynamic, multi-season, or multi-species occupancy models. Similarly, one can readily extend the hierarchical models to incorporate individual or environmental covariates affecting population demographics, survey effort, or the detection processes. The approach we have used is quite general, and widely applicable to ecological models where population heterogeneity is a consideration.

We have used the nimble R package to specify the hierarchical models described herein, and to fit these models to data using MCMC. nimble provides degrees of flexibility which are not available in other software packages. Specifically, we leverage nimble's ability to use custom-written likelihood distributions in a hierarchical model, specifically the ecological likelihood distributions provided in the nimbleEcology package. Further, nimble provides the ability to specify the sampling algorithms used by the MCMC, and even to oneself write customized sampling algorithms for use in the MCMC. Indeed, the MCMC sampling algorithms used for fitting our models -specifically those used for the CRP concentration parameter α, and for the CRP-distributed group membership indicators g1:N -are themselves custom sampling algorithms written for precisely these non-parametric motifs, and added into nimble's MCMC repertoire of algorithms. That said, nimble does not attempt to provide "canned" algorithms, nor any particular pre-written model structures, but rather an environment for writing custom functions, statistical algorithms, and distributions, and the application of these to generallyspecified hierarchical model structures. The goal of nimble is to provide a flexible model and algorithmic programming environment to facilitate highly efficient analysis of models and complex data.

It is common that heterogeneity will be present to some degree in the detection process of ecological models. In practice, this may be detected by goodness-of-fit tests [START_REF] Jeyam | Assessing heterogeneity in transition propensity in multistate capture-recapture data[END_REF], or perhaps based on prior expert knowledge. When detection heterogeneity is known or suspected to be present, and suitable covariates are not available to accurately model this heterogeneity, we recommend using a BNP modelling approach. This approach alleviates the necessity of selecting the number of components used in a finite mixture model, which is an inherently difficult and oftentimes subjective process. No less, the exact number of mixture components is generally not the primary inferential focus. Use of a BNP modelling approach, as demonstrated herein, accounts for whatever degree of heterogeneity may be present while requiring no subjective choices or guesswork. This provides an effective approach to reducing bias in the resulting demographic inferences.

Figure 1 .

 1 Figure 1. Prior probability of number of subgroups induced by the CRP(α) prior distribution for different values of α, when N = 60. A color version is available in the electronic version of this article.

Figure 2 .

 2 Figure 2. Capture-recapture simulation results for survival probability φ using different models for heterogeneity in detection probability. N = 1, 600 individuals are simulated between two population subgroups, one with probability of detection p0 = 0.8, and the other with probability of detection p. True survival probability is fixed at φ = 0.7. Solid lines show posterior median estimates of φ from each model, and dashed lines show upper and lower limits of a 95% Bayesian credible interval. A color version is available in the electronic version of this article.

Figure 3 .

 3 Figure 3. Occupancy simulation results for occupancy ψ using different models for heterogeneity in detection probability. N = 4, 000 individuals are simulated between two population subgroups, one with probability of detection p0 = 0.8, and the other with probability of detection p. True occupancy is fixed at ψ = 0.7. Solid lines show posterior median estimates of ψ from each model, and dashed lines show upper and lower limits of a 95% Bayesian credible interval. A color version is available in the electronic version of this article.

Figure 4 .

 4 Figure 4. Posterior distribution of the number of population subgroups under the non-parametric model, for the capture-recapture example.

Figure 5 .

 5 Figure 5. Posterior distribution of the number of population subgroups under the non-parametric model, for the occupancy model example.

Table 1

 1 Capture-recapture example results. Posterior inferences are for individual survival probability φ, and WAIC values indicate the goodness-of-fit of each model.

	Model	Median 95% BCI WAIC
	Homogeneous	.80	(.70, .89)	237.5
	Finite Mixture 2	.91	(.79, .99)	209.1
	Finite Mixture 3	.91	(.80, .99)	201.9
	Finite Mixture 4	.92	(.80, .99)	198.3
	Finite Mixture 5	.92	(.80, .99)	198.3
	Finite Mixture 6	.91	(.80, .99)	199.9
	Finite Mixture 7	.91	(.79, .99)	200.6
	Finite Mixture 8	.90	(.79, .99)	201.4
	Finite Mixture 9	.90	(.79, .98)	202.1
	Finite Mixture 10	.90	(.78, .98)	202.8
	Non-Parametric	.92	(.81, .99)	199.8
	data), and the lowest estimates of φ. Posterior inferences from
	all mixture models and the non-parametric model are nearly
	indistinguishable, with median posterior values for φ around
	0.91, and 95% credible intervals of approximately (0.80, 0.99).

Table 2

 2 Occupancy model example results. Posterior inferences are for site occupancy proportion ψ, and WAIC values indicate the goodness-of-fit of each model.

	Model	Median	95% BCI	WAIC
	Homogeneous	.063	(.054, .073) 2237.0
	Finite Mixture 2	.079	(.063, .101) 2195.8
	Finite Mixture 3	.087	(.067, .122) 2178.3
	Finite Mixture 4	.092	(.068, .139) 2168.7
	Finite Mixture 5	.089	(.068, .132) 2172.4
	Finite Mixture 6	.090	(.069, .130) 2172.1
	Finite Mixture 7	.089	(.068, .125) 2173.7
	Finite Mixture 8	.087	(.068, .122) 2174.9
	Finite Mixture 9	.087	(.068, .119) 2175.9
	Finite Mixture 10	.086	(.069, .119) 2175.4
	Non-Parametric	.090	(.066, .359) 2166.7

Acknowledgements

Support for DT was provided by Fulbright Research Scholarship Award 9183-FR and also the Williams College Class of 1945 World Fellowship. Support for CW was partially provided by award nsf-dms 1622444. Support for OG was provided by a grant from CNRS and "Mission pour l'interdisciplinarité" through its "Osezl'interdisciplinarité call." We warmly thank the French Office of Biodiversity (OFB) for sharing the wolf datasets.

Appendix Simulation Code library(nimble) library(nimbleEcology)

## number of individuals N <-800 ## number of observation periods T <-8 ## time period of first capture first <-rep(1, N) ## length of observation history from first capture len <-T -first + 1 ## survival probability phi <-0.7 ## detection probability pVec <-rep(c(0.2, 0.8), each=N/2) ## simulate z (alive / dead status), ## and y (encounter histories) set.seed(0) z <-matrix(NA, nrow=N, ncol=T) y <-matrix(NA, nrow=N, ncol=T) for(i in 1

## Homogeneous Capture-Recapture Model code <-nimbleCode({ phi ~dunif(0, 1) p ~dunif(0, 1) for(i in 1:N) { y [i,first[i]:T] ~dCJS_ss(phi, p, len=len[i]) } }) constants <-list(N=N, T=T, first=first, len=len) data <-list(y=y) inits <-list(phi=0.5, p=0.5) Rmodel <-nimbleModel(code, constants, data, inits)

## fixed number of groups constants <-list(N=N, T=T, first=first, len=len, K=K) data <-list(y=y, one=rep(1,K-1)) inits <-list(phi=0.5, pi=rep(1/K,K), p=rep(0.5,K), g=rep(1,N)) Rmodel <-nimbleModel(code, constants, data, inits)

## maximum number of subgroups constants <-list(N=N, T=T, first=first, len=len, M=M) data <-list(y=y) inits <-list(phi=0.5, alpha=1, xi=rep(1,N), p=rep(0.5,M)) Rmodel <-nimbleModel(code, constants, data, inits)

## number of sites N <-4000 ## number of observation periods T <-6 ## probability of occupancy pOcc <-0.7 ## detection probability pVec <-rep(c(0.2, 0.8), each=N/2) ## simulate z (occupied status), ## and y (encounter histories) set.seed(0) z <-rep(NA, N) y <-matrix(NA, nrow=N, ncol=T) for(i in 1

pOcc ~dunif(0, 1) p ~dunif(0, 1) for(i in 1:N) { y[i,1:T] ~dOcc_s(pOcc, p, len=T) } }) constants <-list(N=N, T=T) data <-list(y=y) inits <-list(pOcc=0.5, p=0.5) Rmodel <-nimbleModel(code, constants, data, inits)

## fixed number of groups constants <-list(N=N, T=T, K=K) data <-list(y=y, one=rep(1,K-1)) inits <-list(pOcc=0.5, pi=rep(1/K,K), p=rep(0.5,K), g=rep(1,N)) Rmodel <-nimbleModel(code, constants, data, inits) ## 3-Group Finite Mixture Occupancy Model code <-nimbleCode({ pOcc ~dunif(0, 1) for(k in 1:K

## fixed number of groups constants <-list(N=N, T=T, K=K) data <-list(y=y, one=rep(1,K-1)) inits <-list(pOcc=0.5, pi=rep(1/K,K), p=rep(0.5,K), g=rep(1,N)) Rmodel <-nimbleModel(code, constants, data, inits)

## maximum number of subgroups constants <-list(N=N, T=T, M=M) data <-list(y=y) inits <-list(pOcc=0.5, alpha=1, xi=rep(1,N), p=rep(0.5,M)) Rmodel <-nimbleModel(code, constants, data, inits)

## compile model and MCMC Cmodel <-compileNimble(Rmodel) Cmcmc <-compileNimble(Rmcmc, project=Rmodel) set.seed(0) samplesList <-runMCMC (Cmcmc, niter=10000, nchains=3)