Oksana Grentea

Thibault Saubussea

Olivier Gimenezb

Eric Marboutinc

Christophe Duchampa

Wolf depredation hotspots in France: Clustering analyses adjusting for livestock availability

Keywords: Human-carnivore conflict, depredation hotspot, conservation resource allocation, gray wolf, clustering analysis

Areas exhibiting high levels of predations on livestock generate conflicts between humans and large carnivores. Managers generally seek to identify these hotspots, in order to diagnose the causes that lead to hotspot formations and to provide financial or technical support to the involved livestock owners. When locating depredation hotspots, previous studies have not adjusted for livestock availability, making it difficult for managers to discriminate hotspots resulting from underlying livestock clusters from those due to other factors such as environmental factors. We studied hotspots of wolf depredation on sheep in France from the beginning of the natural wolf recolonisation in 1994 up to 2018. For each year, we applied the Ripley's K-function and Ripley's Kinhom to determine the general depredation spatial pattern and the Kulldorff statistic to locate depredation hotspots. We showed that omitting livestock availability in these analyses led to flawed inference about the depredation pattern, and resulted in a substantial number of unidentified hotspots, including pastoral surfaces with low sheep availability. Our methodology provides reliable information for managers to understand the depredation pattern over space and time and to allocate resources.

Introduction

Predation on domestic prey by large carnivores, hereafter depredation, is the main driver of conflicts between humans and large carnivores [START_REF] Lute | Conservation professionals agree on challenges to coexisting with large carnivores but not on solutions[END_REF]. The financial and social costs associated to depredations reduce acceptance of these species and lead to retaliatory killings [START_REF] Abade | Assessing the relative importance of landscape and husbandry factors in determining large carnivore depredation risk in Tanzania's Ruaha landscape[END_REF]. Understanding the spatial and temporal patterns of depredation is a major challenge in large carnivore conservation to mitigate current conflicts and prevent future ones [START_REF] Miller | Mapping attack hotspots to mitigate human-carnivore conflict: approaches and applications of spatial predation risk modeling[END_REF]. To do so, two complementary approaches can be applied, namely the spatial correlation analysis and risk modelling. The spatial correlation analysis quantifies the level of clustering or regularity (i.e. repulsion) of a depredation pattern [START_REF] Baddeley | Spatial Point Patterns: Methodology and Applications with R[END_REF][START_REF] Hoffmann | Spatial Pattern Analysis Reveals Randomness Among Carnivore Depredation of Livestock[END_REF] and enables the identification of hotspots or coldspots where events are unusually aggregating or scarce, respectively. Risk modelling aims at estimating the spatial risk of depredation, by quantifying the relationship between ecological or nonecological features and depredation numbers or occurrences [START_REF] Miller | Mapping attack hotspots to mitigate human-carnivore conflict: approaches and applications of spatial predation risk modeling[END_REF].

