J Villaret-Cazadamont

N Butin

J.-F Martin

M Tremblay-Franco

C Canlet

E Jamin

R Gautier

D Zalko

N J Cabaton

F Bellvert

N Poupin

Fast-Quenching

Comparison of two sampling

Introduction

Aim  Compare two metabolites sampling methods to perform targeted and untargeted metabolomics on adherent cells.

Material and methods

Results

Discussion and Conclusion

Untargeted HRMS

-n=9 -Data were generated by reverse phase liquid chromatography and using positive and negative electrospray ionization for mass spectrometry (QToF). Data were pre-processed with XCMS using W4M (detection, integration, alignment, export table). Data were then filtered using quality control samples. A total of 1099 and 696 features were obtained for positive and negative ionization respectively.

NMR

-n=9 -1H-NMR spectra were obtained on a Bruker 600 MHz spectrometer. Spectra were acquired at 300 K using NOESYPR1D (512 scan). Fixe size bucketing was performed with a 0.02 width from 0 to 10 ppm. Solvants pics were excluded from the final bucket table.

Experimental design Sample preparation and collection

Instrumental analyses Data processing Statistical analysis  Samples obtained with the fast-quenching method showed less variability than samples from the standard protocol (fig. 3, fig. 4 and fig. 6). Regarding the score plots, Biological batchs appear to be less important in the fast-quenching protocol. This is due to the fast stop of metabolism which limits cell stress and prevents the unphysiologycal transformation of metabolites. Cells' metabolites can be observed in more physiological state.  Samples without cells appear to be less polluted with residual metabolites from media with the fast-quenching method (fig 2. and fig. 5.). This observation is very important for targeted metabolomics to prevent overestimation of metabolite concentrations and also for untargeted metabolomics during the blank filtration in order to preserve a maximum of variables. Cell culture on coverglass allows the extraction of metabolites in a new plate with a minimum of residual metabolites from medium and more efficient washing.  Strong separation observed in PCA score plots between both protocols (fig. 4 and fig. 6) is explained by major differences in metabolites extraction (fig. 7). These differences of metabolites recovery could be explained by the different extraction solutions used in this study. As already demonstrated in the literature, each extraction solution is linked with specific advantages and disadvantages as regards metabolites characterization, due to their different physicochemical properties [6]. However, the extraction solution can easily be adapted depending of the targeted metabolites.

Fast-Quenching Standard

  [1] Z.Leon et al. 2013, Electrophoresis, Mammalian cell metabolomics:Experimental design and samplepreparation.[2] H.Bi et al. 2013, Analytical and Bioanalytical Chemistry, Optimization of harvesting, extraction, and analytical protocols for UPLC-ESI-MS-based metabolomic analysis of adherent mammalian cancer cells. [3] S. Dietmair et al. 2010, Analytical Biochemistry, Towards quantitative metabolomics of mammalian cells: Development of a metabolite extraction protocol. [4] N. J. Cabaton et al. 2018, Frontiers in Endocrinology, An Untargeted Metabolomics Approach to Investigate the Metabolic Modulations of HepG2 Cells Exposed to Low Doses of Bisphenol A and 17β-Estradiol. [5] G. Martano et al. 2014, Nature Protocols, Fast sampling method for mammalian cell metabolic analyses using liquid chromatography-mass spectrometry. [6] C. A. Sellick et al. 2010, Metabolomics, Evaluation of extraction processes for intracellular metabolite profiling of mammalian cells: matching extraction approaches to cell type and metabolite targets. * Joran.Villaret-Cazadamont@inra.fr There are numerous ways to perform metabolomics on adherent cells [1]. The main steps for sample preparation and collection are :  Cell Wash  must be efficient to avoid contamination by metabolites from media.  Quenching  metabolism arrest has to be quick and reproductible to minimize unphysiological changes of metabolites (temperature change, cell stress …)  Cell detachment  appeared to be in favor of scrapping due to metabolites leakages with trypsination [2].  Extraction solvent  should maximize metabolites recovery but impossible to observe the entire range of metabolites by a single extraction and analytical method [3dilution standards, absolute concentration of 54 metabolites using the 12 C/ 13 C ratio were obtained : -Amino acids (LC-MS positive ionization Qexactive)  19 variables. -Major carbon metabolites from central metabolism (intermediates of glycolysis, pentose phosphate pathway or tricarboxylic acid cycle) (IC-MS negative ionization LTQ orbitrap velos).  35 variables

Fig. 7 .

 7 Fig. 7. Venn diagram representing metabolites recovery between both protocolsData were obtained in negative mode of ionization. CAMERA was used to associate variables to specific metabolites adducts.  differences in metabolites extraction due to the use of different solvent between both protocols. More metabolites were obtained with the fast-quenching protocol.

Fig. 3 .

 3 Fig. 3. Comparison of coefficients of variation obtained with absolute quantification between both protocols.Coefficients of variation were obtained using the 12 C/ 13 C ratio for energetic metabolites (A) and for amino acids (B). Medians and interquartile range are represented for both protocols. Median were compared using a Mann & Whitney test  lower coefficient of variation for the fast-quenching protocol for energetic metabolites (pvalue=0.0002) (A). No significant differences for amino acids (B).

Fig. 1 Fig. 4 .

 14 Fig. 1 Experimental design, sample collection and intrumental analyses used to compare both extraction protocols.

Fig. 5 .

 5 Fig. 5. Comparison of blank pollution between both protocolesSpectra were cut with a 0.02 ppm width without normalization. Bucket areas were sumed for each spectra. Median were compared using a Mann & Whitney test).  less pollution in blank (pvalue = 0.0188) due to more efficient washing for the fastquenching protocol.

Fig. 6 .

 6 Fig. 6. PCA score plot of mass spectrometry spectra for fast-quenching and standard protocol. Score plots were obtained using a principle component analysis in negative mode of ionization. Log-pareto scaling was perform for multivariate analyses. Biological batchs are represented with the label B1/B2/B3.  less variability for the fast-quenching protocol.

Fig. 2 .

 2 Fig. 2. Comparison of metabolites levels on blank (o) vs cell samples () in both protocols.12 C/ 13 C ratio were used for these dotplots.  less pollution from cell culture media due to more efficient washing in fast-quenching protocol.


Pablo Rodríguez-Mier

Nathalie Poupin

Carlo De Blasio

Laurent Le Cam

Jourdan Fabien

Improving prediction of essential genes using context-specific metabolic network ensembles

Motivation Dataset curation (essential genes under aerobic conditions for S. Cerevisiae)

• Context-specific constraint-based modelling techniques are useful to reconstruct the metabolic network for a particular tissue/condition and make predictions.

• Current methods generate only one network per condition using gene expression, although there is usually no single network that better fits the data.

• We hypothesize that using an ensemble of alternative optimal networks could improve on average the predictions made by a single network.

• To test this hypothesis, we curated a dataset of essential aerobic genes that are directly implicated in the metabolism of the Saccharomyces Cerevisiae.

• Then, we reconstructed thousands of metabolic networks using a modification of iMAT 1 to enumerate alternative solutions based on two MIP cut methods.

• Results using different ensemble approaches to combine the predictions show that ensembles can be used to improve the classification error for essential genes.

Reconstruction & enumeration of metabolic networks

Reconstruction of the yeast aerobic metabolic sub-network and enumeration of alternative optimal sub-networks

Evaluation & conclusions

• We used three different schemes to combine the predictions of the ensembles: or (essential if any metabolic network predicts essential), and (essential if all networks predict essential), and union (essential if the single network generated from the union of the alternative networks predicts essential). • We observed that the "or ensemble" strategy increases the True Positive Rate without increasing too much the False Positive Rate.

• Around 50 essential genes could only be detected using the "or-ensemble" of networks, from which YNL256W (FOL1) was the most common detected essential gene and YFR030W (MET10), YOL064C (MET22), YDR007W (TRP1) among others (see chart below on the right side) are among the least detected genes.

Yeast gene expression data 2 from aerobic condition Base metabolic network (Yeast 6 3 ) Saccharomyces Genome Deletion Project [START_REF] Winzeler | Functional characterization of the S. cerevisiae genome by gene deletion and parallel analysis[END_REF] Yeast 6 metabolic genes Manual curation

  PABLO RODRIGUEZ-MIER a , NATHALIE POUPIN a , CARLO DE BLASIO b , LAURENT LE CAM b , FABIEN JOURDAN a a Toxalim (Research Center in Food Toxicology), Université de Toulouse, INRAE, ENVT, INP-Purpan, UPS, 31300 Toulouse, France b Institut de Recherche en Cancérologie de Montpellier, Campus Val d'Aurelle, 34298 Montpellier, France


Jérôme A Dupuis
email: dupuis@math.ups-tlse.fr

On the Bayesian estimation of species richness and related quantities from quadrat sampling

Keywords: Bayesian estimate, Biodiversity, Missing data, Quadrat sampling, Species richness

We consider the problem of estimating the number of species of a biological community located in a region R divided into J quadrats. Recently, two parametric approaches have been developed which both take into account in a same modeling framework the detectability and the occurence of species in the quadrats. One developed by Dorazio and Royle (2005, J.A.S.A. 100, 389-398) ignores the unsampled part of R (thus also J) and models the occurence of species only in the sampled quadrats. We show that this approach can be used only if J is large, which limits its use in practice, since the value of J can strongly vary from one study to another. The other developed by Dupuis and Joachim (2006, Biometrics 62, 706-712) does not have this limitation but it applies only in the presence of prior information. In this paper, we propose a new approach which extends these two approaches since it can be used in non informative setups, and applies without limitation on J. We develop our approach within a simple model, which assumes that the species population is homogeneous. This constitutes the suitable framework for examining the effect that ignoring J has on the Bayesian estimate of species richness. In particular, a simulation study is undertaken which shows that the approach of Dorazio and Royle generates an error which can be important for small or moderate values of J, when species are spatially rare or hard to detect.

Introduction

The species richness of a community of animals or plants -that is the number of species present within this community -is a basic and fundamental measure of its bio-diversity [START_REF] Huston | Biological diversity[END_REF]. Estimating the species richness (denoted by S) of a biological community located in some specified region, called afterwards R, often relies on quadrat sampling [START_REF] Krebs | Ecological Methodology[END_REF]. Typically, the region R under investigation is first divided into J quadrats, and then a random sample of T quadrats from J is taken. From quadrat sampling data, the knowledge of the biodiversity of the region R can be increased by estimating quantities closely related to S, such as the number S a of species present in any subregion R a of R.

For biological motivations, see e.g. [START_REF] Chao | Estimating the number of shared species in two communities[END_REF].

Nonparametric methods have been proposed for estimating the species richness when quadrat sampling is used. They include the Jacknife and the boostrap estimates; for a complete review see [START_REF] Bunge | Estimating the number of species: A review[END_REF] and [START_REF] Chao | Species richness estimation and applications[END_REF]. Recently, two parametric approaches have been developed which both take into account in a same modeling framework the two biological processes which underlie the data: the one related to occurence of species in the quadrats, and the other related to species detection. One approach has been developed by [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF], the other by [START_REF] Dupuis | Bayesian estimation of species richness from quadrat sampling data in the presence of prior information[END_REF].

Both are preferable to the popular capture recapture methods, which are actually rarely relevant for estimating the species richness of a biological community, as underlined by [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF]. However, both suffer from serious limitations.

We consider that the approach developed in [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF], and called afterwards the DR approach, needs to be discussed on three points.

• First, the number of quadrats in which R has been divided plays no part in the DR approach, while it is expected that the estimation of S depends both on T and J; see e.g. [START_REF] Mingoti | Estimating the total number of distinct species using presence and absence data[END_REF], [START_REF] Hass | An estimator of number of species richness from quadrat sampling[END_REF]. Furthermore, we note that no indication or discussion is provided concerning the use of the DR approach with regard to the sampled fraction of R. One of the objectives of this paper is precisely to clarify this question.

• Second, [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF] model the occurence of species in the sampled quadrats, by T independent Bernoulli outcomes. We show in this paper that this assumption of independence is not consistent with the fact that their approach is conditional, which means that they model the occurence of species effectively present in R.

• Third, their approach is not fully Bayesian but rather ad hoc. Indeed, credible posterior intervals are provided for S, but no prior distribution is actually put on S. By way of contrast, in this paper we pay particular attention to the prior adopted for this parameter. Dupuis andJoachim (2003, 2006) developed a Bayesian parametric model to estimate S. In contrast to [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF], their approach, called afterwards the DJ approach, is unconditional, in that it models the occurence of species liable to be present in the region R. It applies without limitation on J, but it is limited to situations where a list L of species liable to be present in the region R can be drawn up, and requires prior information on the probabilities of presence of each species of L not detected, which significantly limits the use of this approach. In particular, the approach of Dupuis andJoachim (2003, 2006) cannot be used in unknown regions. We note that [START_REF] Dorazio | Estimating species richness and accumulation by modeling species occurence and detectability[END_REF] have also developed an unconditional approach. It presents similarities with the DJ approach. In particular, these authors introduce a supercommunity of species which is supposed to encompass the species population of interest. Its implementation is relatively easy when the region R is known, but is problematic in the opposite case. Furthermore, this approach can be used only if J is large, simply because [START_REF] Dorazio | Estimating species richness and accumulation by modeling species occurence and detectability[END_REF] model the occurence of species present in R as in [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF].

In this paper we propose a new approach which does not require the above list L, and can be used with totally unknown regions. It also provides answers to the reservations formulated regarding the DR approach. As [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF] we adopt a conditional approach, but the key difference is that we model the occurence of species in the J quadrats (not only in the T quadrats). Therefore J and T are included in our model. We show that the DR approach is valid only asymptotically (with respect to J), which significantly limits its use (since the range of possible values of J is particularly large in practice, as illustrated in Section 2.1). In fact, the DR approach can be viewed as a particular case of ours, more precisely as a limiting case (namely J -→ ∞). Since J is finite (even though it might be very large in some studies), using the DR approach thus induces an error which needs to be evaluated.

In this paper we focus on a simple model, called M 0 , which assumes that the species population is homogeneous; it can be viewed as the correct version of the DR model (in its homogeneous version). We emphasize that the methodology developed in this paper also applies to heterogeneous populations. (However, the heteregeneous case requires specific and rather lengthy developments and is beyond the scope of this paper.) The model M 0 , which includes a small number of parameters, is actually a suitable framework for pointing out the limitations of the DR approach, as well as for examining the extent to which the error entailed by the use of the DR model (instead of M 0 ) is liable to affect the Bayesian estimation of S. However, the model M 0 exhibits some statistical difficulties which have to be overcome before undertaking this analysis. They are due to the presence of randomly missing data: a species can be present in R but, for practical reasons, not have been detected. There are two types of statistical difficulty. First, obtaining the Bayesian estimates of S and S a involves computational difficulties, which are overcome by implementing a Markov Chain Monte Carlo (MCMC) algorithm which uses a suitable partitioning of the missing data. Second, some parameters can be non identifiable (as it is often the case in missing data models). In our context this is effectively the case if each sampled quadrat is visited only once; but we prove that all the parameters are identifiable if repeated observations are made in sampled quadrats.

The paper is organized as follows. The experimental protocol and the missing data structure are described in Section 2. Section 3 is devoted to the model M 0 . The approach of [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF] is discussed in Section 4. The links existing between conditional and unconditional approaches are clarified in Section 5. We conclude in Section 6.