The spatial correlation analysis is often presented as a preliminary step to risk modelling and as of limited interest if applied alone [START_REF] Gastineau | Spatiotemporal depredation hotspots of brown bears, Ursus arctos, on livestock in the Pyrenees, France[END_REF][START_REF] Hoffmann | Spatial Pattern Analysis Reveals Randomness Among Carnivore Depredation of Livestock[END_REF]. This may explain why risk modelling is more popular (e.g. [START_REF] Bradley | Assessing factors related to wolf depredation of cattle in fenced pastures in Montana and Idaho[END_REF][START_REF] Fowler | Non-linear relationships between human activities and wolf-livestock depredations[END_REF] than spatial correlation analyses in the scientific literature. However, risk models only reflect the depredation-factor relationships at a given point in time and space. These relationships may evolve because predator-prey systems are dynamic and carnivores can adapt their predation behaviour according to environmental or livestock management changes [START_REF] Miller | Mapping attack hotspots to mitigate human-carnivore conflict: approaches and applications of spatial predation risk modeling[END_REF]. This is especially true for (re)colonising large carnivore populations of plastic species such as the grey wolf (Canis lupus) which can embrace a large range of habitats [START_REF] Llaneza | Insights into wolf presence in humandominated landscapes: the relative role of food availability, humans and landscape attributes[END_REF][START_REF] Mech | Wolves: Behavior, Ecology, and Conservation[END_REF]. Risk models developed for this species were either restricted to a unique region and type of habitat (e.g. [START_REF] Clark | Predicting Spatial Risk of Wolf-Cattle Encounters and Depredation[END_REF] or showed a low predictive power when applied to habitats that differ from those originally used in the analysis [START_REF] Hanley | Forecasting cattle depredation risk by recolonizing gray wolves[END_REF]. Therefore, the outcomes of risk models are hardly transposable to other areas or on long-term periods, which makes their adoption by managers difficult [START_REF] Miller | Mapping attack hotspots to mitigate human-carnivore conflict: approaches and applications of spatial predation risk modeling[END_REF]. In contrast, information about current or past hotspot locations provided by spatial correlation analyses can directly help managers when allocating conservation resources such as subsidies of preventive measures or when investigating the causes of depredation in situ. Prioritising conservation efforts on depredation hotspots is indeed recommended as recurrences of high levels of attacks in the same areas are likely to trigger negative attitudes [START_REF] Stahl | Predation on livestock by an expanding reintroduced lynx population: long-term trend and spatial variability[END_REF].

So far, several analyses of spatial correlation have been developed to identify depredation hotspots. The simplest one consists in defining an arbitrary threshold for a selected indicator like the number of depredations which, if reached, turns the spatial unit into a hotspot. Units can be administrative areas (e.g. [START_REF] Dhungana | Livestock depredation by leopards around Chitwan National Park, Nepal[END_REF] or simple shapes like circles [START_REF] Stahl | Predation on livestock by an expanding reintroduced lynx population: long-term trend and spatial variability[END_REF]. To avoid arbitrary thresholds, some studies have applied spatial statistical analyses (e.g. [START_REF] Gastineau | Spatiotemporal depredation hotspots of brown bears, Ursus arctos, on livestock in the Pyrenees, France[END_REF][START_REF] Hoffmann | Spatial Pattern Analysis Reveals Randomness Among Carnivore Depredation of Livestock[END_REF][START_REF] Packer | Speciesspecific spatiotemporal patterns of leopard, lion and tiger attacks on humans[END_REF] where depredation events are considered as a spatial point pattern, i.e. a dataset of observed spatial locations of a biological process [START_REF] Baddeley | Spatial Point Patterns: Methodology and Applications with R[END_REF]. The areas exhibiting significantly more depredation events than expected under Complete Spatial Randomness (CSR) are statistically identified as hotspots. To gain understanding of the biological process, the spatial statistical analyses can also be conducted to determine the summarised spatial structure of the pattern, i.e. the range of distances over which the pattern generally exhibits clustering, randomness or regularity [START_REF] Baddeley | Spatial Point Patterns: Methodology and Applications with R[END_REF], as in [START_REF] Kushnir | Using landscape characteristics to predict risk of lion attacks on humans in south-eastern Tanzania[END_REF].

While studies on depredation risk modelling generally controlled for spatio-temporal livestock availability (e.g. [START_REF] Balbuena-Serrano | Hotspots of livestock depredation by pumas and jaguars in Brazil: a biome-scale analysis[END_REF][START_REF] Fowler | Non-linear relationships between human activities and wolf-livestock depredations[END_REF], quantitative studies of spatial correlation did not. The reason was that the statistical method they used did not allow it or because the relevant data were not available. Yet the spatial pattern of depredation is inherent to livestock risk exposure. Not only are conflicts exclusive to areas where carnivore and livestock distributions overlap, but they are also dependent on livestock availability, which is itself a by-product of livestock counts and time spent in a specific area.

If livestock availability is ignored, it is difficult to disentangle the hotspots based on high levels of livestock availability from the hotspots where livestock densities are low but predation rates disproportionally high. In the latter cases, conflicts may be favoured by factors linked to, e.g., the environment, the predator or the husbandry practices which can be identified and managed to some extent.

The grey wolf has been naturally recolonising France since the early 1990s from the Italian Apennines. Since then, wolves have been expanding and densifying in the south-east part of France [START_REF] Louvrier | Mapping and explaining wolf recolonization in France using dynamic occupancy models and opportunistic data[END_REF]. The recolonisation has come along with a gradually increasing number of depredations, with more than 3,000 attacks on livestock in 2018, 90% of wolf attacks being on sheep (Ovis aries). French authorities have so far focused on hotspots to manage the conflict by subsidising preventive measures according to the level and recurrence of depredation events at the town level (Ministère de l'Agriculture et de l'Alimentation, 2019), but without accounting for sheep availability at risk to identify hotspots.

Here, we aimed at conducting spatial statistical correlation analyses of wolf predations on sheep in south-east France, while accounting for sheep risk exposure. First, we determined the summarised spatial structure of the depredation pattern and second, we located the significant depredation hotspots. We compared how the adjustment for sheep availability in space and time affected the results. We conducted annual analyses, from 1994 to 2018, to study the temporal variations in depredation patterns and in hotspot locations throughout wolf recolonisation. We also explored the depredation patterns at a local scale, within the historic area of wolf recolonisation. The aims of the local scale analysis were to increase perceptibility of wolf territorial behaviour in the results, and to observe if years of wolf presence changed the depredation pattern over time compared to the regional scale where the colonisation process was still ongoing. Finally, we discuss the management implications of our results.

Materials and methods

Study area and study period

The study area covered the two south-east regions of France, Provence-Alpes-Côte-d'Azur We analysed the summarised spatial structure of the depredation pattern at two scales of the study area. First, the analyses were applied over the whole study area, i.e. the regional scale.

Second, the analyses were applied over a sub-area of the study area, located in the very southeast of the study area (2 146 In accordance with wolf biology, we defined the annual periods over which the analyses were applied as biological years y, starting from the 1 st of April of year y-1 (i.e. wolf pup birth) to the 31 st March of year y. The whole study period covered biological years (hereafter, 'years') from 1995 to 2018.

Depredation records

Approximately three quarters of the depredations occurred between June and October, when most of sheep flocks were grazing in pastures. Up to 2019, the French Ministry of Ecology compensated for any killed farmed animal for which wolf responsibility could not be discarded, regardless of preventive measures. Each claim was controlled and checked in the field by an accredited governmental agent who used a standardized protocol [START_REF] Duchamp | A dual frame survey to assess time-and space-related changes of the colonizing wolf population in France[END_REF]. Therefore, most if not all depredations were reported if noticed, and the risk of false claims was low. We restricted the clustering analyses to depredations on sheep only as they constituted 90% of depredations on livestock. An event of depredation corresponded to at least one killed or wounded sheep by wolves. Because grazing activity may extend outside the official pastoral limits, or because of geolocation approximations, some depredation events were not located inside the geolocated pastoral surfaces (22% of the annual dataset on average). We excluded the depredations farther away from 500 meters of any pastoral surface (10% of the annual dataset on average) and assigned the remainder to their nearest pastoral surface.

Wolf distribution

The Wolf-Lynx French Network managed by the French Biodiversity Agency is in charge of the national wolf monitoring in France. Its field experts are trained to opportunistically collect presence signs all year along, such as biological samples (mainly faeces, combined with genetic analysis), tracks or direct observations, which are geolocated and validated according to a standardized protocol [START_REF] Duchamp | A dual frame survey to assess time-and space-related changes of the colonizing wolf population in France[END_REF][START_REF] Louvrier | Mapping and explaining wolf recolonization in France using dynamic occupancy models and opportunistic data[END_REF]. The annual wolf distribution was determined by combining the annual presence signs and all wolf depredations (i.e. on any type of livestock, including all those located outside the pastures) within the study area. Wolf presence during year y was reported on a 10x10 km cell grid, and defined as the collection of at least one presence sign or one depredation in the cell from y-1 to y-3 and from y-2 to y (see [START_REF] Marboutin | On the effects of grid size and shape when mapping the distribution range of a recolonising wolf (Canis lupus) population[END_REF]. The use of depredation data to estimate the wolf distribution was necessary to identify the pastures that were at the edges of the wolf recolonisation area. Indeed, depredations (along with visual observations) are generally the first detections of wolves that settle in newly recolonized areas in France [START_REF] Duchamp | A dual frame survey to assess time-and space-related changes of the colonizing wolf population in France[END_REF].

Spatio-temporal sheep availability

Around 1 500 000 sheep are bred for meat or milk production each year in the study area [START_REF] Idele | Ovins 2017 Productions lait et viande, Les chiffres clés du GEB[END_REF]. During the summer period, a large part of these flocks becomes transhumant and moves to high-altitude pastures in the Alps. Otherwise, sheep are grazing in low-or midaltitude pastures around farms. Sheep are in sheds during the whole winter except in the south where climate is mild (Gervasi et al., submitted).

Two georeferenced censuses of the pastoral surfaces in the study area were carried out in 1996-1997 and in 2012-2014 by the National Research Institute of Science and Technology for the Environment and Agriculture (IRSTEA). For each pastoral surface, information about the grazing livestock count per species and the annual number of grazing days was included.

One or several herds, of one or multiple livestock species, could feed on a pastoral surface.

We restricted the dataset to pastoral surfaces with grazing sheep, regardless of whether other species of livestock were present (Fig. between the north and the south of the regional scale due to climate differences and resulting grass phenology. The pastoral surfaces of the first census and of the second census were respectively 2.7 km² and 1.9 km² on average, and they hosted 606 sheep (sd=649) and 516 sheep (sd=547), during 146 days (sd=96) and 202 days (sd=96) per year on average.

For each year y, we identified the set of pastoral surfaces at depredation risk for both scales by selecting the pastoral surfaces which overlapped, even partially, with the wolf distribution.

For years between 1995 to 2005, we used the 6 488 pastoral surfaces with sheep from the oldest census, and for years between 2006 to 2018, we used the 12 438 pastoral surfaces with sheep from the latest census. The oldest census did not include pastures in the west of the study area (i.e. outside the Alpine area), which explains the lower number of pastures compared to the latest census. However, the wolf distribution from 1995 to 2005 was mainly restricted to the east of the study area, within the Alpine area (Fig. B.2), where pastoral information in the oldest census was available. For example, using the latest census instead of the oldest census only provided a difference of 6.7% on average in the numbers of pastoral surfaces at depredation risk at the regional scale from 1995 to 2005. Thus, we considered this gap was of little impact.

IRSTEA did not carry out the census process every year for two reasons. First, the census process was the result of a long period of work, exceeding a year of work. The field work, which consisted in meetings with local contact persons, corresponded to 750 full working days on its own for the second census (F. Bray, pers. com.). Second, the use of pastoral surfaces by farmers was relatively stable over time, because the environmental characteristics (area, topography, type of vegetation, weather…) of a pasture restricted the herd size the pasture could host and the time of grazing. Moreover, the large but limited number of pastures caused hard competition between livestock owners for access to pastoral resources. Livestock owners tended to secure their access by renting or buying these surfaces for a long period of time (C. Duchamp, pers. com.). For all these reasons, we assumed the information on spatiotemporal sheep distribution of a census was valid for several years before and after the census.

Analysis 1: Summarised spatial structure of wolf depredations

For each year y, we analysed at the regional and local scales the summarised spatial structure of the wolf depredation pattern considering sheep availability by using Ripley's Kinhom [START_REF] Baddeley | Spatial Point Patterns: Methodology and Applications with R[END_REF]. We used the package spatstat [START_REF] Baddeley | spatstat: An R Package for Analyzing Spatial Point Patterns[END_REF] in R 4.0.2 (R Core Team, 2019) to perform the analysis.

First, we simulated the expected depredation pattern under CSR considering sheep availability. The simulated pattern was composed of n points called 'controls', in opposition to the n annual observed depredations called 'cases'. The distribution of controls was simulated on a grid of 200x200 m cells over the pastoral surfaces at depredation risk, according to an inhomogeneous Poisson Point Process (PPP). In this model, the number of controls falling into a cell g was expected to be equal to the intensity !(#) defined as:

!(#) = exp(a + b × log (121 ! ) + c × log (4567 ! ))
where 121 ! and 4567 ! were the numbers of sheep and of grazing days of the pastoral surface u covered by g, and a, b and c were model parameters to be estimated. Covariates were logtransformed because of their skewed distributions. Thus, the controls were more likely to fall into the cells covering pastoral surfaces with high numbers of 121 ! and 4567 ! . For each year, we simulated 499 control patterns through this model.

Second, we computed Kinhom for the pattern of n cases. This function drew a circle of radius r around each depredation i covering more or less neighbouring depredations {j}. Then, the function summed the values 1 (!(# " )!9# # :) ⁄ for all pairs of depredations i-j within this circle, !(# " ) and !(# # ) being the intensity values of the cell(s) containing i and j respectively.

Because the unit of measurement of ! was in ind/cell (i.e. ind/0.04km²), Kinhom returned a surface for each tested r. If the pattern was randomly distributed considering sheep availability, the surface should equal the surface of the circle, <= $ . Consequently, at a specific r, if the Kinhom value of the observed depredation pattern was higher or lower than <= $ , we could conclude that the depredations were respectively forming hotspots or coldspots of radius r. A large set of continuous radiuses r was tested.

However, the smaller the dataset, the larger the likelihood of differences from the expected number of points within a circle, even if the dataset has a random pattern. Therefore, a Kinhom value which differed from <= $ was not necessarily the result of a non-randomly point pattern, but could be due to stochastic effects. To take this into account, we also computed Kinhom for the 499 control patterns and selected the 50 th -lowest and 50 th -highest values to build a 'control envelope'. For a specific r, if the observed Kinhom was higher (or lower) than the upper (or lower) envelope limit, the depredations were significantly clustered (or regular, i.e. tended to avoid each other) at this distance considering sheep availability. The larger the differences with the envelope, the more aggregated or regular the pattern.

We also annually applied the Ripley's K-function at the regional scale. The PPP of the Kfunction was homogeneous. Hence, the number of controls falling into the cells was expected to be equal over the whole analysis area [START_REF] Baddeley | Spatial Point Patterns: Methodology and Applications with R[END_REF]. We tested two possibilities to define the grid of 200x200 m cells. First, we simulated the same grid over the pastoral surfaces at depredation risk as in the Kinhom analysis, in order to simulate a situation for which sheep distribution was known, but sheep availability was not. Second, we simulated the grid over the wolf distribution, to simulate a situation for which pastoral data were totally unavailable.

Analysis 2: Wolf depredation hotspot location

For each year y, we located at the regional scale the significant wolf depredation hotspots by applying the first version of the Kulldorff statistic [START_REF] Kulldorff | A spatial scan statistic[END_REF] which adjusts for spatially inhomogeneous population at risk. We used the package SpatialEpi [START_REF] Kim | SpatialEpi: Methods and Data for Spatial Epidemiology[END_REF] in R 4.0.2 (R Core Team, 2019) to perform the analysis.

First, we defined all the possible zones Z which could be later identified as hotspot by the Kulldorff statistic. We defined the pastoral surface as the smallest possible zone Z. Then, a circle was put on the centroid of each pastoral surface, with a progressively increasing radius.

Each time the growing circles encompassed the centroid of an adjacent pastoral surface, a new zone Z was defined. The sheep availability of a zone Z was calculated as the sum of the sheep availability >(?) of all the pastoral surfaces u whose centroids were in Z. We determined the limit size of a zone Z as following: a zone Z could not include more than 5% of the sum of >(?) of all pastoral surfaces. Otherwise, the hotspots were too vast to be informative [START_REF] Kulldorff | A spatial scan statistic[END_REF]. Because information on the pastoral area was omitted in this analysis, we included this information in the definition of the sheep availability of a pastoral surface >(?),

as >(?) = 121 ! × 4567 ! × @=7@ ! , where @=7@ ! was the surface u in km². Therefore, the larger a pastoral surface, the larger its risk of depredation. We indeed considered that two flocks of the same size and grazing time should have a different depredation risk depending on the size of their pastures. Farmers generally used the entire pasture for their flock, either by dividing the flock into smaller flocks or by moving the flock within the pasture. Increasing the pasturing area should make the flock more exposed to depredation by increasing the risk of encounters with wolves. In addition, we restricted the analysis to the depredated pastoral surfaces. The inclusion of pastoral surfaces without depredation into the analysis would indeed have increased sheep availabilities of certain zones, especially zones where depredated pastoral surfaces were surrounded by pastoral surfaces without depredation. This would have reduced the capacity of the analysis to identify such zones as hotspots.

Second, the Kulldorff statistic modelled the total number of observed depredations n as:

A ~ Poisson(1 % × >(? ∈ H) + I % × >(? ∉ H))
with the probabilities 1 % and I % that a depredation event occurred within or outside Z, respectively, and >(? ∈ H) and >(? ∉ H) the sum of >(?) inside and outside Z, respectively.

For each Z, the Kulldorff statistic tested the null hypothesis 1 % = I % (M0) and the alternative hypothesis 1 % > I % (MA). It calculated the likelihood ratio LM(H) as

L(H) L & ⁄ . L(H)
corresponded to the result of the likelihood function L of the model under MA for the zone Z.

L & corresponded to the result of the likelihood function L of the model under M0, which was the same for all zones because under M0 the model can be reduced to A ~ Poisson(1 × ∑ >(?)). The zones Z for which MA was the most likely had the highest LM(H).

Third, we simulated 499 sets of randomly distributed n controls over the depredated pastoral surfaces, proportionally to their sheep availabilities. The second step was repeated for each control set, to eventually compute the control distribution of the highest likelihood ratio

LM(H).

If the observed highest LM(H) were among the top 5% of this control distribution, then the corresponding zones were significant hotspots.

We repeated this analysis with a simulated homogeneous sheep availability to show the discrepancies when ignoring prey availability: each depredated pastoral surface had then a sheep availability equal to the mean observed sheep availability of the depredated pastoral surfaces.

Results

Trends in depredation risk and observed depredations

We collected 25 220 wolf presence signs and 18 764 verified depredations on any type of

Analysis 1: Summarised spatial structure of wolf depredations

At the regional scale, the Kinhom function identified significant aggregation of the depredations for all years, with some Kinhom values higher than the control envelope (Fig. 1; Fig. B.4). The control envelope was thinner with time as the depredation pattern included more and more events. However, depending on years, the range of radiuses r for which aggregation was significant varied, as well as the aggregation intensity (i.e. difference between the observed and control envelope values at the same r). We identified four time periods.

The first period ranged from 1995 to 1999. In 1995, aggregation was significant up to 20 km.

The depredation pattern exhibited randomness beyond. From 1996, aggregation was significant for the whole range of computed radiuses r, up to 15-30 km around depredations.

Aggregation tended to intensify with r, especially in 1996 and 1998. This suggested that the

Analysis 2: Wolf depredation hotspot location

The number of hotspots identified by the Kulldorff statistic adjusting for sheep availability tended to increase over time in the study area, from 2 hotspots in 1995 to 32 in 2018 (Fig. 2A). The annual proportion of depredated pastoral surfaces into hotspots was higher for the period 2012-2018 (mean=29.2%, sd=3.8%) than for the period 1995-2011 (mean=23.5%, sd=5.9%) (Fig. 2B). Hotspots were of various sizes but generally did not exceed 50 km² of pastoral area (Fig. 2C). The average number of pastoral surfaces within hotspots tended to increase with time, from 3.0 in 1995 (sd=0) to 7.8 in 2018 (sd=11.5). Hotspots could be located at different places within the study area (Fig. 3A In contrast to the analysis adjusting for sheep availability (hereafter 'first Kulldorff analysis'), the results of hotspot identification changed in the analysis with simulated homogeneous sheep availability (hereafter 'second Kulldorff analysis'). For almost all years, the numbers of hotspots were lower in the second Kulldorff analysis than in the first one (Fig. 2A). The annual proportion of depredated pastoral surfaces into hotspots was always lower than in the second Kulldorff analysis and exceeded 20% only in 2014 and 2017 (Fig. 2B). The pastoral surfaces identified as hotspots in the second Kulldorff analysis were generally identified as hotspots in the first Kulldorff analysis (e.g.

Discussion

To reduce the number of depredations in an area, managers seek to act on the factors driving depredation. Depredations primarily rely on the rate of encounters between predators and livestock [START_REF] Mech | Wolves: Behavior, Ecology, and Conservation[END_REF], and therefore on livestock availability [START_REF] Balbuena-Serrano | Hotspots of livestock depredation by pumas and jaguars in Brazil: a biome-scale analysis[END_REF][START_REF] Bradley | Assessing factors related to wolf depredation of cattle in fenced pastures in Montana and Idaho[END_REF]. The room for maneuver on livestock availability is however generally limited for managers. Other factors can favour depredations, such as environmental, husbandry or predator-linked factors, or a combination of these factors [START_REF] Fowler | Non-linear relationships between human activities and wolf-livestock depredations[END_REF]. There is therefore a strong interest for managers to easily identify the areas where the high number of depredations is not only caused by the primary livestock availability factor, but also by other factors.

Depredation spatial risk modelling studies have generally integrated livestock availability.

However, the produced risk maps do not allow managers to clearly disentangle between areas with significant depredation clustering given livestock availability (e.g. [START_REF] Balbuena-Serrano | Hotspots of livestock depredation by pumas and jaguars in Brazil: a biome-scale analysis[END_REF][START_REF] Fowler | Non-linear relationships between human activities and wolf-livestock depredations[END_REF]. If adjusted for livestock availability, the depredation hotspot analysis can allow for such identification. In our case study, we showed that adjusting or not for livestock availability modified the results, whether it was for the summarised depredation pattern or the hotspot map. These two analyses are complementary to help managers to identify the causes favouring depredations.

The hotspot map provided the locations of areas where investigations in situ can be led to understand the unexpected high level of depredations. The most straightforward factor to investigate is husbandry. Generally, detailed husbandry practices or implementation of preventive measures are not available for spatial statistical analyses, as in our case study. Because a regional inventory of preventive measures is often difficult if not impossible to put in place, especially in the long term, it can be necessary to go in the field to check if preventive measures are implemented, but also suited to the local context. The involved farmers can receive technical support provided by the public authorities in the implementation of appropriate preventive tools [START_REF] Littlewood | Terrestrial Mammal Conservation: Global Evidence for the Effects of Interventions for Terrestrial Mammals Excluding Bats and Primates[END_REF]. Other factors can encompass larger spatial scales than the sole pasture (e.g. environmental factors as type of vegetation, or wolves). Investigation of such factors is difficult if the hotspot definition is based on each spatial unit without considering the other units, as it is currently done by the French administration. Our approach has the benefit of identifying groups of spatial units belonging to the same hotspot, and therefore to facilitate the investigation of large-scale factors.

The analysis of the summarised depredation pattern can also help to understand the spatial scale of the depredation pattern, and therefore the scale relevant to management. The hotspot radiuses were estimated between 12 and 50 km in our Kinhom analysis at the local and regional scales. This did not match with the 7-8 km wolf territory radius estimated in France through telemetry and genetic tracking [START_REF] Duchamp | A dual frame survey to assess time-and space-related changes of the colonizing wolf population in France[END_REF]. Therefore, these results suggest that most of the hotspots encompassed more than one pasture, but also more than one wolf territory. In other words, the depredation clustering was not linked to one 'problematic' wolf or pack, or one isolated unprotected flock. Instead, these results suggest that a local context could favour depredations, e.g. a lack of wild preys or environmental conditions making an efficient guarding difficult (as bushy environments). The Kinhom analysis also showed that depredation pattern tended to randomise through time. This randomisation pattern did not seem to be linked to the increase of the area recolonised by wolves, as the randomisation occurred at the local scale long after wolves recolonised the MNP. The depredation risk homogenisation among pastoral surfaces could result from the increasing number of farmers who protected their flocks against wolf depredations during wolf recolonisation [START_REF] Mtes | État des lieux et diagnostic de la situation dans la perspective de la rédaction du futur plan national d'actions 2018-2023 sur le loup et les activités d'élevage[END_REF], as suggested by [START_REF] Meuret | Missing shots: has the possibility of shooting wolves been lacking for 20 years in France's livestock protection measures? Rangel[END_REF].

Besides the identification of the causes favouring depredation, the results provided by the hotspots maps can help in resource allocation. For example, France currently allocates its highest level of subsidies for preventive measures in depredation hotspots. Each year hotspots are defined as the municipality with a mean of at least 15 depredations over the last three years, and all its adjacent municipalities. This definition accounts neither for the surface of the municipality nor for the surface of pastures. This definition creates biases, large municipalities having higher chances to be identified as hotspots than others (DREAL Auvergne-Rhône-Alpes et al., 2020). Moreover, relying on three years to define hotspots prevents the detection of emerging hotspots in time. In our case study, the methodology we used corrected these biases. Correcting for livestock availability also allowed us to detect hotspots of pastoral surfaces with low sheep counts and grazing time. The current administrative definition of hotspots does not account for the size of flocks or the time spent on pastures. Therefore, hotspots with low livestock availability can remain undetected by the administration. Yet, the consequences of depredations for the involved livestock owners are high in these situations because a high proportion of the flock is lost to predation [START_REF] Stahl | Predation on livestock by an expanding reintroduced lynx population: long-term trend and spatial variability[END_REF].

Thus, the spatial correlation analyses are powerful tools for large carnivore conservation.

They help to better define priority conservation areas. These areas can receive a specific attention from managers, whether it is for preventive measures subsidies or for investigating the effective implementation of these measures, as recommended by [START_REF] Bautista | Large carnivore damage in Europe: Analysis of compensation and prevention programs[END_REF]. The methodology we presented could then optimize the cost-effectiveness of public funds. The For both the Kinhom and the Kulldorff analyses, information about livestock availability was of paramount importance to detect the depredation pattern and hotspot locations. This information is usually hard to obtain, even more at a regional scale. This explains why the studies on depredation hotspots generally do not adjust for livestock information (e.g. [START_REF] Hoffmann | Spatial Pattern Analysis Reveals Randomness Among Carnivore Depredation of Livestock[END_REF] or only for livestock distribution (e.g. [START_REF] Gastineau | Spatiotemporal depredation hotspots of brown bears, Ursus arctos, on livestock in the Pyrenees, France[END_REF]. Here, our information on livestock availability was limited in time to two censuses which in total covered a period of 5 years. Limiting our study period to this 5-year period would have been possible and still informative to show the discrepancies in the results of Kinhom and of the Kulldorff statistic when adjusting or not for livestock availability. However, the extension of the use of livestock information to years surrounding the census' periods brought valuable information about the dynamics of the depredation pattern. Indeed, grazing time and flock size are highly constrained by the climate that differs greatly between the north and the south of the studied area. If it is possible that flocks may have grazed more or less before or after each census, we do not expect these differences to overcome the benefits of using the pastoral censuses that allowed the integration of seasonal differences between pastures in the analysis.

Conclusions

We encourage the use of clustering analyses adjusted for livestock availability in the context of large predator conservation and of the management of their predations on livestock.

Contrary to the conventional wisdom which relegates clustering analyses to a preliminary step in risk modelling, the estimation of hotspot locations is a full-fledged decision-making tool.

Hotspots are the first targets for allocation of conservation resources, or for investigation in situ of the causes favouring depredation. The livestock owners in hotspots can be supported, financially or technically, to quickly adapt their husbandry practices according to the pastures they use and the predator they face. Using spatial statistical analyses with Kinhom or with the Kulldorff statistic combined with livestock availability data would substantially improve the reliability of the estimates of depredation pattern and hotspot locations, and therefore the use of these quantities for conservation.

Figures

Color should be used for all figures in print.

  and Auvergne-Rhône-Alpes, and included all the French Alps and the east part of Massif Central mountains (102 483 km², Fig. A.1). Apart from the highest altitudes which are without forests, habitats range from bush to coniferous forest under Mediterranean climate in the south, to mixed forest in the north and the east.

  km², Fig. A.1), i.e. the local scale. The local scale overlapped the Mercantour National Park and its surroundings (hereafter 'MNP'), which is a mountainous area characterised by a succession of alpine vegetation levels along a wide altitudinal range, from 600 to 3 200 m, and by an important pastoral activity for meat production. The MNP is the oldest place where wolves and pastoralism are co-occurring in the country since wolf recolonisation in the 1990's from Italy.

  B.1). For both censuses, the spatio-temporal sheep availability was inhomogeneous across the study area (Fig. B.1), particularly for grazing time

Fig B. 2 ).

 2 Fig B.2).The annual proportion of pastoral surfaces at depredation risk within the total set of

  depredation pattern showed no particular spatial structure, except in 1995. During the second period from 2000 to 2003, aggregation intensity decreased at the largest computed r values.From 2001, the Kinhom values merged with the control envelope between 40 and 50 km, suggesting a locally emerging clustered pattern. During the third period from 2004 to 2012, a weak significant aggregation was only observed over the lowest r values, and turned into significant repulsion between 15 and 35 km. This was typical of a clustered pattern, with hotspots producing coldspot interzones[START_REF] Hoffmann | Spatial Pattern Analysis Reveals Randomness Among Carnivore Depredation of Livestock[END_REF]. Finally, the last period from 2013 to 2018 showed significant aggregation over 25 to 40 km followed by randomness (except 2014 with continuous aggregation). In certain years, a weak repulsion was observed at the largest r values. Aggregation intensity was greater in 2013 compared to 2012, but it tended to decrease the following years, which brought the depredation pattern closer to randomness with time.Using Kinhom or the K-function provided comparable results only during the first period. But, from 1999-2000, the differences between the analysis of Kinhom and the two analyses of the Kfunction increased (Fig.1; Fig. B.4). The first analysis with the K-function which considered sheep availability as homogeneous across pastoral surfaces identified the typical hotspot structure (i.e. aggregation followed by repulsion) only in 2009 and 2010 with a radius of 30 and 40 km respectively. The other years, aggregation was significant for all values of r and always more intense than the one obtained through Kinhom. When the analysis grid was only based on the wolf distribution, the K-function analysed the pattern over larger values of r than the previous analyses because isolated distribution cells in the west of the study area could sporadically appear without overlapping with pastoral surfaces. This analysis only identified significant aggregation which, in most years, intensified with r. At the local scale, the Kinhom function provided two main findings. First, from 1995 to 2011, the analysis identified significant aggregation over 12 to 17 km, occasionally followed by randomness at large radiuses (Fig. B.4). Repulsion was never observed, except in 2004 where depredations seemed clustered into three single zones (Fig. B.2). Second, the aggregation intensity was generally decreasing with time from 2013, even though the control envelopes were of equivalent thickness. The pattern tended to randomness for all values of r from 2015, with almost complete randomness observed in 2016. Therefore, the depredations were increasingly widespread at the local scale and were no longer significantly clustering in the last years of the study period, long after the full recolonisation of the MNP by wolves.

  ; Fig. B.5). They could be sporadic, appearing for one year. But others could persist over time, such as in the MNP where hotspots were present during the whole study period, and in the south for the period 2013-2018.

  Fig. 3B; Fig. B.5). The unidentified hotspots in the second Kulldorff analysis generally encompassed small pastoral surfaces where sheep availability was low as well as the number of depredations. Thus, annual means of sheep counts, of grazing time and of depredations of the pastoral surfaces within hotspots were significantly larger in the second Kulldorff analysis than in the first one (Wilcoxon tests: P<0.001, α=0.05).

  user interface we developed to present the hotspot maps (Fig B.5) should also facilitate the adoption of the tool by managers.

Figure 1 :

 1 Figure 1: Results at the regional scale of the Kinhom function (solid yellow line) and of the K-

Figure 2 :

 2 Figure 2: Descriptive results of the Kulldorff statistic used with the observed heterogeneous

Figure 3 :

 3 Figure 3: Locations of depredated pastoral surfaces identified as hotspots (dark colour) or not

Fig. B. 3 ̶

 3 Fig. B.3 ̶ Annual distributions of the number of wolf depredations on sheep per pastoral surface, at the regional and local scales between biological years 1995 and 2018. https://oksanagrente.shinyapps.io/Hotspots-SuppFig3/ Fig. B.4 ̶ Annual results between biological years 1995 and 2018 of the Kinhom analysis and of the two K-function analyses (one with the analysis grid simulated over the pastoral surfaces at depredation risk, the other only over the wolf distribution) at the regional and local scales. https://oksanagrente.shinyapps.io/Hotspots-SuppFig4/ Fig. B.5 ̶ Annual locations between biological years 1995 and 2018 of depredated pastoral surfaces identified or not as hotspots, according to the Kulldorff statistic used with the observed heterogeneous sheep availability or the simulated homogeneous sheep availability at the regional scale. https://oksanagrente.shinyapps.io/Hotspots-SuppFig5/

Acknowledgements

This work was partly funded by a grant from 'Mission pour l'interdisciplinarité' of CNRS, through its 'Osez l'interdisciplinarité' call. We are indebted to the governmental agents in charge of recording depredations and pastoral information and to the field experts from the Wolf-Lynx French Network who collected wolf presence signs. We are grateful to IRSTEA and the regional administration of Auvergne-Rhône-Alpes for the access to the pastoral and depredation databases. We thank A. Baddeley and E. Gabriel for valuable methodological advice, G. Loucougaray, F. Bray and M. Meuret from INRAE and L. Garde from CERPAM for advice on the pastoralism overview and diagnostic. We thank N. Drouet-Hoguet and S.

Bauduin for comments on the first draft of the manuscript.

Supplementary material Appendix A -Figures


Olivier Gimenez

Maëlis Kervellec

Jean-Baptiste Fanjul

Anna Chaine

Lucile Marescot

Yoann Bollet

Christophe Duchamp

J.-B Fanjul

Trade-off between deep learning for species identification and inference about predator-prey co-occurrence: Reproducible R workflow integrating models in computer vision and ecological statistics

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Computer vision is a field of artificial intelligence in which a machine is taught how to extract and interpret the content of an image [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF]. Computer vision relies on deep learning that allows computational models to learn from training data -a set of manually labelled images -and make predictions on new data -a set of unlabelled images [START_REF] Baraniuk | The Science of Deep Learning[END_REF][START_REF] Lecun | Deep Learning[END_REF]. With the growing availability of massive data, computer vision with deep learning is being increasingly used to perform tasks such as object detection, face recognition, action and activity recognition or human pose estimation in fields as diverse as medicine, robotics, transportation, genomics, sports and agriculture [START_REF] Voulodimos | Deep Learning for Computer Vision: A Brief Review[END_REF].

In ecology in particular, there is a growing interest in deep learning for automatizing repetitive analyses on large amounts of images, such as identifying plant and animal species, distinguishing individuals of the same 1 arXiv:2108.11509v1 [stat.AP] 25 Aug 2021 or different species, counting individuals or detecting relevant features [START_REF] Christin | Applications for Deep Learning in Ecology[END_REF][START_REF] Lamba | Deep Learning for Environmental Conservation[END_REF][START_REF] Weinstein | A Computer Vision for Animal Ecology[END_REF]. By saving hours of manual data analyses and tapping into massive amounts of data that keep accumulating with technological advances, deep learning has the potential to become an essential tool for ecologists and applied statisticians.

Despite the promising future of computer vision and deep learning, there are challenging issues toward their wide adoption by the community of ecologists (e.g. [START_REF] Wearn | Responsible AI for Conservation[END_REF]. First, there is a programming barrier as most, if not all, algorithms are written in the Python language while most ecologists are versed in R [START_REF] Lai | Evaluating the Popularity of R in Ecology[END_REF]. If ecologists are to use computer vision in routine, there is a need for bridges between these two languages (through, e.g., the reticulate package [START_REF] Allaire | Reticulate: R Interface to Python[END_REF] or the shiny package [START_REF] Tabak | Improving the Accessibility and Transferability of Machine Learning Algorithms for Identification of Animals in Camera Trap Images: Mlwic2[END_REF]). Second, recent applications of computer vision via deep learning in ecology have focused on computational aspects and simple tasks without addressing the underlying ecological questions [START_REF] Sutherland | Identification of 100 Fundamental Ecological Questions[END_REF], or carrying out statistical data analysis to answer these questions [START_REF] Gimenez | Statistical Ecology Comes of Age[END_REF]. Although perfectly understandable given the challenges at hand, we argue that a better integration of the why (ecological questions), the what (automatically labelled images) and the how (statistics) would be beneficial to computer vision for ecology (see also [START_REF] Weinstein | A Computer Vision for Animal Ecology[END_REF].

Here, we showcase a full why-what-how workflow in R using a case study on the structure of an ecological community (a set of co-occurring species) composed of the Eurasian lynx (Lynx lynx) and its two main preys. First, we introduce the case study and motivate the need for deep learning. Second we illustrate deep learning for the identification of animal species in large amounts of images, including model training and validation with a dataset of labelled images, and prediction with a new dataset of unlabelled images. Last, we proceed with the quantification of spatial co-occurrence using statistical models.

Collecting images with camera traps

Lynx (Lynx lynx) went extinct in France at the end of the 19th century due to habitat degradation, human persecution and decrease in prey availability [START_REF] Vandel | Distribution Trend of the Eurasian Lynx Lynx Lynx Populations in France[END_REF]. The species was reintroduced in Switzerland in the 1970s [START_REF] Breitenmoser | Large Predators in the Alps: The Fall and Rise of Man's Competitors[END_REF], then re-colonised France through the Jura mountains in the 1980s [START_REF] Vandel | Distribution Trend of the Eurasian Lynx Lynx Lynx Populations in France[END_REF]. The species is listed as endangered under the 2017 IUCN Red list and is of conservation concern in France due to habitat fragmentation, poaching and collisions with vehicles. The Jura holds the bulk of the French lynx population.

To better understand its distribution, we need to quantify its interactions with its main preys, roe deer (Capreolus capreolus) and chamois (Rupicapra rupicapra) [START_REF] Molinari-Jobin | Variation in Diet, Prey Selectivity and Home-Range Size of Eurasian Lynx Lynx Lynx in Switzerland[END_REF], two ungulate species that are also hunted. To assess the relative contribution of predation and hunting, a predator-prey program was set up jointly by the French Office for Biodiversity, the Federations of Hunters from the Jura, Ain and Haute-Savoie counties and the French National Centre for Scientific Research.

Animal detections were made using a set of camera traps in the Jura mountains that were deployed in the Jura and Ain counties (see Figure 1). We divided the two study areas into grids of 2.7 × 2.7 km cells or sites hereafter [START_REF] Zimmermann | Optimizing the Size of the Area Surveyed for Monitoring a Eurasian Lynx (Lynx Lynx) Population in the Swiss Alps by Means of Photographic Capture-Recapture[END_REF]) in which we set two camera traps per site (Xenon white flash with passive infrared trigger mechanisms, model Capture, Ambush and Attack; Cuddeback), with 18 sites in the Jura study area, and 11 in the Ain study area that were active over the study period (from February 2016 to October 2017 for the Jura county, and from February 2017 to May 2019 for the Ain county). The location of camera traps was chosen to maximise lynx detection. Camera traps were checked weekly to change memory cards, batteries and to remove fresh snow after heavy snowfall.

In total, 45563 and 18044 pictures were considered in the Jura and Ain sites respectively after manually droping empty pictures and pictures with unidentified species. Note that classifying empty images could be automatised with deep learning (Norouzzadeh et al. 2021;[START_REF] Tabak | Improving the Accessibility and Transferability of Machine Learning Algorithms for Identification of Animals in Camera Trap Images: Mlwic2[END_REF]. We identified the species present on all images by hand (see Table 1) using digiKam a free open-source digital photo management application (https://www.digikam.org/). This operation took several weeks of labor full time, which is often identified as a limitation of camera trap studies. To expedite this tedious task, computer vision with deep learning has been identified as a promising approach (Norouzzadeh et al. 2021;[START_REF] Tabak | Machine Learning to Classify Animal Species in Camera Trap Images: Applications in Ecology[END_REF][START_REF] Willi | Identifying Animal Species in Camera Trap Images Using Deep Learning and Citizen Science[END_REF].

Deep learning for species identification

Using the images we obtained with camera traps (Table 1), we trained a model for identifying species using the Jura study site as a calibration dataset. We then assessed this model's ability to automatically identify species on a new dataset, also known as transferability, using the Ain study site as an evaluation dataset.

Training -Jura study site

We selected at random 80% of the annotated images for each species in the Jura study site for training, and 20% for testing. We applied various transformations (flipping, brightness and contrast modifications; [START_REF] Shorten | A Survey on Image Data Augmentation for Deep Learning[END_REF]) to improve training (see Appendix). To reduce model training time and overcome the small number of images, we used transfer learning [START_REF] Yosinski | How Transferable Are Features in Deep Neural Networks?[END_REF][START_REF] Shao | Transfer Learning for Visual Categorization: A Survey[END_REF] and considered a pre-trained model as a starting point. Specifically, we trained a deep convolutional neural network (ResNet-50) architecture [START_REF] He | Deep Residual Learning for Image Recognition[END_REF]) using the fastai library (https://docs.fast.ai/) that implements the PyTorch library [START_REF] Paszke | PyTorch: An Imperative Style, High-Performance Deep Learning Library[END_REF]. Interestingly, the fastai library comes with an R interface (https://eagerai.github.io/fastai/) that uses the reticulate package to communicate with Python, therefore allowing R users to access up-to-date deep learning tools. We trained models on the Montpellier Bioinformatics Biodiversity platform using a GPU machine (Titan Xp nvidia) with 16Go of RAM. We used 20 epochs which took approximately 10 hours. The computational burden prevented us from providing a full reproducible analysis, but we do so with a subsample of the dataset in the Appendix. All trained models are available from https://doi.org/10.5281/zenodo.5164796.

We calculated three metrics to evaluate our model performance at correctly identifying species (e.g. [START_REF] Duggan | An Approach to Rapid Processing of Camera Trap Images with Minimal Human Input[END_REF]. Specifically, we relied on accuracy the ratio of correct predictions to the total number of predictions, recall a measure of false negatives (FN; e.g. an image with a lynx for which our model predicts another species) with recall = TP / (TP + FN) where TP is for true positives, and precision a measure of false positives (FP; e.g. an image with any species but a lynx for which our model predicts a lynx) with precision = TP / (TP + FP). In camera trap studies, a strategy [START_REF] Duggan | An Approach to Rapid Processing of Camera Trap Images with Minimal Human Input[END_REF]) consists in optimizing precision if the focus is on rare species (lynx), while recall should be optimized if the focus is on commom species (chamois and roe deer).

We achieved 85% accuracy during training. Our model had good performances for the three classes we were interested in, with 87% precision for lynx and 81% recall for both roe deer and chamois (Table 2).

Transferability -Ain study site

We evaluated transferability for our trained model by predicting species on images from the Ain study site which were not used for training. Precision was 77% for lynx, and while we achieved 86% recall for roe deer, 3).

To better understand this pattern, we display the results under the form of a confusion matrix that compares model classifications to manual classifications (Figure 2). There were a lot of false negatives for chamois, meaning that when a chamois was present in an image, it was often classified as another species by our model.

Overall, our model trained on images from the Jura study site did poorly at correctly predicting species on images from the Ain study site. This result does not come as a surprise, as generalizing classification algorithms to new environments is known to be difficult [START_REF] Beery | Recognition in Terra Incognita[END_REF]. While a computer scientist might be disappointed in these results, an ecologist would probably wonder whether ecological inference about the interactions between lynx and its prey is biased by these average performances, a question we address in the next section.

Spatial co-occurrence

Here, we analysed the data we acquired from the previous section. For the sake of comparison, we considered two datasets, one made of the images manually labelled for both the Jura and Ain study sites pooled together (ground truth dataset), and the other in which we pooled the images that were manually labelled for the Jura study site and the images that were automatically labelled for the Ain study site using our trained model (classified dataset).

We formatted the data by generating monthly detection histories, that is a sequence of detections (Y sit = 1) and non-detections (Y sit = 0), for species s at site i and sampling occasion t (see Figure 3).

To quantify spatial co-occurrence betwen lynx and its preys, we used a multispecies occupancy modeling approach [START_REF] Rota | A Multispecies Occupancy Model for Two or More Interacting Species[END_REF][START_REF] Clipp | A Penalized Likelihood for Multi-Species Occupancy Models Improves Predictions of Species Interactions[END_REF]) using the R package unmarked [START_REF] Fiske | unmarked: An R Package for Fitting Hierarchical Models of Wildlife Occurrence and Abundance[END_REF].

The multispecies occupancy model assumes that observations y sit , conditional on Z si the latent occupancy state of species s at site i are drawn from Bernoulli random variables Y sit |Z si ∼ Bernoulli(Z si p sit ) where p sit is the detection probability of species s at site i and sampling occasion t. Detection probabilities can be modeled as a function of site and/or sampling covariates, or the presence/absence of other species, but for the sake of illustration, we will make them only species-specific here.

The latent occupancy states are assumed to be distributed as multivariate Bernoulli random variables (Dai, Ding, and Wahba 2013). Let us consider 2 species, species 1 and 2, then Z i = (Z i1 , Z i2 ) ∼ multivariate Bernoulli(ψ 11 , ψ 10 , ψ 01 , ψ 00 ) where ψ 11 is the probability that a site is occupied by both species 1 and 2, ψ 10 the probability that a site is occupied by species 1 but not 2, ψ 01 the probability that a site is occupied by species 2 but not 1, and ψ 00 the probability a site is occupied by none of them. Note that we considered species-specific only occupancy probabilities but these could be modeled as site-specific covariates. Marginal occupancy probabilities are obtained as Pr(Z i1 = 1) = ψ 11 + ψ 10 and Pr(Z i2 = 1) = ψ 11 + ψ 01 . With this model, we may also infer potential interactions by calculating conditional probabilities such as for example the probability of a site being occupied by species 2 conditional of species 1 with Pr(Z i2 = 1|Z i1 = 1) = ψ 11 ψ 11 + ψ 10 .

Detection probabilities were indistinguishable whether we used the ground truth or the classified dataset, with p lynx = 0.51(0.45, 0.58), p roe deer = 0.63(0.57, 0.68) and p chamois = 0.61(0.55, 0.67).

We also found that occupancy probability estimates were similar whether we used the ground truth or the classified dataset (Figure 4). Roe deer was the prevalent species, but lynx and chamois were also occurring with high probability (Figure 4). Because marginal occupancy probabilities were high, probabilities of co-occurrence were also estimated high (Figure 5). Our results should be interpreted bearing in mind that co-occurrence is a necessary but not sufficient condition for actual interaction. When both preys were present, lynx was more present than when they were both absent (Figure 5). Lynx was more sensitive to the presence of roe deer than that of chamois (Figure 5).

Discussion

In this paper, we aimed at illustrating a reproducible workflow for studying the structure of an animal community and species spatial co-occurrence (why) using images acquired from camera traps and automatically labelled with deep learning (what) which we analysed with statistical occupancy models accounting for imperfect species detection (why). Overall, we found that, even though model transferability could be improved, inference about the potential interactions between lynx and its preys was similar whether we analysed the ground truth data or classified data.

This result calls for further work on the trade-offs between time and resources allocated to train models with deep learning and our ability to correctly answer key ecological questions with camera-trap surveys.

In other words, while a computer scientist might be keen on spending time training models to achieve top performances, an ecologist would rather rely on a model showing average performances and use this time to proceed with statistical analyses if, of course, errors in computer-annotated images do not make ecological inference flawed. The right balance may be found with collaborative projects in which scientists from artificial intelligence, statistics and ecology agree on a common objective, and identify research questions that can pick the interest of all parties.

Our demonstration remains however empirical, and we encourage others to try and replicate our results. We also see two avenues of research that could benefit the integration of deep learning and ecological statistics.

First, a simulation study could be conducted to evaluate bias and precision in ecological parameter estimators with regard to errors in image annotation by computers. The outcome of this exercise could be, for example, guidelines informing on the confidence an investigator may place in ecological inference as a function of the amount of false negatives and false positives. Second, annotation errors could be accomodated directly in statistical models. For example, single-species occupancy models account for false negatives when a species is not detected by the camera at a site where it is present, as well as false positives when a species is detected at a site where it is not present due to species misidentification by the observer [START_REF] Miller | Improving Occupancy Estimation When Two Types of Observational Error Occur: Non-Detection and Species Misidentification[END_REF]. Pending a careful distinction between ecological vs. computer-generated false negatives and false positives, error rates could be added to multispecies occupancy models [START_REF] Chambert | Two-Species Occupancy Modelling Accounting for Species Misidentification and Non-Detection[END_REF]) and informed by recall and precision metrics obtained during model training [START_REF] Tabak | Improving the Accessibility and Transferability of Machine Learning Algorithms for Identification of Animals in Camera Trap Images: Mlwic2[END_REF].

With regard to the case study, our results should be seen only as preliminary. First, we aim at quantifying the relative contribution of biotic (lynx predation on chamois and roe deer) and abiotic (habitat quality) processes to the composition and dynamic of this ecological community. Second, to benefit future camera trap studies of lynx in the Jura mountains, we plan to train a model again using manually annotated images from both the Jura and the Ain study sites. These perspectives are the object of ongoing work.

With the rapid advances in technologies for biodiversity monitoring (Lahoz-Monfort and Magrath 2021), the possibility of analysing large amounts of images makes deep learning appealing to ecologists. We hope that our proposal of a reproducible R workflow for deep learning and statistical ecology will encourage further studies in the integration of these disciplines, and contribute to the adoption of computer vision by ecologists.

Appendix: Reproducible example of species identification on camera trap images with CPU

In this section, we go through a reproducible example of the entire deep learning workflow, including data preparation, model training, and automatic labeling of new images. We used a subsample of 467 images from the original dataset in the Jura county to allow the training of our model with CPU on a personal computer. We also used 14 images from the original dataset in the Ain county to illustrate prediction.

Training and validation datasets

We first split the dataset of Jura images in two datasets, a dataset for training, and the other one for validation. We use the exifr package to extract metadata from images, get a list of images names and extract the species from these. )) %>% # when pix has no tag select(SourceFile, FileName, Keywords) %>% mutate(Keywords = fct_recode(Keywords, "chat" = "chat forestier", "lievre" = "lièvre", "vehicule" = "véhicule", "ni" = "Non identifié")) %>% filter(!(Keywords %in% c("ni", "wo_tag")))

Then we pick 80% of the images for training in each category, the rest being used for validation.

# training dataset pix_train <-labels %>% select(SourceFile, FileName, Keywords) %>% group_by(Keywords) %>% filter(between(row_number(), 1, floor(n()*80/100))) # 80% per category # validation dataset pix_valid <-labels %>% group_by(Keywords) %>% filter(between(row_number(), floor(n()*80/100) + 1, n()))

Eventually, we store these images in two distinct directories named train and valid.

# create dir train/ and copy pix there, organised by categories dir.create( pix/train ) # create training directory for (i in levels(fct_drop(pix_train$Keywords))) dir.create(paste0( pix/train/ ,i)) # create dir for labe for (i in 1:nrow(pix_train)){ file.copy(as.character(pix_train$SourceFile[i]), paste0( pix/train/ , as.character(pix_train$Keywords[i]))) # copy pix in corresp dir } # create dir valid/ and copy pix there, organised by categories. dir.create( pix/valid ) # create validation dir for (i in levels(fct_drop(pix_train$Keywords))) dir.create(paste0( pix/valid/ ,i)) # create dir for labe for (i in 1:nrow(pix_valid)){ file.copy(as.character(pix_valid$SourceFile[i]), paste0( pix/valid/ , as.character(pix_valid$Keywords[i]))) # copy pix in corresp dir } # delete pictures in valid/ directory for which we did not train the model to_be_deleted <-setdiff(levels(fct_drop(pix_valid$Keywords)), levels(fct_drop

(pix_train$Keywords))) if (!is_empty(to_be_deleted)) { for (i in 1:length(to_be_deleted)){ unlink(paste0( pix/valid/ , to_be_deleted[i])) } }
What is the sample size of these two datasets? bind_rows("training" = pix_train, "validation" = pix_valid, .id = "dataset") %>% group_by(dataset) %>% count(Keywords) %>% rename(category = Keywords) %>% kable(caption = "Sample size (n) for the training and validation datasets.") %>% kable_styling()

Transfer learning

We proceed with transfer learning using images from the Jura county (or a subsample more exactly). We first load images and apply standard transformations to improve training (flip, rotate, zoom, rotate, ligth transform).

dls <-ImageDataLoaders_from_folder( path = "pix/", train = "train", valid = "valid", item_tfms = Resize(size = 460), bs = 10, batch_tfms = list(aug_transforms(size = 224, min_scale = 0.75), # transformation Normalize_from_stats( imagenet_stats() )), num_workers = 0, ImageFile.LOAD_TRUNCATED_IMAGES = TRUE)

Then we get the model architecture. For the sake of illustration, we use a resnet18 here, but we used a resnet50 to get the full results presented in the main text. ----------------------------------------------------#

Transferability

In this section, we show how to use our freshly trained model to label images that were taken in another study site in the Ain county, and not used to train our model. First, we get the path to the images.

fls <-list.files(path = "pix/pixAin", full.names = TRUE, recursive = TRUE)

Then we carry out prediction, and compare to the truth.

predicted <-character(3) categories <-interp$vocab %>% str_replace_all ("[[:punct:]]", " ") %>% str_trim() %>% str_split(" ") %>%

Figure 2 :Figure 3 :

 23 Figure2: Confusion matrix comparing automatic to manual species classifications. Species that were predicted by our model are in columns, and species that are actually in the images are in rows.

Figure 4 :

 4 Figure 4: Marginal occupancy probabilities for all three species, lynx, roe deer and chamois). Parameter estimates are from a multispecies occupancy model using either the ground truth dataset (in red) or the classified dataset (in blue-grey).

Figure 5 :

 5 Figure 5: Lynx occupancy probability conditional on the presence or absence of its preys (roe deer and chamois). Parameter estimates are from a multispecies occupancy model using either the ground truth dataset (in red) or the classified dataset (in blue-grey).

  learn <-cnn_learner(dls = dls, arch = resnet18(),

Table 1 :

 1 Species identified in the Jura and Ain study sites with samples size (n). Only first 10 species with most images are shown.

					B. Ain county
					46.10°N
	A. Jura county
	46.80°N				46.05°N
	46.75°N				46.00°N
	46.70°N
					45.95°N
	46.65°N
					45.90°N
	46.60°N
					45.85°N
	5.8°E	5.9°E	6.0°E	6.1°E
					45.80°N
					5.35°E 5.40°E 5.45°E 5.50°E 5.55°E 5.60°E 5.65°E 5.70°E
		Figure 1: Study area, grid and camera trap locations.
		Species in Jura study site	n Species in Ain study site	n
		human			31644 human	4946
		vehicule			5637 vehicule	4454
		dog			2779 dog	2310
		fox			2088 fox	1587
		chamois			919 rider	1025
		wild board			522 roe deer	860
		badger			401 chamois	780
		roe deer			368 hunter	593
		cat			343 wild board	514
		lynx			302 badger	461

Table 2 :

 2 Model training performance. Images from the Jura study site were used for training.

	species	precision recall
	badger	0.78	0.88
	red deer	0.67	0.21
	chamois	0.86	0.81
	cat	0.89	0.78
	roe deer	0.67	0.81
	dog	0.78	0.84
	human	0.99	0.79
	hare	0.32	0.52
	lynx	0.87	0.95
	fox	0.85	0.90
	wild boar 0.93	0.88
	vehicule	0.95	0.98

Table 3 :

 3 Model transferability performance. Images from the Ain study site were used for assessing transferability.

		precision recall
	badger	0.71	0.89
	rider	0.79	0.92
	red deer	0.00	0.00
	chamois	0.82	0.08
	hunter	0.17	0.11
	cat	0.46	0.59
	roe deer	0.67	0.86
	dog	0.77	0.35
	human	0.51	0.93
	hare	0.37	0.35
	lynx	0.77	0.89
	marten	0.05	0.04
	fox	0.90	0.53
	wild board	0.75	0.94
	cow	0.01	0.25
	vehicule	0.94	0.51

our model performed poorly for chamois with 8% recall (Table

Table 4 :

 4 Species considered, and number of images with these species in them.

	Keywords	n
	humain	143
	vehicule	135
	renard	58
	sangliers	33
	chasseur	17
	chien	14
	lynx	13
	chevreuil	13
	chamois	12
	blaireaux	10
	chat	8
	lievre	4
	fouine	1
	cavalier	1

Table 5 :

 5 Sample size (n) for the training and validation datasets. Now we are ready to train our model. Again, for the sake of illustration, we use only 2 epochs here, but used 20 epochs to get the full results presented in the main text. With all pictures and a resnet50, it took 75 minutes per epoch approximatively on a Mac with a 2.4Ghz processor and 64Go memory, and less than half an hour on a machine with GPU. On this reduced dataset, it took a bit more than a minute per epoch on the same Mac. Note that we save the model after each epoch for later use.

	dataset	category	n
	training	humain	114
	training	vehicule	108
	training	chamois	9
	training	blaireaux	8
	training	sangliers	26
	training	renard	46
	training	chasseur	13
	training	lynx	10
	training	chien	11
	training	chat	6
	training	chevreuil	10
	training	lievre	3
	validation humain	29
	validation vehicule	27
	validation chamois	3
	validation blaireaux	2
	validation sangliers	7
	validation renard	12
	validation chasseur	4
	validation lynx	3
	validation chien	3
	validation fouine	1
	validation chat	2
	validation chevreuil	3
	validation lievre	1
	validation cavalier	1

  We may dig a bit deeper in training performances by loading the best model, here model_1.pth, and display some metrics for each species.We may extract the categories that get the most confused.

	## ##	(downsample): Sequential( (1): Sequential(
	## ##	(0): Conv2d(64, 128, kernel_size=(1, 1), stride=(2, 2), bias=False) (0): AdaptiveConcatPool2d(
	## ##		(1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (ap): AdaptiveAvgPool2d(output_size=1)
	## ##		) (mp): AdaptiveMaxPool2d(output_size=1)
	## ##	)	)
	## ##	(1): BasicBlock( (1): Flatten(full=False)
	## ##	(conv1): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (2): BatchNorm1d(1024, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
	## ##	(bn1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (3): Dropout(p=0.25, inplace=False)
	## ##	(relu): ReLU(inplace=True) (4): Linear(in_features=1024, out_features=512, bias=False)
	## ##	(conv2): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (5): ReLU(inplace=True)
	## ##	(bn2): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (6): BatchNorm1d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
	## ##	) (7): Dropout(p=0.5, inplace=False)
	## ##	) (8): Linear(in_features=512, out_features=12, bias=False)
	# 0 ## ## ## 1 ## ## ) one_cycle (6): Sequential( 2.447199 ) 1.736201 (0): BasicBlock( 0.850911 0.666055 ## ## (conv1): Conv2d(128, 256, kernel_size=(3, 3), stride=(2, 2), padding=(1, 1), bias=False) 0.760417 0.239583 01:24 0.781250 0.218750 01:25 ## interp <-ClassificationInterpretation_from_learner(learn) (bn1): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ## interp$print_classification_report() (relu): ReLU(inplace=True) epoch train_loss valid_loss accuracy error_rate ## 1 0 ## (conv2): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) 2.447199 0.8509111 0.7604167 0.2395833 ## 2 1 1.736201 0.6660554 0.7812500 0.2187500 ## (bn2): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ## (downsample): Sequential( interp %>% most_confused()
	learn$load("model_1") ## Sequential( ## (0): Sequential( ## (0): Conv2d(3, 64, kernel_size=(7, 7), stride=(2, 2), padding=(3, 3), bias=False) ## (1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ## (2): ReLU(inplace=True) ## (3): MaxPool2d(kernel_size=3, stride=2, padding=1, dilation=1, ceil_mode=False) ## (4): Sequential( ## (0): BasicBlock( ## (conv1): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) ## (bn1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ## (relu): ReLU(inplace=True) ## (conv2): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) ## (bn2): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## (bn2): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ## ) (conv2): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) ## ) (relu): ReLU(inplace=True) ## ) (bn1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ## (bn2): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (conv1): Conv2d(64, 128, kernel_size=(3, 3), stride=(2, 2), padding=(1, 1), bias=False) ## (conv2): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) (0): BasicBlock( ## (relu): ReLU(inplace=True) (5): Sequential( ## (bn1): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ) ## (conv1): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) ) ## (1): BasicBlock( (bn2): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ## ) (conv2): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) ## ) (relu): ReLU(inplace=True) ## (1): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) (bn1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ## (0): Conv2d(256, 512, kernel_size=(1, 1), stride=(2, 2), bias=False) (conv1): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) ## (downsample): Sequential( (1): BasicBlock( ## (bn2): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ) ## (0): Conv2d(128, 256, kernel_size=(1, 1), stride=(2, 2), bias=False) ## ## V1 V2 V3 (1): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ## ## humain vehicule 4 ) ## ## chasseur humain 3 ) ## ## blaireaux renard 1 (1): BasicBlock( ## ## blaireaux sangliers 1 (conv1): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) ## ## chat renard 1 (bn1): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ## ## chat sangliers 1 (relu): ReLU(inplace=True) ## ## chevreuil chien 1 (conv2): Conv2d(256, 256, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) ## ## chevreuil renard 1 (bn2): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ## ## chien chamois 1 ) ## ## chien sangliers 1 ) ## ## humain chasseur 1 (7): Sequential( ## ## lievre renard 1 (0): BasicBlock( ## ## renard blaireaux 1 (conv1): Conv2d(256, 512, kernel_size=(3, 3), stride=(2, 2), padding=(1, 1), bias=False) ## ## renard sangliers 1 (bn1): BatchNorm2d(512, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True) ## ## sangliers renard 1 (relu): ReLU(inplace=True) ## (conv2): Conv2d(512, 512, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False) ## vehicule humain 1

Table 6 :

 6 Comparison of the predictions vs. ground truth. %>% as.integer() -> index # extract relevant info predicted[i] <-categories[index + 1] # match it with categories } data.frame(truth = c("lynx", "roe deer", "wild boar"), prediction = predicted) %>% kable(caption = "Comparison of the predictions vs. ground truth.") %>% kable_styling()

	truth	prediction
	lynx	chevreuil
	roe deer	chamois
	wild boar sangliers
	unlist()
	for (i in 1:length(fls)){
	result <-learn %>% predict(fls[i]) # make prediction
	result[[3]] %>%
	str_extract("\\d+")

Acknowledgments

We warmly thank Mathieu Massaviol, Remy Dernat and Khalid Belkhir for their help in using GPU machines on the Montpellier Bioinformatics Biodiversity platform, Julien Renoult for helpful discussions, Delphine Dinouart and Chloé Quillard for their precious help in manually tagging the images, and Vincent Miele for having inspired this work, and his help and support along the way. We also thank the staff of the Federations of Hunters from the Jura and Ain counties, hunters who helped to find locations for camera traps and volunteers who contributed in collecting data. Last, we thank Auvergne-Rhône-Alpes Region, Ain and Jura departmental Councils, The French National Federation of Hunters, French Environmental Ministry based in Auvergne-Rhone-Alpes and Bourgogne Franche-Comté Region and the French Office for Biodiversity for funding the Lynx Predator Prey Program. This work was also partly funded by the French National Research Agency (grant ANR-16-CE02-0007).

files) is available from https://github.com/oli viergimenez/computo-deeplearning-occupany-lynx

Session information

## R version 4. 1.0 (2021-