2. The experimental protocol and the missing data structure 2.1 The experimental protocol and data description Before describing the protocol, we stress that the region R under investigation has to be of course bounded, or else the quantity S is not well defined, and estimating S is an ill-posed problem (as pointed out by [START_REF] Dorazio | Estimating species richness and accumulation by modeling species occurence and detectability[END_REF]. Quadrat sampling can be done in two distinct variants.

-The first variant is described in details for example in [START_REF] Dorazio | Estimating species richness and accumulation by modeling species occurence and detectability[END_REF], cf the Section Protocol for Sampling Communities. It proceeds as follows. The region R is divided into J spatial units, called quadrats for convenience, though they may have different shapes.

In this paper we assume that these quadrats are of equal area. A sample of T quadrats is then taken, and the sampled quadrats are numbered from 1 to T . The draw is usually performed at random so as to have a sample representative of the whole region R. Finally, an experimenter visits K times each sampled quadrat and records the species detected in each. As in [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF], we assume that K ≥ 2. Detections are typically based on visual or aural recognitions; we assume that species are correctly identified.

-The second variant slightly differs from the one above described. It proceeds as follows:

T quadrats of equal area, say A, are placed (typically at random) in the region R; then, they are explored as above indicated (see e.g. [START_REF] Mingoti | Estimating the total number of distinct species using presence and absence data[END_REF]. Contrary to the previous protocol, the region R is not beforehand divided in quadrats. However, we emphasize that the model M 0 also applies to this procotcol, provided the unsampled part of R can be divided (even virtually) in T * quadrats of area A (which is generally the case, in practice); so we have J = T + T * .

When K = 4 and T = 6, a possible record (or history) for a species s present in R is: y s = (3 0 0 0 4 0). Such a record means that species s has been detected in quadrat 1 during three visits, and detected in quadrat 5 during each visit. Moreover, its presence has not been detected in quadrat 2, 3, 4, 6.

The problem is to estimate the number S of species present in R, from the data formed by the records of species whose the presence has been detected at least once. Sometimes, all the J quadrats are explored (thus T = J). Although the methodology developed in this paper focuses on the case J < T , we emphazise that it also applies to the case T = J, with straightforward changes. This remark is of importance, since when J is not too large, sampling the whole region R is not rare (see below).

The size of J is a key element in the discussion concerning the validity of the DR approach (cf Section 4.3). In practice, the range of J is particularly large, since its value can strongly vary from one survey to another. Let us illustrate this variability through a few examples, by limiting ourselves to ornithological surveys made in the south of France.

In [START_REF] Joachim | Peuplements d'oiseaux nicheurs des forêts pyrénéennes centrales; comparaison du peuplement de différentes essences[END_REF], ten forests have been divided in J ≤ 15 quadrats; in this study the values of J are thus small. The size of J can be moderate, as in [START_REF] Dupuis | Bayesian estimation of species richness from quadrat sampling data in the presence of prior information[END_REF] where the forest of Montech has been divided in J = 40 quadrats; or as in [START_REF] Decamps | The importance for birds of the riparian woodlands within the alluvial corridor of the rive Garonne, s.w. France[END_REF] where the Bouconne forest has been divided in J = 98 quadrats. It can be large, as in [START_REF] Lauga | Modeling the effects of forest fragmentation on certain species of forest-breeding birds[END_REF], where R includes J = 5865 quadrats. The distinction between small, moderate and large J is of course rather arbitrary; it actually rests on some observations made in the framework of our simulation study (cf Section 4.3). The data analysed by [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF], as well as by [START_REF] Dorazio | Estimating species richness and accumulation by modeling species occurence and detectability[END_REF], concern the North American breeding bird survey; it includes around 4000 roadsides, and each route contains 50 quadrats. In this survey, J is thus particularly large.

2.2 A specific missing data structure.

To specify the missing data structure inherent in quadrat sampling data, we view the record y s = (y sj ; j = 1, . . . , T ) as the result of two processes: one is related to the presenceabsence process, and the other is related to the detection process. Such a formalism also allows us to formulate rigorously the biological assumptions made and to introduce, in a natural way, the parameters of biological interest.

-For s = 1, . . . , S and j = 1, . . . , J, we denote by z sj the indicator of presence of species s in quadrat j. The vector (z sj ; j = 1, . . . , J) is denoted by z s .

-For a species s present in a sampled quadrat j, we denote by x sj the number of times that species s has been detected in quadrat j during the K visits. Note that x sj is thus defined only conditionally on z sj = 1; the vector formed by the x sj 's is denoted by x s . Note also that x sj = 0 and y sj = 0 do not have the same meaning (see below).

Missing data can occur in different circumstances. First, when a species s has not been detected in quadrat j, it is clear that z sj is missing; this event covers in fact two exclusive situations: either species s is present in quadrat j but has not been detected, or it is not present in quadrat j (and cannot have been detected). Formally, one has the equivalence (y sj = 0) ⇐⇒ (z sj = 1 and x sj = 0) or (z sj = 0). Conversely, when 1 ≤ k ≤ K, one has (y sj = k) ⇐⇒ (z sj = 1 and x sj = k). Secondly, z sj is missing, when quadrat j is not a part of the sampled quadrats. Thirdly, a species s present in R and undetected is such that the whole vector z s is missing. The set of the missing z sj 's is denoted by z m .

The above observations lead to partitioning z m into four blocks, according to the status of quadrat j (sampled or not) and the status of species s (detected or not); these blocks are denoted by z [i] m ; i = 1, .., 4 and defined as follows. The blocks z [1] m and z [2] m are formed by the missing z sj 's where s designates any detected species, and j denotes either a sampled quadrat (i = 1) or an unsampled quadrat (i = 2). The blocks z [3] m and z [4] m are formed by the missing z sj 's where s designates any undetected species, and j denotes either a sampled quadrat (i = 3) or an unsampled quadrat (i = 4).

3. The homogeneous model: M 0 3.1 Notation, assumptions and parameters.

Notation is basically the one adopted in Dupuis andJoachim (2003, 2006). p(.) denotes a probability mass function, and π(.) the density of any prior distribution. The null vector is denoted by 0, and the vector which has all its components equal to 1, by 1. We denote by Z the set {0, 1} J , and by Z * the set Z\{ 0}. Let v be a vector; we denote by |v| the sum of all its components. 1I (C) denotes an indicator function that takes the value 1 when the condition C is true and zero otherwise. The Bayesian estimate of S is denoted by S 0 .

We denote by z sj the vector z s from which the j-th component has been removed, and by z * s the vector formed by the sampled z sj ; thus z * s = (z sj ; j = 1, . . . , T ). Finally, we set x = (x s ; s = 1, . . . , S) and z = (x s ; s = 1, . . . , S).

We assume that the species population is closed, in that the number of species present in R is constant throughout the quadrat sampling experiment. Biological assumptions related to the occurence of species in the J quadrats are supported by the z s 's, and those related to the detections in the T quadrats by the x s 's.

Assumption A1. We assume that:

p(z) = S s=1 p(z s ) and p(x|z) = S s=1 p(x s |z s ).
Assumption A2. We assume that the probability of detecting species s in quadrat j does not depend on its (possible) detections in the other quadrats.

Assumption A3. We assume that x sj |z sj = 1 ∼ Binomial(K, q). Assumptions A1, A2 and A3 are standard (they are also present in the DR's model).

A1 means that the species present in R do not interact relative to their detectability and their presence (in the quadrats). We thus exclude predator-prey relationships between species. A3 means that, s and j being fixed, the detections of species s during the K visits in quadrat j are independent. Note that q represents the probability of detecting species s in quadrat j during any visit, given that it is present in quadrat j.

Assumption A4. For any species s present in R, we assume that:

p(z s |ϕ) = ϕ |zs| (1 -ϕ) J-|zs| 1 -(1 -ϕ) J (3.1)
where ϕ ∈]0, 1[, and |z s | represents the number of quadrats in which species s is present.

Assumption A4 thus models the occurence of species in the J quadrats; this step of the modeling is absent from the DR approach which model their presence only in the T sampled quadrats (see Section 4.1 for details). This is a key difference between the two approaches.

Note that z s takes its values in Z * since z s cannot take the value 0 (recall that s designates a species present in R) and that p(.) given by (3.1) effectively defines a probability mass function on Z * . We now provide two results concerning the distribution of z s , subsequently useful to establish the distribution of the missing data conditionally on the (observed) data.

Proposition 3.1. Let z s (a) be a vector extracted from z s of length J a . We have: proof. See Appendix A1.

p(z s (a)|ϕ) = ϕ |zs(a)| (1 -ϕ) Ja-|zs(a)| 1 -(1 -ϕ) J if z s (a) = 0, and p(z s (a)) = [(1 -ϕ s ) Ja -(1 -ϕ) J ] / [1 -(1 -ϕ) J ] otherwise. Let z s (a),
By using the proposition (3.1) above, it is easy to verify that pr(z sj = 1|z sj = 0) = ϕ; therefore, ϕ represents the probability that species s is present in quadrat j, given that it is present in at least one other quadrat. Note that the z sj 's are not independent, although a certain form of conditional independence between the z sj 's holds. Indeed, from the proposition (3.1), one deduces that if i, j, l denote three distinct quadrats, z s,i and z s,j are independent, conditionally on z s,l = 1. In other words, the presence of species s in quadrat i does not affect its presence in quadrat j (on condition that it is present in a third quadrat). This biological assumption is standard; it has been discussed in Dupuis andJoachim (2006, 2003).

Note that we do not make the parameters q and ϕ depend on s and j, which constitutes the characteristics of the homogeneous model M 0 . We stress that the model M 0 deals only with species s present in R; the approach adopted in this paper is thus conditional. To avoid needless discussions, we assume that at least one species has been detected by the quadrat sampling, so that the data set is not empty. Concerning S, this implies that we assume that S ≥ 1.

3.2 Likelihood under the model M 0 .

Recall that each record (or history) is described by a vector of length T and having all its components in {0, 1, . . . , K} (cf Section 2.1). Histories are numbered from h = 0 to

H = K T -1;
the history numbered 0 being associated with the record 0. For convenience, the history numbered h is afterwards refered as the history h.

ρ h = T j=1 K b j .
Due to the presence of missing data, writing down the likelihood in a closed form is not immediate, but performing this task is greatly facilitated by the formalism introduced in Section 2.2.

Proposition 3.2. The likelihood of the data y under the model M 0 is:

L 0 (θ; y) = p(y|S, ϕ, q) = S! (S -d)! H h=1 n h ! β S-d 0 H h=1 β n h h (3.2)
where

β h = ρ h q w h (1 -q) Kv h -w h ϕ v h [(1 -q) K ϕ + 1 -ϕ] T -v h 1 -(1 -ϕ) J (3.3)
if h = 0, and

β 0 = [(1 -q) K ϕ + 1 -ϕ] T -(1 -ϕ) J 1 -(1 -ϕ) J . (3.4) proof. See Appendix A2.
We set V = H h=1 v h and W = H h=1 v h . From the above proposition, we deduce that:

L 0 (θ; y) ∝ S! (S -d)! ϕ V q W (1 -q) KV -W F (ϕ, q) dT -V [F (ϕ, q) T -(1 -ϕ) J ] S-d [1 -(1 -ϕ) J ] d
where F (ϕ, q) = (1q) K ϕ + 1ϕ and the constant of proportionality depends only on

y. It is thus clear that Y = (d, V, W
) is a sufficient statistic. This remark is used later, in Section 4.3. Moreover, we emphasize that the expressions (3.3) and (3.4) also hold when T = J (for brevity, details are omitted ).

3.3

The prior distributions.

We assume that S, ϕ, q are a priori independent; thus, π(θ) = π(S)π(ϕ)π(q). Uniform distributions are placed on q and ϕ. Of course, beta distribution can be used if some prior information is available. The negative binomial distribution is usually the distribution adopted in practice for incorporating some prior information on an integer parameter; see, eg [START_REF] King | On the Bayesian estimation of population size[END_REF], in a capture-recapture set-up where the parameter of interest is the size of the animal population. Since the negative binomial distribution can be parametrized by its mean and its variance (cf Appendix 3), adopting such a distribution is thus very flexible; in particular, a large variance for S will reflect a poor prior on S.

When no prior information on an integer parameter is available, the improper Jeffreys prior distribution, is usually adopted: that is π(S) ∝ 1 S in our context. Now, to our knowledge, no motivation exists in the literature concerning the Jeffreys prior, as already pointed out by Kass and [START_REF]The selection of prior distributions by formal rules[END_REF]; these authors simply note that its extends the standard non informative prior put for a real parameter α > 0 (namely, π(α) ∝ 1/α) to the case of an integer parameter. Our motivation for the use of the Jeffreys prior as a non informative prior rests on that it coincides with the limiting case of a negative binomial distribution in which the prior variance tends to ∞ (the prior mean being fixed). Details appear in Appendix A3. The main alternative to the Jeffreys prior is to take π(S) ∝ 1 (eg Casteldine, 1981). In the next Section, we explain why these two non informative priors should give, in practice, very similar estimations of S.

Identifiability issues.

As mentioned in the introduction, the parameter θ of the model M 0 is not identifiable if each sampled quadrat is visited only once (that is when K = 1); for brevity, the proof is omitted. When K ≥ 2, we have established an identifiability result which holds under a very slight restriction.

Theorem 3.3. We assume that S ≥ 2. The parameters S, ϕ and q are all identifiable.

proof. See Appendix A4.

MCMC algorithms for estimating S and related quantities

This Section is organized as follows. At first, we consider the problem of estimating S alone. To obtain E[S|y], two MCMC algorithms are considered: one uses the missing data structure exhibited in Section 2.2, the other uses only the likelihood of the observed data (3.2). The first algorithm is called Alg01, and the second one Algo2. Both are described in the next Section. Next, we show how to obtain the Bayesian estimates of some quantities related to S; only Alg01 is able to do this (with slight modifications according to which quantity we aim to estimate).

Estimating S alone.

We first describe the MCMC algorithm Algo1. It is implemented on (θ, z m ); such a strategy is actually standard in missing data models (see eg [START_REF] Robert | Monte-Carlo Statistical Methods[END_REF].

Algo1 is a Metropolized Gibbs sampling algorithm. The parameter ϕ is up-dated via a Hastings-Metropolis step. We up-date q and ξ = (S, z m ) via Gibbs steps, as follows:

q ∼ π(q|ξ, ϕ, y) and ξ ∼ p(ξ|ϕ, q, y), where ξ = (S, z m ) is simulated as follows:

S ∼ π(S|ϕ, q, y) and z m ∼ p(z m |S, ϕ, q, y).

Before indicating how to implement each step, we provide the expression of the complete data likelihood denoted by L 0 (θ; y, z m ). By using assumptions A1, A2, A3 and A4, and by observing that (y, z m ) and (x, z) provide the same information on θ, it is easy to check that:

L 0 (θ; y, z m ) ∝ ϕ V ′ (1 -ϕ) JS-V ′ [1 -(1 -ϕ) J ] S q W (1 -q) K(V [1] m +V [3] m ) (3.7)
where

V ′ = V + V m , V m = V [1] m + V [2] m + V [3] m + V [4]
m , and V [i] m denotes the sum of the z sj 's belonging to z [i] m .

Updating S. S is simulated according to the distribution of S|q, ϕ, y. Let π(θ|y) denote the density of the posterior distribution. We have π(θ|y) ∝ L 0 (θ; y)π(θ). Taking into account the expression of L 0 (θ; y), we have:

π(S|q, ϕ, y) ∝ S! (S -d)! λ S-d 0 π(S) 1I (S≥d) (3.8)
where λ 0 is given by (3.4). The indicator 1I (S≥d) expresses the fact that, conditionnally on y, S is necessarily greater than d. From (3.8), it is straightforward to deduce that, S -d|q, ϕ, y follows a NegBin(d, 1

-β 0 ) distribution if π(S) ∝ 1/S, a NegBin(d + 1, 1 -β 0 )
distribution if π(S) = 1. Therefore, these two non informative priors should give very close estimates of S, as long as 1 is small compared with d (which is the case in most studies).

Finally, S -d|q, ϕ, y follows a NegBin

(d + r, 1 -(1 -b)β 0 ) if S ∼ NegBin(r, b).
Updating z m . Due to the form of the complete data likelihood we simply have to simulate

the

V [i] m 's. Simulating V [1] m and V [2]
m directly uses the proposition 3.5 proposition 3.4. Conditionally on (θ, y), the two blocks z [1] m and z [2] m are independent.

V [1] m |y, ϕ, q ∼ Binomial(dT -V, γ) where γ = [ϕ(1 -q) K ]/[ϕ(1 -q) K + (1 -ϕ)], and V [2] m |y, ϕ, q ∼ Binomial((J -T )d, ϕ).
For any undetected species s we have:

p(z s |y s = 0, ϕ, q) = (1 -q) K|z * s | ϕ |zs| (1 -ϕ) J-|zs| [(1 -q) K ϕ + 1 -ϕ] T -(1 -ϕ) J
(3.9)

proof. See Appendix A5.

Simulating V [3] m and V [4] m needs special attention, because these two blocks are not, conditionally on (θ, y), independent, due to the constraint z s = 0. Recall also that the species s involved in the simulation of V [3] m and V [4] m are those which have not been detected (that is such that y s = 0), and for which the whole vector z s is missing and has to be simulated according to (3.9). Simulating V [3] m and V [4] m proceeds as follows. Set V [3] m = 0 and V [4] m = 0, and repeat Sd times the two following steps: 1. Simulate z s according to (3.9).

Set V

[3] m = V [3] m + |z * s | and V [4] m = V [4] m + |z • s |, where z • s = (z sj ; j = T + 1, J).
Updating q. Simulate q according to a Beta 1 + W,

1 + K(V [1] m + V [3] m ) .
Updating ϕ. This is done via a Hastings-Metropolis step, as follows. The proposal ϕ ′ is accepted with probability:

min 1, π(ϕ ′ |S, q, z m , y) π(ϕ|S, q, z m , y) × g(ϕ|ϕ ′ ) g(ϕ ′ |ϕ)
where ϕ, q, S, z m represent the current values, and g(.|.) denotes the density of the instrumental distribution (typically an uniform distribution). Note that π(ϕ|S, q, z m , y) ∝ L 0 (θ; y, z m ) which is given by (3.7).

The other algorithm Algo2 does not use the missing data structure of quadrat sampling data; it only uses the observed data likelihood given by (3.2). Parameter S is updated exactly as in Algo1. ϕ and q are updated via Hastings-Metropolis steps. Details are omitted since this is straightforward.

3.5.1 Estimating some quantities related to S.

We denote by M j the number of species present in the sampled quadrat j, and by M the total number of species present in the sampled part of R. More generally, we consider the quantity S a , which represents the number of species present in some subregion R a composed both of sampled and unsampled quadrats. Only the quantities M j (where j designates a sampled quadrat) and M have been considered by [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF], and we have doubts about the validity of their estimates of M j and M. Explanations are postponed until the end of the Section 5.

Estimating M, M j and S a , uses the fact that each of these quantities is a (simple) function of S and z. For example, S a = S s=1 1I (|zs(a)|≥1) where z s (a) denotes the vector (z sj ; j ∈ R a ). Now, recall that, once the data y are available, only a part of z is known, the other part being missing. Thus, conditionally on y, estimating any of these quantities comes down to estimating its missing part. It is the reason why the missing data approach is particularly well suited to estimating M j , M and S a .

• For obtaining M, we use Algo1 without any modification. The Bayesian estimate of M is obtained by applying the ergodic theorem, as follows:

d + 1 L L l=1 h(S (l) , z (l) m ) -→ E[M|y] = M (L -→ +∞)
where (l) indicates the step of the algorithm, and h(S, z m ) = S s=d+1 1I (|z * s |≥1) .

• For obtaining M j , in addition to simulating V [1] m as indicated in the previous Section, we simulate separately each missing z sj ∈ z [1] m (j being fixed); it is done according to the distribution of z sj |y s , that is according to a Bernoulli(γ) distribution, where γ is given by proposition 3.5. The simulation of V [2] , V [3] and V [4] are without change. The Bayesian estimate of S j is then obtained by applying the ergodic theorem:

d j + 1 L L l=1 h(S (l) , z (l) m ) -→ E[M j |y] = M j (L -→ +∞)
where d j denotes the number of species detected in the sampled quadrat j, and

h(S, z m ) = d s=1 1I (z sj =1, y sj =0) + S s=d+1
1I (z sj =1) .

• We denote by J a the number of quadrats (among the J quadrats) located in R a , and by T a the number of sampled quadrats (among the J a ). We denote by y a the part of the data collected in these T a quadrats, and by d a the number of species detected in R a . Note that it will not be correct to perform inference on S a , only on the basis of y a , because y a and y b (where y b = y\y a ) are not marginally independent (they are independent only conditionally on θ). For obtaining the species richness S a of R a , we complete the simulation -Concerning (V [3] m , V [4] m ), set F a = 0 and add to the steps 1 and 2 of the paragraph untitled Updating z m of the Section 3.5.!1, the following step: F a = F a + 1I (|zs(a)|≥1) .

of (V [1] m , V [2] m ) and (V [3] m , V [4] m ), as follows. -Concerning (V [1] m , V [2] m ), set E a =
Note that E a represents the number of species present in R a , not detected in R a , but detected in R b , where R b = R\R a . F a represents the number of species present in R a but undetected. S a is obtained by applying the ergodic theorem:

d a + 1 L L l=1 [E (l) a + F (l) a ] -→ E[S a |y] = S a (L -→ +∞)
4. Discussion of the model of Dorazio and Royle.

This Section is organized as follows. In Section 4.1, we explain why the way of which [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF] model the occurence of species in the quadrats is not consistent with the fact that their approach is conditional. In Section 4.2, we first provide the likelihood of the model of [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF] in its homogeneous version (called afterwards the DR model). Then, we show that the DR approach is valid only asymptotically (with respect to J). In Section 4.3, a simulation study is performed to quantify the error resulting from the use of the DR model instead of M 0 .

We stress that, in [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF], J denotes the number of sampled quadrats, while this quantity is denoted by T in this paper, as in Dupuis andJoachim (2003, 2006).

4.1 The model of Dorazio and Royle and the constraint z s = 0

First, recall that the approach of [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF] is conditional; that is they model the occurence of species present in R. Consequently the constraint z s = 0 applies to their approach. It is of interest to point out that this constraint does not exist when an unconditional approach is adopted since it models the occurence of species liable to be present in R, and it is quite possible that z s = 0, for some s (see Section 6 for details).

Recall also that [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF] do not model the occurence of species s in the J quadrats, but only its occurence in the T sampled quadrats, by assuming that z s,1 , . . . , z s,T are independent outcomes of a Bernoulli random variable of parameter ψ (cf line 6 of the right part of the page 391 of their paper). The J -T remaining quadrats are thus ignored in their modeling, as well as the constraint z s = 0. From the independence assumption of [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF], called afterwards the DR assumption, it follows that:

p(z * s ) = ψ |z * s | (1 -ψ) T -|z * s | , (4.1)
where |z * s | represents the number of sampled quadrats in which species s is present. The dependence of ψ on s, adopted by these authors, is without importance in this discussion, and has been removed for convenience. The dependence of ψ on j, equally adopted by these authors, is considered later. The expression of p(z * s ) given by (4.1) has to be compared with the one calculated under the model M 0 , and deduced from the proposition (3.1):

p(z * s ) = [ϕ |z * s | (1 -ϕ) T -|z * s | ] / [1 -(1 -ϕ) J ] if z * s = 0 (4.2) and [(1 -ϕ) T -(1 -ϕ) J ] / [1 -(1 -ϕ) J ] if z * s = 0. Note that p(z * s ) in (4.
2) depends on J and T (contrary to 4.1). We suggest that the modeling adopted by [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF] is not correct: the problem comes from the fact that the DR assumption is not consistent with the constraint z s = 0. We support this assertion by means of two arguments.

Argument 1. This concerns the case T = J (that is when the whole region R has been sampled). Assuming that z s,1 , . . . , z s,T are independent implies that pr(z s = 0) = (1ψ) J ; hence the contradiction, since it gives a positive probability to the event z s = 0, which is not consistent with the constraint z s = 0. In fact, pr(z s = 0) = 0 implies ψ = 1, that is z s = 1 with probability one, which means that species s is present in all the quadrats of R (with probability one). Now, this situation has been discarded by [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF] as being not tenable from a biological point of view; recall that it is precisely this latter point which motivated their paper (cf the introduction).

Argument 2. The distribution of z s should (of course) not depend on the number T of sampled quadrats. Now, this is not the case when the DR assumption is adopted, as illustrated by a simple example. We assume that J = 5. When all the quadrats have been sampled (that is when T = J = 5), we have pr(z s = 00001) = (1ψ) 4 ψ, by using (4.1). We now asssume that 80% of the quadrats have been sampled (thus T = 4), and that the sampled quadrats are the quadrats numbered from 1 to 4. We can write

pr(z s = 00001) =pr(z s5 = 1|z s1 = z s2 = z s3 = z s4 = 0)×pr(z s1 = z s2 = z s3 = z s4 = 0). Due to the constraint z s = 0 we have pr(z s5 = 1|z s1 = z s2 = z s3 = z s4 = 0) = 1; moreover due to (4.1) we have pr(z s1 = z s2 = z s3 = z s4 = 0) = (1 -ψ) 4 ; hence pr(z s = 00001) = (1 -ψ) 4
which differs from the value obtained when T = 5. Therefore, when one adopts the DR assumption, the distribution of z s depends on T , which is not acceptable. Note that the two values of pr(z s = 00001) are in fact equal if and only if ψ = 1.

If we now make ψ depend on j, that does not modify the conclusions. Details appear in Appendix A6.

The likelihood under the DR model

The model considered by Dorazio and Royle meets the Assumptions A1, A2, and A3.

The model M 0 and the DR model (that is the homogeneous version of the model considered in [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF] actually differ only from the way of modeling the occurence of species in the quadrats. The DR model is parametrized by ω = (S, ψ, q), where ψ denotes the probability that z sj = 1. Establishing the expression of the likelihood under the DR model, denoted afterwards by L DR (ω; y), is straigthforward, contrary to L 0 (θ; y). (The fact of ignoring the unsampled quadrats and of assuming the independance between the sampled quadrats greatly facilitates the calculations.) If we let µ h =pr(y s = h|q, ψ)

L DR (ω; y) = S! (S -d)! H h=1 n h ! µ S-d 0 H h=1 µ n h h where µ h = ρ h q w h (1 -q) Kv h -w h ψ v h [(1 -q) K ψ + 1 -ψ] T -v h if h = 0,
where ρ h has been defined in Section 3.2 and

µ 0 = [(1 -q) K ψ + 1 -ψ] T .
If we make J -→ +∞ (ϕ, q, T and K being fixed), we observe that the likelihood of the DR model and ours coincide. Indeed, (1ψ) J -→ 0 when J tends to +∞, and the two parameters ϕ and ψ are confounded asymptotically. To check this second point, we introduce the parameter ν equal to the probability that z sj = 1 under the model M 0 ; by using the proposition 3.1, we have

ν = ϕ/[1 -(1 -ϕ) J ].
We note that the parameters ψ and ν do not take their values in the same space; indeed, ψ ∈]0, 1[ while ν ∈]1/J, 1[; the verification is immediate. But, when J -→ +∞, the three parameters ϕ, ν and ψ are confounded, since ν = ϕ/[1 -(1ϕ) J ] -→ ϕ and ν and ψ take now their values in the same set. The approach of Dorazio and Royle is thus valid asymptotically with respect to J, and it is therefore expected that, for large J, the Bayesian estimation of S under M 0 and M DR should be close. The question is to know from which value of J, it is effectively the case. This issue is examined in the next Section, using a wide sample of data sets. We now propose another reading of the fact that the DR likelihood and ours coincide when J -→ +∞. Using the DR model, and thus ignoring the number of quadrats -called J * for this remark -in which R has been divided, comes down actually to using the model M 0 with J -→ ∞, while in fact J = J * ; that clearly entails an error which is quantified in the next Section.

4.3 A simulation study.

It is expected that, for some given data y, more or less important differences appear exhibit some examples for which some significant differences exist, and to examine when they occur (with regard to J, T , K, ϕ and q). The differences e 1 and e 2 , both together, may be globally interpreted as the error (in a statistical sense) one makes concerning S, when one uses the DR model, instead of M 0 . Of course, other alternatives could be taken to appreciate the proximity between the posterior distributions of S under M 0 and M DR .

We have chosen to examine the values taken by e 1 and e 2 from average data sets.

To clarify what we mean by average data set, recall that the statistic average data sets, it makes sense to compare the errors to each other (when q, ϕ, J, T , K vary).

Y = (d, V, W ) is sufficient. J,
Throughout this study S is fixed: S = 100. For different values of J, T , K, ϕ and q, we have calculated the corresponding sufficient statistics Y by classical Monte Carlo methods (by simulating 1000 data sets similar to y). The resulting average data sets have been numbered from n = 1 to 18. Then, for each average data set, we have calculated S 0 by implementing the MCMC algorithm Algo2 (based on the observed likelihood L 0 ). deviation σ 0 appearing in Table 1. J being fixed, we note that, not surprisingly, σ 0 decreases with T , K, ϕ and q (the other factors remaining each time fixed). When T , K, ϕ and q are fixed, it of interest to observe that σ 0 increases with J, which is not surprising (though that is -maybe -less intuitive than the previous observations); indeed, as J increaseswhile T remains fixed -the fraction of the unsampled part increases, which creates more and more uncertainty. See the cases: 3 and 7; 9 and 12; 11, 15, 16 and 17.

Table 1 q * , (d, V, W ), S 0 , σ 0 , S DR , σ DR , e 1 , e 2 for different values of J, T, K, q, ϕ (S = 100) n J T K ϕ q q * d V W S 0 σ 0 S DR σ DR e 1 e
• We now comment on the magnitude of the errors e 1 and e 2 , with respect to J, T , K, ϕ and q by distinguishing two cases, according to J ≤ 20 or J ≥ 40.

Case 1: J ≤ 20. We observe that, when q = ϕ = 0.1, the error e 1 is important if J = 10 (whatever T and K) and smaller if J = 20. We also note that, T , K, q and ϕ being fixed, σ DR is markedly larger than σ 0 . Moreover, we observe that the estimation of S produced by the DR model may exhibit surprising behaviors. So, e 1 is not modified (even increases)

when K increases (cf the cases 3 and 4, as well as the cases 7 and 8); e 1 increases when q increases (cf the cases 1 and 3, as well as the cases 5 and 7); and e 1 increases when T increases (cf the cases 7 and 9). Finally, when ϕ = 0.3, we note that the errors e 1 and e 2 are relatively small, even null (see the cases 2 and 6).

Case 2: J ≥ 40. We note that, when ϕ = 0.1 and K = 4, the errors e 1 and e 2 are now small, even very small: see the cases 10 and 11. But, if we now consider for ϕ values smaller than 0.1, as 0.05, the errors may be particular high: see the case 14 (which has to be compared with the case 10). Keeping the value 0.05 for ϕ and increasing J (cf the cases 15,16,17, 18), we observe that S 0 and S DR , as well as σ DR and σ 0 , practically coincide from about J = 200 (cf the case 18). If we now consider species both relatively spatially rare (ϕ = 0.1) and hard to detect (q = 0.1 and K = 2), the values e 1 = 6 and e 2 = 7 (cf the case 12) are not negligible; the case 12 has to be compared with the case 10 where e 1 = e 2 = 2. When J = 40, we have not observed the undesirable behaviors of S DR above mentioned (that is when J = 10, 20).

• To sum up, two striking facts emerge from this simulation study. Firstly, the errors e 1 and e 2 due to the use of the DR model (instead of M 0 ) globally decrease when J increases (as forecasted in Section 4.2); but, the value of J, beyond which S DR and S 0 , as well as σ 0 and σ DR , will practically coincide (T , K, φ and q being fixed) may be rather large.

Furthermore, this value depends on T , K, ϕ, and q (as well as on the true value of S which has been fixed in our study) which prevents from indicating a threshold beyond which the DR approach will yield numerically acceptable results (although Table 1 gives some orders of magnitude). Secondly, we observe that, for small or moderate values of J (≤ 20), the estimations yielded by the DR model may exhibit undesirable behaviors, with respect to K, T and q.

• Our study is limited to homogeneous species populations, but it can reasonably be expected that similar observations will be made in practically homogeneous populations (in the sense that the variability within the ϕ s 's and the q s 's is small); for heterogeneous populations, the situation is difficult to apprehend intuitively, and additional simulation studies -out of the scope of our paper -will be considered in the future.

Conditional and unconditional approaches

We first clarify the links existing between the unconditional approach of Dupuis andJoachim (2003, 2006) and the conditional approach developed in this paper. Then, we briefly comment the unconditional approach of [START_REF] Dorazio | Estimating species richness and accumulation by modeling species occurence and detectability[END_REF].

The DJ approach first assumes that a list L of species liable to be present in R is available. The idea is that in a known region we may have a precise idea of which species are liable to be present in R; this prior information is typically based on previous studies made in R (or encompassing R), or else in regions similar as R. For any species s of the list L we introduce the indicatrice ξ s which takes the value 1 if species s is present in the region R and 0 otherwise. Note that: S = L s=1 ξ s where L denotes the size of L, and that S ≤ L. The probability mass function (pmf) of z s is as follows:

p(z s |λ s , ϕ s ) = λ s ϕ |zs| s (1 -ϕ s ) J-|zs| 1 -(1 -ϕ s ) J 1I (zs = 0) + (1 -λ s )1I (zs= 0) ,
where λ s represents the probability that ξ s = 1. It is clear that p(.) effectively defines a pmf on Z, that p(.) gives a positive probability to the event z s = 0 (namley 1λ s ), and that

p(z s |ξ s = 1, ϕ s ) = ϕ |zs| s (1 -ϕ s ) J-|zs| 1 -(1 -ϕ s ) J .
since z s = 0 is equivalent to ξ s = 1. Note that p(z s |ξ s = 1) defines a pmf on Z * . We emphasize that it coincides with the pmf given in (3.1); in other words, the occurence of any species s present in R is modeled as in the conditional approach. An analogous characteristics, is also present in the unconditional approach developed by [START_REF] Dorazio | Estimating species richness and accumulation by modeling species occurence and detectability[END_REF]; see below. [START_REF] Dorazio | Estimating species richness and accumulation by modeling species occurence and detectability[END_REF] have also developed an unconditional approach. As mentioned in the introduction, it presents similarities with the DJ approach. In particular, a supercommunity of species is introduced, which is supposed to include the species population of interest located in the region R. [START_REF] Dorazio | Estimating species richness and accumulation by modeling species occurence and detectability[END_REF] uses the notion of maximal species list to make this idea more concrete. The size of this supercommunity (afterwards denoted by L) is assumed to be known; it has to be fixed so that L ≥ S, while S is actually unknown (a similar difficulty also exits in the DJ approach). In known regions, one may have some information on the size of S, which facilitates this task. In unkown regions, drawing up the above mentioned maximal species list is not faisable, and fixing L is clearly problematic (even illusive). Furthermore, the strategy consisting in choosing a very high value for L does not actually solve the difficulty, since high values of L imply high computational costs, as stressed by [START_REF] Dorazio | Estimating species richness and accumulation by modeling species occurence and detectability[END_REF]. (A way of overcoming these difficulties would be probably to put a non informative prior distribution on L.) In short, one can consider that this unconditional approach (as it is) suffers from the same limitation as the DJ approach.

Furthermore, we observe that [START_REF] Dorazio | Estimating species richness and accumulation by modeling species occurence and detectability[END_REF] have modeled the occurence of species present in R as in [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF]; cf the Section untitled Modeling heterogeneity in occurence and detection of species. In particular, the independence assumption between the sampled z sj 's appears, page 845, above the equation (1). Consequently, the unconditional approach developed by these authors suffers from the same limitation as the (conditional) DR approach: that is, it can be used only when J is large.

Recall that M j represents the number of species present in quadrat j. As for the estimate of M j calculated by [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF], we have doubts about the validity of the estimate of M j calculated by [START_REF] Dorazio | Estimating species richness and accumulation by modeling species occurence and detectability[END_REF]. In the first paper j designates a sampled quadrat, and an unsampled quadrat in the second one. The problem comes from the fact that these authors ignore the status of s (detected or not) and of j (sampled or not),

when they indicate the distribution of a missing z sj conditional on the data: concerning Bernoulli (ϕ) (ψ with their notation) whatever the status of s and of j, while, on the contrary, the distribution of z sj |y depends on these status (cf the Proposition 3.4). For concision, we do not comment more the unconditional approach of [START_REF] Dorazio | Estimating species richness and accumulation by modeling species occurence and detectability[END_REF].

Conclusion

In this paper we have developed a new conditional approach to estimate the species richness S of a biological community located in a specified region R divided in J quadrats. It takes into account both the occurence and the detectability of species in the quadrats, as does the DR approach. But, contrary to the latter, it models the occurence of species in the J quadrats (not only in the sampled quadrats).

We have shown that using the DR model, and therefore ignoring the unsampled part of R, comes down to use the model M 0 with J -→ ∞. From a practical point of view, this means that the DR approach can be used only for large values of J. Concerning the data set analysed by [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF], the use of the DR approach poses of course no problem, since J is particularly large. For small or moderate values of J (≤ 100 to give an order of magnitude) our simulation study shows that the error resulting from the use of the DR model (instead of M 0 ) may be important; it is especially the case if the population is composed of spatially rare species and/or hard to detect. We note that our conclusions agree, in a way, with those of [START_REF] Mingoti | Estimating the total number of distinct species using presence and absence data[END_REF] who observed that approaches which ignore J, such as the standard Jacknife estimate, can perform poorly. [START_REF] Hass | An estimator of number of species richness from quadrat sampling[END_REF] made similar comments. Our study has thus pointed out some limits of the DR approach. However, it is difficult to provide a general threshold beyond which the DR model could be used (though our study provides some orders of magnitude for some values of T , K and S). Finally, we are of the same opinion as [START_REF] Hass | An estimator of number of species richness from quadrat sampling[END_REF], [START_REF] Mingoti | Estimating the total number of distinct species using presence and absence data[END_REF], when they assert that any statistical analysis should a priori include J in the model (except if J is particularly large, as in the the data set analysed by [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF] Compared now with the DJ approach, the one developed here has the advantage of applying, to both informative as to non-informative settings, while the DJ approach is strictly limited to situations where a list of species liable to be present in R can be drawn up. In other words, on constrat to the DJ approach, ours applies to both known and unknown regions.

Being able to model the occurence of species in the whole region R allows us to estimate quantities of biological interest which were out of reach before, such as the number of species present in any subregion of R. In fact, the MCMC algorithm implemented in this paper, which takes advantage of the missing data structure of quadrat sampling data, allows us to tackle more complex problems, such as estimating the number of species shared by two (or more) distinct subregions of R. This problem, said to be difficult, has been solved by [START_REF] Chao | Estimating the number of shared species in two communities[END_REF], but only when data consist of a sample of individuals (not of a sample of quadrats) and when the region R has been divided into two subregions.

The homogeneous model M 0 is a basic one, which has turned out to be a suitable framework for showing the importance of introducing J in species richness modeling. Of course, subsequent work will have to focus on more complex models. From this point of view, the model M 0 can be taken as a starting point. For example, the methodology developed in this paper easily extends to models including heterogeneity at a species and/or spatial level. Other extensions of biological interest are equally possible. Current models assume that species occupy the quadrats independently, a challenging problem would thus be to develop a model for estimating the species richness of an animal population within which predator-prey relationships exist.

θ 1 = θ 2 . Let us notice, that conditionally on θ 1 and θ 2 , a data set y, in which d species have been detected, has to satisfied d ≤ inf(S 1 , S 2 ). This remark concerns only the fourth data set denoted by y ′′ (see later).

• We first prove that the parameter q is identifiable. For that, we consider two particular data sets. In the first one, called y * , only one species has been detected, and it has been detected at each visit (therefore |y * | = KT ). In the second one, called y * * , again only one species has been detected, and it has been detected at each visit apart from once (therefore |y * * | = KT -1). Using the proposition 3.2, it is easy to check that:

p(y * |θ 1 ) = S 1 (ϕ 1 q K 1 ) T [(ϕ 1 (1 -q 1 ) K + (1 -ϕ 1 )) T -(1 -ϕ 1 ) J ] S 1 -1 [1 -(1 -ϕ 1 ) J ] S 1
and that:

p(y * * |θ 1 ) = S 1 K ϕ T 1 (1q 1 )q KT -1 After simplification, we obtain: q 1 1q 1 = q 2 1q 2 which implies q 1 = q 2 .

1 [ϕ 1 (1 -q 1 ) K + 1 -ϕ 1 ] T -[1 -ϕ 1 ] J S 1 -1 [1 -(1 -ϕ 1 ) J ] S
• We now prove that ϕ is identifiable. We introduce a third data set, called y ′ , in which only one species s has been detected and y sj = K for all the quadrats j, except for one quadrat in which species s has not been detected; otherwise |y ′ | = K(T -1). Now, using the proposition 3.2, we have:

p(y ′ |θ 1 ) = S 1 [ϕ 1 q K 1 ] T -1 [ϕ 1 (1 -q 1 ) K + 1 -ϕ 1 ] [ϕ 1 (1 -q 1 ) K + 1 -ϕ 1 ] T -[1 -ϕ 1 ] J S 1 -1 [1 -(1 -ϕ 1 ) J ] S 1 .
After simplification, we obtain:

p(y * |θ 1 ) p(y ′ |θ 1 ) = ϕ 1 ϕ 1 (1q 1 ) K + 1ϕ 1 .

Using now the hypothesis that p(y * |θ 1 ) = p(y * |θ 2 ) and p(y ′ |θ 1 ) = p(y ′ |θ 2 ) yields

ϕ 1 ϕ 1 (1 -q 1 ) K + 1 -ϕ 1 = ϕ 1 ϕ 2 (1 -q 2 ) K + 1 -ϕ 2 ,
from which we deduce that ϕ 1 = ϕ 2 (since q 1 = q 2 ).

• To prove that S is identifiable, we use a proof by contradiction. We assume that S 2 = S 1 ; for example, S 2 > S 1 . We introduce a new data set, called y ′′ , in which d = S 1 species have been detected, and all these species have been detected during each visit. We have:

p(y ′′ |θ 1 ) = [ϕ 1 q K 1 ] T S 1 [1 -(1 -ϕ 1 ) J ] S 1 . and p(y ′′ |θ 2 ) = S 2 S 1 [ϕ 2 q K 2 ] T S 1 [ϕ 2 (1 -q 2 ) K + 1 -ϕ 2 ] T -[1 -ϕ 2 ] J S 2 -S 1 [1 -(1 -ϕ 2 ) J ] S 2
where the presence of the terms [[ϕ 2 (1q 2 ) K + 1ϕ 2 ] T -[1ϕ 2 ] J ] S 2 -S 1 and S 2 S 1 in p(y ′′ |θ 2 ) is an immediate consequence of our assumption S 2 > S 1 .

Considering that q 1 = q 2 and ϕ 1 = ϕ 2 , it is easy to verify that p(y ′′ |θ 1 ) = p(y ′′ |θ 2 ) implies that

S 2 S 1 = 1 β 0 S 2 -S 1 (1)
where

β 0 = [ϕ(1 -q) K + 1 -ϕ)] T -(1 -ϕ) J 1 -(1 -ϕ) J
in which q = q 1 = q 2 and ϕ = ϕ 1 = ϕ 2 .

Moreover, it is straightforward to check that p(y * |θ 1 ) = p(y * |θ 2 ) implies that

S 2 S 1 = 1 β 0 S 2 -S 1
.

(2)

From ( 1) and ( 2), we deduce that:

S 2 S 1 = S 2 S 1 (3)
We set n 1 = S 1 -1 and n 2 = S 2 -1; note that n 1 and n 2 are such that n 2 > n 1 ≥ 1 since S 2 > S 1 ≥ 2. It is clear that (3) is equivalent to S 2 S 1 = 1. Now, S 2 S 1 is always > 1 whatever n 2 > n 1 ≥ 1 (the verification is immediate); hence the contradiction.

(Note that S 2 S 1 is not > 1 when n 1 = 0, hence the condition S ≥ 2 in the terms of our proposition.) It is clear that starting with S 2 < S 1 leads to the same contradiction.

Therefore, we conclude that our starting assumption (namely S 2 = S 1 ) is false, and that S 1 = S 2 is thus true.

Appendix A5

• Let s be a detected species; first note that z obs s is not empty and is necessarily different from 0. Moreover, given y s , we can partition z s into three parts: z obs s , z mis s (1) = {z mis sj |j ∈ T } and z mis s (2) = {z mis sj |j / ∈ T } where T = {1, . . . , T }. Due to proposition 3.1, z mis s (1) and z mis s (2) are independent conditionally on (y s , θ) (since z obs s = 0). By using now the first part of Assumption A1, we deduce that the two blocks z [1] m and z [2] m are independent conditionally on (y, θ). Using the same arguments, we deduce that the z sj 's of z [2] m are independent conditionally on (y, θ). Thus V [2] m |y, ϕ, q ∼ Bin((J -T )d, ϕ) as the sum of (J -T )d independent Bernoulli (ϕ) r.v.s.

• Using again the same arguments, we deduce that the missing z sj 's of z [1] m are independent conditionally on (y, θ). We now calculate the probability γ = pr(z sj = 1|y s , ϕ, q), where y s is such that y s = 0 and y sj = 0 since s represents a species which has not been detected in quadrat j. By appling the Bayes formula, we have: pr(z sj = 1|y s ) = p(y s |z sj = 1) pr(z sj = 1) p(y s |z sj = 1) pr(z sj = 1) + p(y s |z sj = 0) pr(z sj = 0) where A i denotes the event: z sj = 0 for all j = i. Moreover, due to the constraint z s = 0, we have pr(z si = 1|A i ) = 1; due to the DR assumption of independence, we have pr(A i ) = j =i (1ψ j ), from which we deduce that: p(z [i] s ) = j =i (1ψ j ). We now assume that T = J; due to DR assumption of independence, we have p(z [i] s ) = ψ i j =i (1ψ j ) which differs from the value obtained when T = J -1.

  z s (b), z s (c) denote three vectors extracted from z s and pairwise disjoint. If z s (c) = 0, thus z s (a) and z s (b) are independent, conditionally on z s (c).

  For h = 0 to H = K T -1, we denote by n h the number of species having the history h. Moreover we denote by d the number of detected species; note that n 0 = Sd. Data are denoted by y and formed by the H counts {n 1 , . . . , n H } which clearly constitutes a sufficient statistic. The probability pr(y s = h) is denoted by β h . Let s be any species having h = 0 as its history, we denote by v h the number of quadrats in which species s is detected, and by w h the total number of visits during which it is detected. Let h ∈ {1, . . . , H} and (b 1 , . . . , b j , . . . , b T ) be the history h, thus we set

  0 and repeat dd a times the two following steps: 1. Simulate |z * s (a)| ∼ Binomial(T a , γ) and |z • s (a)| ∼ Binomial(J a -T a , ϕ), where z * s (a) = (z sj (a); j = 1, T a ) and z • s (a) = (z sj (a); j = T a + 1, J a ).2. Set E a = E a + 1I (|zs(a)|≥1) .

  between E[S|y, M 0 ] and E[S|y, M DR ], as well as between σ[S|y, M 0 ] and σ[S|y, M DR ], where σ[S|y, M] = var[S|y, M]. These quantities are afterwards respectively denoted by S 0 , S DR , σ 0 and σ DR (dropping y in this notation, for convenience). The differences | S DR -S 0 | and |σ 0σ DR | are respectively denoted by e 1 and e 2 . Our objective is to

  T , K and θ being fixed, we call an average data set, any data set y having Y = (a, b, c) as sufficient statistic, where a is the integer closest to E[d|θ], b the integer closest to E[V |θ], and c the integer closest to E[W |θ], the expectation being taken with respect to the distribution of y given θ, and under M 0 . Choosing average data sets is, in a way, neutral, compared with arbitrary (or simulated) data sets. Finally, when one uses

  [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF], see Section 3.4, page 393, line 1; and concerning[START_REF] Dorazio | Estimating species richness and accumulation by modeling species occurence and detectability[END_REF], see page 847, Section untitled Predicting species accumulation as a function of species occurence, second paragraph, line 14. More precisely, they consider that z sj |y ∼

  1 . Now, by hypothesis, p(y * |θ 1 ) = p(y * |θ 2 ) and p(y * * |θ 1 ) = p(y * * |θ 2 ) ; hence p(y * |θ 1 ) p(y * * |θ 1 ) = p(y * |θ 2 ) p(y * * |θ 2 ) .

Acknowledgements

The author is grateful to Professor Emmanuelle Cam and to Jean Joachim , for helpful discussions.

Concerning S DR , we have used the same algorithm (except that L 0 has been replaced by L DR ). The convergence of the Markov chain has been diagnosed by using standard techniques. Independent replications of the simulation run for five million iterations (with the first 10% discarded as burn-in) and from different starting points (ϕ (0) , q (0) ) in ]0, 1[×]0, 1[ produced identical results to one unit place.

We mainly focus on small values for ϕ, namely in [0.05, 0.1], that is on populations composed of spatially rare species. It is, in fact, for this type of population that we expect to observe significant differences between S 0 and S DR , since the term (1ϕ) J is as less negligible as ϕ is small (J being fixed). For q, we consider values in [0.1, 0.3], which correspond to species that are relatively hard to detect during a visit. Nevertheless, we emphasize that q * = 1 -(1q) K (which represents the probability of detecting a species s during the quadrat sampling experiment, given that it is present in a quadrat j) can be high even if q is small. For example, q * = 0.19 if q = 0.1 and K = 2, but q * = 0.65 if q = 0.1 and K = 10. These ranges for ϕ and q are not unusual in some animal populations, such as birds (see eg [START_REF] Dorazio | Estimating size and composition of biological communities by modeling occurence of species[END_REF][START_REF] Dupuis | Bayesian estimation of species richness from quadrat sampling data in the presence of prior information[END_REF].

For each fixed value of J, T , K, ϕ and q, we provide in Table 1 below q * , Y, S 0 , S R , σ 0 , σ DR , e 1 and e 2 ; non informative prior distributions having been adopted (the Jeffreys prior for S, and the uniform distribution for q, ϕ and ψ).

• We make two general comments concerning the results. First, we observe that S 0 is close to 100 (for most values of K, J, T , ϕ and q); that is explained by the fact that, for these values of K, J, T , ϕ and q, the bias of S 0 is close to 0; this observation comes from an additional simulation study, not presented in this paper, for brevity. Secondly, we note that S DR is systematically greater than S 0 . That is simply explained by the fact that using the DR model comes down to using the model M 0 with J -→ ∞.

• Let us briefly comment qualitatively on the different values of the posterior standard

Appendix A1

Let z s (a) be a vector extracted from z s of length J a . We partition z s in z s (a) and z s (b); hence p(z s (a)) = zs(b) p(z s (a), z s (b)). By observing that this sum is over all the possible values of z s (b) when z s (a) = 0, and over all the possible values of z s (b), apart from 0, when z s (a) = 0 (due to the constraint z s = 0), it is easy, by using (3.1), to establish that:

Afterwards, the conditioning on ϕ is omitted, for convenience.

Let z s (c) be a vector extracted from z s of length J c and different from 0. Starting from

and by applying the above result to the vectors (z s (a), z s (c)) and z s (c), we have:

] from which we deduce that:

Similarly

Appendix A2

• Due to the assumption A1, the random variables y 1 , . . . , y S are independent (conditionally on θ). They are also identically distributed, since q and ϕ do not depend on s.

Consequently,

(n 1 , . . . , n h , . . . , n H )|S, β ∼ Multinomial(S, β)

where β = (β 1 , . . . , β h , . . . , β H ) and β h =pr(y s = h|ϕ, q); hence (3.2). Afterwards, the conditioning on (ϕ, q) is omitted, for convenience.

which is equal to [ϕ(1q) K + (1ϕ)] T . Finally, we obtain:

Appendix A3

If S ∼ NegBin (r, b), where r ∈]0, +∞[ and b ∈]0, 1[, let us first recall that its probability mass function is such that:

(1)

We now express r and b in terms of E(S) and var(S). This is easily done by using the well known formulae E(S) = r 1-b b and var(S) = r 1-b b 2 , from which we deduce that:

and b = E(S) var(S) .

(2)

If we let var(S) -→ +∞ in (2) it is clear that, for any fixed E(S), b -→ 0 and r -→ 0.

If we now let b and r -→ 0 in the right member of (1), it is easy to verify that it tends to 1/S, since Γ(r + S) -→ Γ(S) = (S -1)! and (1b) S -→ 1.

Appendix A4

Let IN 2 = IN \{0, 1}. Let θ 1 = (S 1 , ϕ 1 , q 1 ) and θ 2 = (S 2 , ϕ 2 , q 2 ) denote any two values

We have to prove the following implication:

where the hyptohesis p(.|θ 1 ) = p(.|θ 2 ) means that: p(y|θ 1 ) = p(y|θ 2 ) whatever the data set y. We can thus choose any data set y (see nevertheless the restriction below); the difficulty being to find suitable and, as far as possible, simple data sets which us allow to prove that where ϕ and q have been omitted in the conditionings (for convenience). Now, as mentioned at the end of Section 3.3, we have: pr(z sj = 1) = ϕ/[1 -(1ϕ] J ]. Moreover, it is easy to verify that: p(y s |z sj = 1) = (1-q) T p(y s |z sj = 0) (recall that y s = 0). Hence, the expression of γ given in proposition 3.5. Moreover, V [1] m |y, ϕ, q ∼ Bin(dT -V, γ) as the sum of dT -V independent Bernoulli (γ) r.v.s.

• Let s be a species that is not detected. We start from: p(z s |y s = 0) = pr(y s = 0|z s )p(z s ) pr(y s = 0) .

where ϕ and q have been omitted in the conditionings (for convenience). If z s = 0 then pr(y s = 0|z s = 0) = 1 else pr(y s = 0|z s ) = (1q) K|z * s | . Note that this formula holds when z s = 0. By replacing p(z s ) by its expression (cf 3.1), and pr(y s = 0) = β 0 by its own (cf 3.4), we obtain the result. estimates of the parameters of the model (in particular S) and those of any related quantities to S (as N j ).

Appendix A6

The Argument 1 still holds when ψ depends on j, due to the following remark which is parameter-free. Whatever j, one has pr(z sj = 1|z s,j = 0) = 1 (due to the constraint z s,j = 0), and pr(z sj = 1|z s,j = 0) = ϕ. Consequently, the distribution of z sj |z s,j depends on z s,j , which implies that the J r.v. z sj 's cannot be independent.

We now examine the Argument 2, and we set pr(z sj = 1) = ψ j . For any i ∈ {1, . . . , J} we introduce the vector z [i] s = (z [i] sj ; j = 1, . . . , J) where z

[i] sj = 1 if j = i and zero otherwise. i being fixed, we assume that one wishes to infer on S from the sample E i which includes all the quadrats, except the quadrat i; therefore, E i includes T = J -1 quadrats. Omitting, for convenience, the conditionings on the ψ j 's, we can write: p(z [i] s ) = pr(z si = 1|A i ) × pr(A i )