Marine Lallemant

Chrystelle Aurélie D'antona

Aude Vidal

Claire Bourtembourg

Malek Toubin

Marie Chehab

Guillaume Vilchez

Rajeev Boiteux

Lionel Ramanah

Pazart

Aurélie D'antona

Chrystelle Vidal

Aude Bourtembourg

Claire Toubin

Malek Chehab

Marie Vilchez

Guillaume Boiteux

Rajeev Ramanah

Lionel Pazart

Didier Riethmuller

Nicolas Mottet
email: n1mottet@chu-besancon.fr

Conservative management versus systematic suture of isolated vaginal or first-degree perineal tears after delivery: A preliminary randomized efficacy trial

Keywords: perineal lacerations, suturing, pain, sexual well-being

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

INTRODUCTION

In the United Kingdom in 2013, the rate of vaginal tears varied with patient age from 30% to 40% [START_REF]NHS Maternity Statistics[END_REF] . In a study published in 2000, Lundquist et al. found 75% of women sustain any perineal or vulvar lesion after vaginal delivery [START_REF] Lundquist | Is it necessary to suture all lacerations after a vaginal delivery?[END_REF] . The viewpoint of health professionals regarding perineal repair has evolved, and women desire perineal protection at birth. To reduce the medicalization of childbirth and adverse outcomes of surgical repair, a new debate concerning whether to suture has emerged. In 1998, Gordon et al. showed a decrease in immediate postpartum pain and dyspareunia at 3 months and 1 year when the skin was not sutured during second-degree perineal tears [START_REF] Gordon | The Ipswich Childbirth Study: 1. A randomised evaluation of two stage postpartum perineal repair leaving the skin unsutured[END_REF][START_REF] Fleming | Does perineal suturing make a difference? The SUNS trial[END_REF] . A Cochrane meta-analysis found only two relevant studies [START_REF] Elharmeel | Surgical repair of spontaneous perineal tears that occur during childbirth versus no intervention[END_REF] , which showed no differences in clinical outcomes up to eight weeks postpartum between surgical and nonsurgical repair of first-and second-degree perineal tears. Although the distribution of perineal tears has changed following decreases in episiotomy rates, no study has attempted to replicate these findings since 2003.

The presence of sutures can reactive edema, and the resorption of suture material may contribute to perineal pain. Digital compression to achieve hemostasis has been proposed as an option for managing perineal tears. We hypothesized that bleeding vessels in superficial tears during delivery could be digitally compressed, avoiding unnecessary sutures.

The first objective of this study was to evaluate the efficacy and safety of conservative management compared with systematic suture in isolated vaginal or firstdegree perineal tears after birth. The secondary objective was to compare sexual well-being at two and six months postpartum.

MATERIALS AND METHODS

Design and setting

We conducted a prospective preliminary efficacy study at a level III maternity unit between March 2017 and September 2018. This study implemented Simon's 2-step plan (interim analysis and final analysis) to test the success rate of conservative management. A randomized controlled group was added to estimate the efficacy of systematic suture. A preliminary and nondefinitive comparison of the standard regimen and experimental regimen was made for secondary criteria. The protocol was approved by the EST II Ethics Committee (reference number 16/521) and was registered as a clinical trial before the initiation of study recruitment.

All methods were performed in accordance with the relevant guidelines and regulations.

This study was registered on clinicaltrials.gov (NCT02870712).

Participants

Women were eligible and recruited into the study during the eighth-month prenatal visit if all the following criteria applied: a primiparous woman, ≥ 18 years of age, singleton pregnancy, and ≥ 32 weeks of gestation. Women were excluded if they did not consent to participate in the study or if any of the following criteria applied: women that did not fulfill all the inclusion criteria, suspicion of a bleeding disorder, anticoagulant treatment during pregnancy, a history of perineal surgery or infibulation, women serving a prison sentence, and woman unable to read or understand French.

Participants were randomized after delivery if the following criteria applied: occurrence of a first-degree perineal tear or an isolated vaginal tear (without muscle involvement), birth at ≥ Women were not randomized if one of the following criteria applied: withdrawal of informed consent, episiotomy, cesarean delivery, assisted vaginal delivery, manual placental delivery and/or uterine exploration, intact perineum, or a second-degree perineal tear or higher.

Randomization and allocation

A computer-generated 1:1 simple randomization table was created using SAS software.

Treatment allocation was written on index cards and concealed in identical, sealed opaque envelopes to ensure balance between the groups. Randomization was performed in blocks of four. After creating the randomization cards, the computer-generated randomization table was deleted. Study allocation envelopes remained sealed until an isolated vaginal or first-degree perineal tear was confirmed by the physician investigator and a second oral consent was obtained.

Study procedures

Eligible women between 32 +0 and 36 +6 weeks of gestation were approached and informed of the study by obstetric physicians or midwives during antenatal visits. The schedules of enrollment and interventions were described verbally to eligible patients, and a written information was also given (Table 1). All participants provided signed consent after completing the informed consent process. When a woman who had previously consented to participate in the study arrived in the delivery room, her wish to continue the study was orally reconfirmed. The quality of sexual activity was assessed using a self-assessment and the French version of the Female Sexual Functional Index (FSFI) questionnaire [START_REF] Rosen | The Female Sexual Function Index (FSFI): a multidimensional self-report instrument for the assessment of female sexual function[END_REF] .

Immediately after childbirth (D0), midwives performed standardized examination of the perineum to define the degree of perineal tear. A physician investigator confirmed the type of tear and confirmed eligibility for randomization. A second oral consent was obtained.

Randomization was performed only after a confirmation of an isolated vaginal or first-degree perineal tear during a spontaneous term, cephalic, vaginal delivery.

In the intervention group (the conservative management group), no suture and no digital compression were performed in the absence of bleeding. If bleeding occurred, the physician investigators or midwives performed a 5-minute digital compression to obtain satisfactory hemostasis without suturing. If bleeding persisted, a suture was performed to obtain appropriate hemostasis with simple interrupted stiches using fast absorbable thread (OptimR 0 and 3/0). In the control group (the systematic suture group), the physician investigators or midwives performed the suture according to the usual practices in the obstetric care unit-i.e., suture with simple interrupted stiches using absorbable suture thread (OptimR 0 and 3/0).

During the 3-day postpartum period, the midwives assessed each morning for analgesic consumption and perineal pain (Table 1). Pain was assessed using a visual analog scale (VAS) presented horizontally to participants [START_REF] Collins | The visual analogue pain intensity scale: what is moderate pain in millimetres?[END_REF] . On the third postpartum day (D3), when leaving the maternity ward, an assessment of perineal pain (VAS score), pain medication consumption and perineal healing was performed. Postpartum healing of the perineum was evaluated using the REEDA (Redness, Edema, Ecchymosis, Discharge and Approximation of skin reparation) scale. A score of 0 to 3 was assigned to the five following components: redness, edema, bruising, oozing, and appreciation of tear healing 8 (Appendix 1). An appointment was set for the tenth (D10) postpartum day (+/-1 day) with one of the doctor investigators, and a self-monitoring logbook was provided, enabling us to record medication consumption and pain scores (VAS score). At D10, the investigator performed a clinical examination to assess the progression of perineal healing using the REEDA scale, identified complications if any, and evaluated the intensity of perineal pain using the VAS score.

The final follow-up visit occurred two months (± 2 weeks) after delivery, coinciding with the systematic postnatal visit. During this consultation, the practitioner collected information on the resumption of sexual intercourse as well as an assessment of perineal pain (VAS score) and healing quality (REEDA scale). A postnatal FSFI questionnaire was given to the patient for a new self-evaluation during this appointment, as well as for use at six months postpartum.

Data collection

The data collected included maternal medical and obstetric history, maternal demographic and clinical characteristics, details about delivery and perineal tearing, intervention procedures, and maternal perceptions at 3 days, 10 days, and 2 months after childbirth.

Data collection was performed in a standardized manner with a paper case report form (CRF). Next, the data were recorded in an electronic research database (CleanWEB ®). An outpatient follow-up booklet recorded pain, any medication taken (analgesic), and notification of any consultations with a practitioner (gynecologist, midwife, general practitioner) for two months.

Outcomes

The primary outcome was the success of the intervention 10 days after delivery defined as pain evaluated by VAS< 3, satisfactory healing defined as a REEDA score ≤ 2 (and at each sub score), and no bleeding or infection.

Secondary outcomes included:

-Success of the intervention on the third postpartum day with satisfactory wound healing defined as a REEDA score ≤ 2, perineal pain with VAS score < 5, and no local complication (bleeding or infection).

-Success of the intervention in the second month postpartum defined as the absence of perineal pain (VAS equal to 0) and resumption of sexual intercourse without dyspareunia (FSFI questionnaire).

-Sexual well-being (FSFI questionnaire) at six months postpartum -Consumption of analgesics on the third day and second month postpartum -Need for medical care during the first two months postpartum We performed a per-protocol analysis to examine our results for consistency after removing protocol violations that may have confounded our results.

Sample size

No former study has assessed the success rate of systematic suturing. Therefore, we arbitrarily decided, according to our professional experience, a success percentage below which conservative management would be rejected as 80% and a success percentage of efficacy of 92%. With a power of 80% and an alpha risk of 5%, 56 women were included in the experimental arm in two stages: a first stage with 22 women and a second stage with 34 additional patients. The size of the reference group (systematic suture) was the same as that of the experimental group (conservative management)-i.e., 56 patients. Therefore, the total number required was 112 patients.

Statistical analysis

The study was planned as an efficacy trial to assess the primary outcome. To verify the absence of perineal non-healing, complications, and pain in the arm where tears were not systematically sutured, an interim analysis was planned before the trial began and performed by the statisticians after the inclusion of 22 women in the conservative management arm and after collection of follow-up information from D10. This analysis focused on the efficacy and safety of conservative management, integrated at the same time: pain assessment (VAS score <3), perineal healing (all REEDA sub scores ≤ 2), and no bleeding or local infection at D3. Of the first 22 women in the conservative management group, according to Simon's plan, the trial could be stopped for ineffectiveness if the number of successes on D10 was ≤ 18. If the preliminary results were effective, 34 additional women were included in this study arm.

Participant and delivery characteristics were compared between the groups using appropriate tests (t test, Wilcoxon rank sum test, χ 2 test, or Fisher exact test). The proportions of success with conservative management were described using a 90% confidence interval. The calculation of the confidence interval considered the sequential nature of the design. The individual components of the primary and secondary outcomes were compared using X 2 and Fisher exact tests. Bonferroni corrections were performed for multiple comparisons. Statistical analyses were performed using SAS® for Windows version 9.4.

RESULTS

Recruitment

During the inclusion period, we performed 2,209 antenatal visits at our center with women between 32 +0 and 36 +6 weeks of gestation (Figure 1). In total, 861 women (39%) met the eligibility criteria, and initial participation consent was obtained for 569 women (66%). After exclusion, 143 women were randomized to the delivery room, with an inclusion rate of 143/170 (84.1%). Overall, 72 patients were allocated to systematic suture, and 71 were allocated to conservative management. For both groups, the initial maternal and neonatal characteristics were comparable (Table 2). A disparity in the distribution of perineal tears between the groups was notable (p=0.02). In the systematic suture group, 28 women (38.89%) and 44 women (61.11%) had a first-degree perineal tear and an isolated vaginal tear, respectively. In the conservative management group, 15 women (21.13%) and 56 women (78.87%) had a first-degree perineal tear and an isolated vaginal tear, respectively.

Overall, 27/143 (18.8%) women were lost to follow-up on the 10 th day postpartum: 17 in the systematic suture group and 10 in the conservative management group. Outcomes were compared between the two groups using a per-protocol analysis: 72 (at first and third day follow-up) and 55 patients (at 10 th day follow-up) in the systematic suture group and 71 (at first and third day follow-up) and 61 patients (at 10 th day follow-up) in the conservative management group. In this last group, 5 (7%) women had no digital compression and no suture, 33 (46.5%) women had digital compression without suture, and 33 (46.5%) women had sutures because of the ineffectiveness of digital compression.

Results of the interim analysis

Among the 22 women examined on the 10 th day postpartum, the success rate at D3 in the conservative management group was 95.45% (n=21). The interim analysis showed that the success rate was higher than 80%. Therefore, the study continued to reach 56 patients in each group.

Primary outcome of the final analysis

The success rate at D10 was 87.8% (90% CI: [70.5-93.54] (42/55)) in the systematic suture group versus 90% (90% CI: [78. 3-93.8] (53/61)) in the conservative management group (Table 3 and4). Pain at D10 was similar in both groups. The REEDA score was significantly higher in the systematic suture group (1.4 versus 0.9; p=0.036). A REEDA score >2 was observed for 8 patients (14%) in the systematic suture group and 4 patients (6%) in the conservative management group (p=0.2).

Secondary outcomes of the final analysis

Perineal pain was significantly higher at D1 in the systematic suture group (2.38 vs 1.69; p= 0.034). Beyond D3, D10 and M2, no significant difference was found in terms of pain.

The success rate at D3 was not significantly different between the two groups (65.28% vs 77.46%; p=NS). The success rate was higher at two months postpartum in the conservative management group (35% vs 19%; p=NS). Thirty-seven of 42 women (88%) had an EVA = 0 at two months postpartum in the conservative management strategy group.

Pain medications were first-level analgesics. No difference was found in their use at D3 (15% vs 14%; p=NS), at D10 (4% vs 3%; p=NS) or at two months postpartum (0% in both groups) between the two groups. No atypical healing progress, infection, or persistent bleeding requiring a medical consultation occurred in either group. FSFI questionnaires were completed for 106 women at inclusion and 82 (57%) at inclusion and 2 months postpartum (39 in the systematic suture group and 43 in the conservative management group). No significant difference was found between the two groups at inclusion or 2 and 6 months postpartum. At inclusion, the total FSFI score was 23.52 (± 9.19) in the conservative management group and 22.35 (± 8.43) in the systematic suture group. At 2 months postpartum, 73 women had resumed sexual intercourse: 35 (64%) and 38 (66%) women in the systematic suture and conservative management groups, respectively. Concerning dyspareunia at 2 months postpartum, 11 of 31 women (35%) in the conservative management group versus 6 of 32 women (19%) in the systematic suture group had almost never or never experienced dyspareunia on penetration (p= NS). By contrast, 10 of 32 women (31%) in the systematic suture group versus 3 of 31 women (10%) in the conservative management group had almost always or always experienced dyspareunia during penetration (p=NS). Four patients (12%) in the systematic suture group versus none in the conservative management group had almost always or always experienced dyspareunia after penetration (p=0.113). At 6 months postpartum, 97% of women in each group had sexual intercourse: 38/39 women in the systematic suturing group and 34/35 women in the conservative management group.

DISCUSSION

With a success rate of 90% at the 10 th -day follow-up, we demonstrated that conservative management of superficial perineal tears after delivery was safe and efficient. Women from the conservative management group had lower pain at the first-day follow-up and lower REEDA scores at the 10 th -day follow-up. Adverse events were not higher with conservative management than with systematic suture. Whether or not the perineal tear was sutured, no local bleeding or infection events were observed at follow-up.

In our study, the success rates of intervention were similar in both groups at 10 days postpartum (pain evaluated by VAS< 3, satisfactory healing defined by a REEDA score ≤ 2 (and at each sub score), and no bleeding or infection), suggesting that conservative management could be a safe option for managing superficial perineal tears. Interestingly, the REEDA score and edema were significantly higher in the systematic suture group at D10. Lower REEDA scores were found in the conservative management group at two months.

These results contrast with those of Lundquist et al. [START_REF] Lundquist | Is it necessary to suture all lacerations after a vaginal delivery?[END_REF] into two groups: non-suture versus systematic suture. They concluded there were no significant differences in the healing process in the early postpartum period up to day three 2 .

The SUNS trial examined the differences in outcomes between primiparous women (n=74) with and without sutures for first-or second-degree tears during spontaneous vaginal birth [START_REF] Fleming | Does perineal suturing make a difference? The SUNS trial[END_REF] .

They reported significantly faster healing associated with better approximation of the wound in the suture group at day ten. However, we noted the lack of measurement standardization for healing in Lundquist et al., while Fleming et al. also used the REEDA score. The inclusion of perineal muscle tears in these studies makes interpretation of the results more complicated.

Perineal pain evaluated using the visual analog scale was significantly higher at day one in the systematic suture group. However, no significant difference was found in terms of pain evaluation beyond D3, D10 and M2. This difference at day one could be explained by the lack of blindness and the patient's satisfaction at not having been sutured. This bias could have faded over time, explaining why a significant difference was not found later in the study. This possibility is consistent with the finding of no significant difference in the use of analgesic medication between the groups. In the SUNS trial, Fleming et al. also found no significant difference in pain between the groups [START_REF] Fleming | Does perineal suturing make a difference? The SUNS trial[END_REF] .

At six months postpartum, 97% of women in each group had sexual intercourse. In our study, 31% of women (n =10) in the systematic suture group had dyspareunia during penetration compared with 10% (n=3) in the conservative management group (p=NS).

Lundquist et al. found similar results: 18% in the suture group versus 8% in the non-suture group [START_REF] Lundquist | Is it necessary to suture all lacerations after a vaginal delivery?[END_REF] . Superficial perineal tears are frequent (more than 50% of primiparous women). After the exclusion of noneligible women after delivery, the randomization rate for our study was 84.1% (143/170). In the SUNS trial, only 74 of 1,314 eligible (6%) women were randomized.

The main reason given for this low randomization rate was nonadherence by the midwives' team to the study protocol. Our high randomization rate can be explained by our using several training sessions and our regular involvement of our team members in research protocols evaluating the prevention of perineal tears with restrictive use of episiotomy and careful delivery of shoulders with active delivery of the anterior arm using Couder's maneuver [START_REF] Forey | Impact of a selective use of episiotomy combined with Couder's maneuver for the perineal protection[END_REF][START_REF] Mottet | Active delivery of the anterior arm and incidence of second-degree perineal tears: a clinical practice evaluation[END_REF] .

Our study showed a higher REEDA score and edema in the systematic suture group at day 10 and a lower REEDA score in the conservative management group at two months postpartum. These results align with our hypothesis that less reactional edema occurs in the absence of sutures. However, no difference was found in perineal pain or the consumption of analgesics after day one. Women reported less dyspareunia at two months postpartum in the conservative management group.

For ethical reasons, we chose not to include second-degree perineal tears in our research protocol because the existing literature concerning the prognosis for more extensive tears are contradictory. The degree of morbidity after muscular tears could directly influence urinary and anorectal functions.

Our study is limited by a lack of blindness to the treatment assignment of patients and the clinical staff that assessed the pain, REEDA scores, and FSFI questionnaire. However, blinding was not possible because the allocated interventions could easily be observed. Despite randomization, isolated vaginal tears were higher and first-degree perineal tears were lower in the conservative management group. Nineteen percent of women were lost to follow-up because the 30-day examination was unusual and considered a burden by women. However, the choice of this postnatal visit seemed essential to carefully monitor healing within a reasonable timeframe and to identify possible complications of each intervention. Its monocentric nature and our policy of the restrictive use of episiotomy limit the external validity of our results. Because of a double examination of the perineal tear by the midwife/doctor pair, we believe that we have limited misidentification of tear degree. Longterm follow-up concerning incontinency and psychological well-being is missing. Finally, a lack of power likely caused the failure to detect meaningful differences between the two groups.

The implementation of our protocol in a single center limits the external validity of the results, which must be confirmed in a larger, multicenter trial. We did not include long-term follow-up, including urinary and fecal incontinences and psychological well-being. Additional studies are required to confirm that digital compression without suturing can be used safely in maternity wards.

This study is currently the largest trial evaluating the efficacy and safety of conservative management compared to systematic suture in cases of superficial perineal tears following spontaneous vaginal delivery. Women easily adhered to the study concept. More than 2,220 pregnant women were screened, and only 5% declined to participate (n= 107) during screening. One percent (n=8) of women declined to participate in the study after delivery.

The number of women included was calculated before the beginning of recruitment (112 patients), and the Simon 2-step plan minimized ineffective interventions. The follow-up was long enough, allowing us to have sufficient hindsight and robust results.

Conclusion

Although systematic suturing remains our usual technique to manage superficial perineal tears, this study suggests that conservative management according to the quality of hemostasis could be a safe and efficient option. Superficial perineal tears may heal spontaneously without affecting the quality and resumption of sexual intercourse. Women often welcome the option of not suturing. We showed that the success of digital compression without suturing is sufficient (≥ 90%) to continue investigation in a controlled randomized multicenter phase 3 study.

  and Fleming et al. (SUNS trial) 4 . Lundquist et al. recruited and randomized 80 women with minor perineal tears (grade I-II)

Figure 1 .

 1 Figure 1. Flow chart n: number of cases, D10: 10 th day postpartum, REEDA: Redness, Edema, Ecchymosis,

Acknowledgments

We thank Grégory Tio for his help with the statical analyses.

Funding

This study is supported by University Hospital of Besançon, APICHU Ref: API/2015/60.

Word account: 3623 words

Funding

This study is supported by University Hospital of Besançon, APICHU Ref: API/2015/60.

Patient consent and approval

Written patient consent was obtained for the research and publication.

Conflict of Interest and Financial Disclosure Statements

The authors report no conflict of interest.

Contributors

All the authors contributed to the conception and design of the study. NM and RR provided the idea for the research or article, created the hypothesis, and wrote the original proposal.


Marie-Pierre Étienne

Éric Parent

Hugues Benoit

Jacques Bernier

Random effects compound Poisson model to represent data with extra zeros

Keywords: EM algorithm, Importance Sampling, Compound Poisson Process, Random Effect Model, zero-inflated Data 2008 MSC: 62F12, 62P12, 62L12, 92D40, 92d50

This paper describes a compound Poisson-based random effects structure for modeling zero-inflated data. Data with large proportion of zeros are found in many fields of applied statistics, for example in ecology when trying to model and predict species counts (discrete data) or abundance distributions (continuous data). Standard methods for modeling such data include mixture and two-part conditional models. Conversely to these methods, the stochastic models proposed here behave coherently with regards to a change of scale, since they mimic the harvesting of a marked Poisson process in the modeling steps. Random effects are used to account for inhomogeneity. In this paper, model design and inference both rely on conditional thinking to understand the links between various layers of quantities : parameters, latent variables including random effects and zero-inflated observations. The potential of these parsimonious hierarchical models for zero-inflated data is exemplified using two marine macroinvertebrate abundance datasets from a large scale scientific bottom-trawl survey. The EM algorithm with a Monte Carlo step based on importance sampling is checked for this model structure on a simulated dataset : it proves to work well for parameter estimation but parameter values matter when re-assessing the actual coverage level of the confidence regions far from the asymptotic conditions.

Introduction

Often data contain a greater number of zero observations than would be predicted using standard, unimodal statistical distributions. This currently happens in ecology (see [START_REF] Martin | Zero tolerance ecology: improving ecological inference by modelling the source of zero observations[END_REF]) when counting species (over-dispersion for discrete data) or recording biomasses (atoms at zero for continuous data). Such data are generally referred to as zero-inflated data and require specialized treatments for statistical analysis [START_REF] Heilbron | Zero-altered and other regression models for count data with added zeros[END_REF]. Common statistical approaches to modeling zero-inflated data make recourse either to mixture models, such as the Dirac function for the occurrence of extra zeros in addition to a standard probability distribution (see for instance [START_REF] Ridout | Models for count data with many zeros[END_REF]), or to two-part conditional models (a presence/absence Bernoulli component and some other distribution for non zero observations given presence such as in [START_REF] Stefansson | Analysis of groundfish survey abundance data: combining the glm and delta approaches[END_REF]). These models are well-known [START_REF] Barry | Generalized additive modelling and zero inflated count data[END_REF] and offer the advantages of separate fits and separate interpretations of each of their components. Parameters are well understood and interpreted as the probability of presence, and the average abundance of biomass if present.

However, a major flaw of those models is their non-additive behavior with regards to variation in within-experiment sampling effort [START_REF] Syrjala | Critique on the use of the delta distribution for the analysis of trawl survey data[END_REF]. Consider for instance the fishing effort measured by the ground surface swept by a bottom-trawl during a scientific survey of benthic marine fauna. If during experiment i, observation Y i is made with some experimental effort corresponding to the harvesting of some area D i and is assumed to stem from a stochastic model with parameters θ(D i ), then the additivity properties of coherence are naturally required: if we consider two (possibly subsequent) independent experiments i and i ′ on the different non overlapping areas D i and D i ′ , we would expect that the random quantity Y i + Y i ′ stems from the same stochastic model with parameters θ(D i ∪ D i ′ ). A compound Poisson distribution is a sum of independent identically distributed random variables in which the number of terms in the sum has a Poisson distribution. Compound Poisson distributions are candidate models purposely tailored to verify the previous desired infinite divisibility property since the class of infinitely divisible distributions coincides with the class of limit distributions of compound Poisson distributions ( [START_REF] Feller | An Introduction to Probability Theory and Its Applications[END_REF], theorem 3 of chapter 27).

Depending on the nature of the term in the random sum, the compound distribution can be discrete or continuous. The construction of such a compound distribution with an exponential random mark for continuous data and with a geometric one for counts is recalled in section 2. This approach is worthwhile for two reasons. The first is parsimony : there is only one parameter for the Poisson distribution plus an additional one for the probability distribution function -pdfof each component of the random sum. Secondly, the compound construction may assist our understanding in cases where the data collection can be interpreted in terms of sampling a latent marked Poisson field. That is to say that the data appear in latent "clumps" that are "harvested" during the experiment, the Poisson parameter being the presence intensity of such clumps. A random variable is used to mimic the quantity (or the number of individuals in the discrete case) independently in each clump. At the upper level of the hierarchy, random effects are added to depict heterogeneous conditions between blocks of experiments.

In section 3, we develop a stochastic version of the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] with a Monte-Carlo step (using importance sampling) for this non Gaussian random effect model with zero-inflated data. Maximum likelihood estimates and the corresponding variance-covariance matrix are derived. The computational task remains rather tractable thanks to simplifying gamma-exponential conjugate properties in the continuous case (and beta-geometric conjugacy in the discrete case).

In section 4, the hierarchical model [START_REF] Wickle | Hierarchical Bayesian spacetime models[END_REF] with compound Poisson distribution for zero-inflated data is exemplified using a real case study with two marine species, urchin and starfish abundance data from a scientific bottom-trawl survey of the southern Gulf of St. Lawrence, Canada. The EM algorithm performs well in obtaining the maximum likelihood estimates of parameters, but for one of the two species we notice some discrepancy between the actual coverage of the confidence intervals and their theoretical levels (as given by the asymptotic normal approximation). Consequently, we further focus on variance covariance matrix estimation in section 5 and investigate via simulation the behavior of coverage level of confidence intervals for various experimental designs, in search of a practical fulfillment of the asymptotic conditions. Finally, we briefly discuss some inferential and practical issues encountered when implementing such hierarchical models for zero-inflated data.

Model construction

We propose a hierarchical construction to represent data with extra zero collected over a non-homogeneous area. The model is divided into two main layers : in the first one, we model the sampling process within a homogeneous sub-area (strata) and in the second layer, we introduce heterogeneity between strata at the top of the hierarchy using random effects. The first subsections detail the hierarchical constructions for continuous data. In the last subsection 2.4, we sketch out an obvious modification to represent count data.

Compound Poisson process to introduce extra zeros

Imagine that data Y are obtained by harvesting an area D and that there are some clumps distributed according to an homogeneous Poisson process : clumps are uniformly distributed with a constant intensity, say µ.

By harvesting an area D, we pick an integer-valued random variable N of clumps. According to Poisson process property N follows a Poisson distribution of parameter µ D. For each clump i the independent random variables X i or marks (with the same probability distribution) represent for instance the possible biomass in each clump to be collected. The final return will consist of the sum over N clumps of the amount contained in each clump. With the convention that Y = 0 if N = 0 , the random sum :

Y = N i=0 X i , (2.1)
is said to follow a compound Poisson distribution. Figure 1 exemplifies a realization of the total amount of a collect (i.e., sum of the marks) in a sampled region D.

The Poisson-based additivity property avoids the drawback of classical models mentioned in the introduction. Generally, D is the area of the sampled area included in R 2 . We assume an homogeneous region µ(D) = µ D, so that the expected number of collected clumps is proportional to the catching effort. The difficulty with the generalization to an inhomogeneous Poisson process lies in the inference step, not in the modeling step. Consequently we used another approach to deal with heterogeneity (see section 2.3). In the following, we mostly omit to index quantities with this catching effort for presentation clarity, explicitly mentioning it only when necessary.

Summary statistics about such compound distribution Y are easily obtained (the characteristic function is given in appendix A) :

E(Y ) = µ D E(X) Var(Y ) = µ D E(X 2 )
Parameter µ rules the occurrence of zero values when assuming P(X = 0) = 0 i.e. that the random mark is non atomic at 0 :

P(Y = 0) = exp (-µ D).

Choice of the random component X for continuous data

For real-valued data with extra zeros, we will concentrate in this paper on the exponential distribution of parameter ρ for component X such that E(X) = ρ -1 , leading to To keep on with an ecological interpretation of the model, assuming that the mark X follows an exponential distribution of parameter ρ, means for the biologist that the probability of finding a large amount of biomass within a clump is exponentially decreasing and that the average quantity in each clump is ρ -1 . When no clump is collected, there occurs a zero for the model Y . We choose the exponential distribution because of parsimony and because of its interesting conjugate property detailed in section 3.1.2.

E(Y ) = µ D ρ and Var(Y ) = 2 µ D ρ 2 .
This compound Poisson distribution was termed law of leaks (LOL) by [START_REF] Bernier | Théorie du renouvellement -application à l'étude statistique des précipitations mensuelles[END_REF], where X represents elementary unobserved leaks occurring at N holes (uniformly located) along a gas pipeline. In summary :

(Y ∼ LOL(µ, ρ)) ⇐⇒    Y = N i=1 X j , N ∼ P(µ), (X 1 , . . . , X N ) i.i.d ∼ E(ρ)    (2.2)
For the discrete case, a similar definition holds with the corresponding geometric distribution for the marks (see section 2.4).

Random effects

Although the previous compound construction could have formally been extended to non-homogeneous Poisson processes, it is easier but still quite realistic to relax the assumption of homogeneity by considering homogeneous blocks (or strata), modeling possible inter-block dispersion using random effects. We consider S blocks ; in a given block s there are I s grouped observations. We denote by Y s = (Y s 1 , . . . , Y s Is ) the random vector in block s and by Y = (Y 1 , . . . , Y S ) the whole vector over the S blocks. The coefficients a and b of the gamma pdf Γ(a, b) for a random variable µ are such that E(µ) = a b and Var(µ) = a b 2 . The random effect model RLOL(a, b, c, d) representing the occurrence of the sample Y is defined by the following set of equations.

Y ∼ RLOL(a, b, c, d) ⇐⇒      (µ 1 , . . . , µ S ) i.i.d ∼ Γ(a, b), (ρ 1 , . . . , ρ S ) i.i.d ∼ Γ(c, d), Y s,1 , . . . , Y s,Is | µ s , ρ s i ∼ LOL (µ s D s,k , ρ s ) ∀s ∈ {1, . . . , S} .
(2.

3) The choice of a gamma distribution for the random effect is motivated by conjugate properties which are useful in the inference of the model. Section 4.1.3 will show that it may also be quite a realistic distribution for some datasets. The hierarchical construction is summed up by the directed acyclic graph (DAG as termed by [START_REF] Spiegelhalter | Computation on Bayesian graphical models (avec discussion)[END_REF]) in Figure 2.

Compound Poisson process for count data

A similar but discrete version to model count data, can be obtained by changing the nature of the random marks of the Poisson process. In this paper, we study a geometric distribution with parameter p = P(X = 1). The core of the model is thus given by the following compound Poisson process with geometric marks :

(Y ∼ DLOL(µ, p)) ⇐⇒    Y = N i=1 X j , N ∼ P(µ), (X 1 , . . . , X N ) i.i.d ∼ G(p)   
To preserve conjugate properties, the gamma distribution for the random effect on the marks is replaced by a beta distribution so that the count data version of the model is given by :

Y ∼ RDLOL(a, b, c, d) ⇐⇒      (µ 1 , . . . , µ S ) i.i.d ∼ Γ(a, b), (p 1 , . . . , p S ) i.i.d ∼ β(c, d), Y s,1 , . . . , Y s,Is | µ s , p s i ∼ DLOL(µ s D s,k , p s ) ∀s ∈ {1, . . . , S} .
(2.4) where DLOL means Discrete version of Law of leaks and RDLOL discrete law of leaks with random effects.

In most of the paper, we will simply state the main results when technical aspects of the proofs are shared between discrete and continuous cases.

Estimation via the EM algorithm with importance sampling

Hierarchical models such as 2.3 or 2.4 cannot be straightforwardly estimated because of the latent variables. The random effects (µ, ρ) and the unknown numbers of clumps N must be integrated out to obtain the likelihood. The likelihood has no closed form and estimators cannot be directly derived. In such a case, a classical strategy is to use Expectation Maximization algorithm ( [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]) to derive max-likelihood estimates. In our case the E step is not analytically accessible. An alternative is to use a stochastic version of this EM algorithm such as Monte-Carlo EM ( MCEM see [START_REF] Mcculloch | Maximum likelihood variance components estimation for binary data[END_REF] or [START_REF] Mcculloch | Maximum likelihood algorithms for generalized linear mixed models[END_REF]) or stochastic approximation of EM (SAEM see [START_REF] Deylon | Convergence of a stochastic approximation version of EM algorithm[END_REF]).

We detail in this section how to implement a MCEM algorithm using Importance sampling to obtain the maximum likelihood estimation and its empirical variance matrix. Similar results concerning count data process are summed up in the last subsection. From this point onwards we will use brackets to denote pdf 's as many conditioning terms will appear in the probabilistic expressions derived from the model fully specified by the set of equations (2.3). The brackets denote either a density or a discrete probability distribution, as in [START_REF] Gelfand | Sampling based approach to calculating marginal densities[END_REF]. Following Bayesian conventions, we will also allow the parameters to appear as conditioning terms (i.e., instead of writing P(X) we will specify [X|a, b, c, d]) so as to help the reader understand which layer of the hierarchical model (2.3) the probability expression refers to (see Fig 2).

Implementation of the MCEM algorithm

In this paper, θ stands for the set of parameters (a, b, c, d) in the RLOL model.

Given the random effects, the data within a block are independent :

L(θ; Y, N, µ, ρ) = S s=1 L s
where L s denotes the complete log-likelihood in block s, i.e. :

L s = L s (θ; Y s , N s , µ s , ρ s ) = i=1 Is ln ([Y s,i |N s,i , ρ s ] [N s,i |µ s ]) + (3.1) ln ([µ s |a, b]) + ln ([ρ s |c, d])
Following [START_REF] Abba | Tools for statistical inference: Methods for the exploration of posterior distributions and likelihood functions[END_REF], the pivotal quantity in the EM algorithm (recalled in appendix D) is the conditional expectation of the complete log-likelihood :

Q(θ, θ ′ ) = E θ ′ (L(θ; Y, N, µ, ρ)|Y)

Maximization step

To maximize Q(θ, θ ′ ) with respect to θ, we focus on the terms that involve θ :

Q(θ, θ ′ ) = C -θ (Y ) + (a -1) × S s=1 E θ ′ ln µ s | Y s + Sa ln b -b S s=1 E θ ′ µ s | Y s -S ln(Γ(a)) + (c -1) × S s=1 E θ ′ ln ρ s | Y s + Sc ln d -d S s=1 E θ ′ ρ s | Y s -S ln(Γ(c)), (3.2)
where C -θ (Y ) denotes a constant which does not depend on θ.

Differentiating with respect to θ, we obtain the set of equations to be satisfied at the maximum argmax θ Q(θ, θ ′ ):

a b = S s=1 E θ ′ µ s | Y s S (3.3) ln a -ψ(a) = ln       S s=1 E θ ′ µ s | Y s S       - S s=1 E θ ′ ln µ s | Y s S (3.4) c d = S s=1 E θ ′ ρ s | Y s S (3.5) ln c -ψ(c) = ln       S s=1 E θ ′ ρ s | Y s S       - S s=1 E θ ′ ln ρ s | Y s S (3.6)
ψ(x) denotes the digamma function defined as the first logarithmic derivative of Γ(x). No analytical expression can be derived for θ as the argument of the maximum of Q(θ, θ ′ ), but a Newton-Raphson algorithm is efficient and easy to implement with a good empirical starting point as indicated in annex B.

Expectation step by conditioning onto the number of clumps

The right-hand side of equations 3.3 to 3.6 involves

E θ ′ µ s | Y s , E θ ′ ln(µ s ) | Y s , E θ ′ ρ s | Y s and E θ ′ ln(ρ s ) | Y s .
To compute these expected values, we will proceed by conditioning onto the hidden number of clumps N. Proposition 3.1 shows that, given N, these four target quantities are simply marginal expectations of the sufficient quantity N s+ , the only necessary function of N that needs to be evaluated within each block s.

In a second step, integration over the number of clumps is performed by recourse to importance sampling within a block s as detailed in proposition 3.

µ s |N, Y, θ ′ ∼ Γ(a ′ + N s+ , b ′ + D s+ ), (3.7)
and The quantities involved in the E step are given by

ρ s |N, Y, θ ′ ∼ Γ(a ′ + N s+ , b ′ + Y s+ ), (3.8
E θ ′ µ s | y s = a ′ + E θ ′ N s+ |y s b ′ + D s+ , (3.9)
E θ ′ ln(µ s ) | y s =E θ ′ ψ(a ′ + N s+ ) y s -ln(b ′ + D s+ ), (3.10)
E θ ′ ρ s | y s = c ′ + E θ ′ N s+ |y s d ′ + Y s+ , (3.11)
E θ ′ ln(ρ s ) | y s =E θ ′ ψ(c ′ + N s+ ) y s -ln(d ′ + Y s+ ). (3.12)
This result merely comes from the conjugacy property between gamma and

Poisson distributions for µ (gamma and exponential distribution concerning ρ).

The moments of gamma and log gamma, beta and log beta distributions are recalled in appendix C.

In order to go one step further into the calculus, we have to perform the integration over N + . Proposition 3.2 gives the distribution of N + |Y + , θ up to a constant.

Subsequently, the integration over N + will make recourse to importance sampling as proposed in [START_REF] Levine | Implementations of the Monte Carlo EM algorithm[END_REF]. This Monte Carlo algorithm is detailed in proposition 3.3.

Proposition 3.2. Assuming Y ∼ RLOL(a, b, c, d) with S strata, and I s records in stratum s, the conditional distribution of N s |θ, y s is given (up to a constant K) by

[N s |θ ′ , y s ] = K Γ(a ′ + N s+ )Γ(c ′ + N s+ ) b ′ + D s+ ) N s+ (d ′ + Y s+ ) N s+ Is i=1, y i >0 y N si si Γ(N si )Γ(N si + 1) Is i=1, y i =0 δ(N si ) (3.13)
To draw a sample according to the rather intricate looking distribution 3.13, an importance sampling based algorithm is detailed in the following proposition for one replicate (often termed particle). In order to obtain a G-sample, this procedure is repeated for each block G times. Proposition 3.3 (Generate one particle in one particular stratum s according to distribution 3.13). A particle g is a vector (N

(g) s+ , N (g) s1 , . . . , N (g)
sIs ) in a particular stratum s. Omitting s to make the reading easier, we may assume with no loss of generality that the first I + terms are non zero and the I -I + followings are the zero ones. The algorithm to generate one particle g runs as follows:

1. Generate N (g) i = 0 wherever y i=0 for i = I -I + + 1, . . . , I. 10

Generate the value of the random sum N (g)

+ according to the importance distribution :

f IS (N + ) ∝ 1 b ′ + D + N + Y + d ′ + Y + N + Γ(a ′ + N + )Γ(c ′ + N + ) I + i=1 Γ y i Y + N + + 1 Γ (N + -I + + 1)
As the one dimensional importance distribution f IS is a quickly decreasing function of N + , its normalizing constant can be easily approximated and a bounded interval is used in practice as the support of N + .

Generate each N (g) i

for i = 1, . . . , I + so that the vector (N

(g) 1 -1, . . . , N (g)
I + -1)
is distributed according to a multinomial distribution M(N (g)

+ -I + , (y 1 /Y + , . . . , y I + /Y + )).
4. Associate to the vector (N

(g) + , N (g) 1 , . . . , N (g)
I ) generated at the previous step, the importance weight :

w (g) = I + i=1 Γ N (g) + y i Y + + 1 Γ(N (g) i + 1)
The proof of this proposition is straightforward from importance sampling theory (see for instance chapter 3 of [START_REF] Robert | Monte Carlo Statistical Methods[END_REF]).

The weighted sample of N + may be used to approximate the expected conditional value defined in equations 3.9 to 3.12. For instance, quantity 3.10 is approximated by :

E θ ′ ln(µ s ) | y s ≈ 1 G g=1 ω (g) G g=1 ω (g) × ψ(a ′ + N (g) s+ ) -ln(b ′ + D s+ ).

Empirical Variance Matrix

This section is devoted to the evaluation of the empirical variance matrix, so as to provide confidence regions. Because of the EM principle, we assume that the algorithm has converged to the maximum likelihood value θ. The empirical Fisher information matrix is then given by proposition 3.4. To explicitly compute this information matrix, we propose to numerically integrate over N thanks to importance sampling as performed for the point estimation step. Technical details are also given in appendix F.

I e (θ) = - ∂ 2 ln [Y|θ] ∂θ i ∂θ j (3.14)
At the maximum likelihood estimator θ, the following equality holds :

I e ( θ, Y) = S     -ψ ′ (â) 1 b 0 0 1 b -â b2 0 0 0 0 -ψ ′ (ĉ) 1 d 0 0 1 d -ĉ d2     + S s=1 (A s + B s ) (3.

15)

with

A s =      E νs (ψ ′ (a ⋆ s )) -1 b ⋆ s 0 0 -1 b ⋆ s Eν s (a ⋆ s ) (b ⋆ s ) 2 0 0 0 0 E νs (ψ ′ (c ⋆ s )) -1 d ⋆ s 0 0 -1 d ⋆ s Eν s (c ⋆ s ) (d ⋆ s ) 2      and B s =       Var νs (ψ(a ⋆ s )) -Covν s (a ⋆ s ,ψ(a ⋆ s )) b ⋆ s Cov νs (ψ(a ⋆ s ), ψ(c ⋆ s )) -Covν s (c ⋆ s ,ψ(a ⋆ s )) d ⋆ s -Covν s (a ⋆ s ,ψ(a ⋆ s ) b ⋆ s Varν s (a ⋆ s ) b ⋆ s 2 -Covν s (a ⋆ s ,ψ(c ⋆ s ) b ⋆ s Covν s (a ⋆ s ,c ⋆ s ) b ⋆ s d ⋆ s Cov νs (ψ(a ⋆ s ), ψ(c ⋆ s )) -Covν s (a ⋆ s ,ψ(c ⋆ s )) b ⋆ s Var νs (ψ(c ⋆ s )) -Covν s (c ⋆ s ,ψ(c ⋆ s )) d ⋆ s -Covν s (c ⋆ s ,ψ(a ⋆ s ) d ⋆ s Covν s (a ⋆ s ,c ⋆ s ) b ⋆ s d ⋆ s -Covν s (c ⋆ s ,ψ(c ⋆ s ) d ⋆ s Varν s (c ⋆ s ) d ⋆ s 2      
,

where a ⋆ s = â + N s+ , b ⋆ s = b + D s+ , c ⋆ s = ĉ + N s+ , d ⋆ s = d + Y s+ and ν s stands for the probability measure of N s+ | θ, Y.
As for the first derivative phase of the EM algorithm detailed in section 3.3, the operations E N| θ,Y and Var N| θ,Y , needed to evaluate A s and B s , can be easily implemented by recourse to the very same Monte-Carlo N + sample that was previously drawn by importance sampling.

Prediction of the random effects

It is of interest to predict the random effects in each stratum, for instance to help illustrate the heterogeneity between units. In a linear mixed model context, the Best Linear Unbiased Estimator is defined by the conditional expectation of the random effect according to the data y and the point estimation. We follow the same avenue of thought and define a predictor of the random effects by the conditional expectation. Using formula 3.9 and 3.11, the random effect predictors are given by :

µ (pred) s = E(µ s |y, θ) = â + E N s+ |y s , θ b + D s+ , (3.16)
and

ρ (pred) s = E(ρ s |y, θ) = ĉ + E N s+ |y s , θ d + Y s+ , (3.17)
The following section aims at highlighting the differences between the continuous case detailed previously and the discrete one.

MCEM algorithm for RDLOL model

Straightforward transposition to the discrete case

The definition of the model designed for the discrete case and called RDLOL model is given by equation 2.4, in this case the pivotal quantity Q(θ, θ ′ ) reads :

Q(θ, θ ′ ) = C -θ (Y ) + (a -1) S s=1 E θ ′ ln µ s | Y s + Sa ln b -b S s=1 E θ ′ µ s | Y s -S ln(Γ(a)) + ln Γ(c + d) Γ(c)Γ(d) + (c -1) S s=1 E θ ′ ln p s | Y s + (d -1) S s=1 E θ ′ ln(1 -p s ) | Y s (3.18)
The equations satisfied at the maximum for (a, b) are again 3.3 and 3.4. Due to the substitution of a gamma pdf into a beta pdf for the random effects governing the geometric discrete marks in the random sum of counts, parameters c and d verify equations 3.19 and 3.20 (equivalent to equations 3.5 and 3.6 in the continuous data model) :

ψ(c + d) -ψ(c) = - S s=1 E θ ′ ln p s | Y s S (3.19) ψ(c) -ψ(d) = S s=1 E θ ′ ln p s 1 -p s Y s S (3.20)
The approach used for the continuous case is reproduced to obtain, in each stratum s the conjugate conditional density of µ, p , so that the analog to propositions 3.1 and 3.2 is :

Proposition 3.5. Assuming Y ∼ RDLOL(θ ′ ) with θ ′ = (a ′ , b ′ , c ′ , d ′ )
, S strata and I s records in stratum s as in 2.4 , then the complete conditional distributions of µ s and p s in one particular stratum s are given by

µ s |N s+ , θ ′ ∼ Γ(a ′ + N s+ , b ′ + D s+ ), (3.21)
and

p s |N s+ , θ ′ ∼ β(c ′ + N s+ , d ′ + Y s+ -N s+ ). (3.22)
Furthermore the conditional distribution function of N s is :

N s |θ ′ , Y ∝     I + i=1 Y si -1 N si -1 D Nsi si N si !       I i = I-I + +1 δ(N si )   Γ(a ′ + N s+ )Γ(N s+ + c ′ )Γ(Y + -N s+ + d ′ ) (b ′ + D s+ ) Ns+ (3.23)
The choice of an efficient importance sampling distribution in the discrete case is not the straightforward adaptation of the continuous gives and a mixture has to be used to obtain an efficient and well behaved algorithm, detailed in appendix H.

The covariance matrix in the discrete case

The covariance matrix in the discrete case benefits from the same conditional independence decompositions and the adaptation of the continuous case is straightforward given the moments of the beta distribution in appendix C; the result is detailed in appendix G. The weighted sample of N + is used to compute the expectations and variance-covariance terms in the matrix components.

Prediction of the random effects

The predictions of the random effects are just given by the conditional expectations. Unsurprisingly, the predictions in the discrete case and in the continuous one look very similar. µ (pred) s is still given by formula 3.16 and

p (pred) s = E(p s |y, θ) = ĉ + E N s+ |y s , θ ĉ + d + Y s+ , (3.24)

Applications

In this section, we apply the EM estimation procedure to two real datasets of ecological interest. We then study the validity of asymptotic assumptions by assessing the coverage level of confidence regions. the hierarchical framework allows some transfer of information between strata, the other data help to predict ρ in these two strata.

Maximum likelihood point estimation

The estimation procedure follows the EM algorithm detailed in appendix D (with a stopping rule when the sixth decimal does not change between iterations) and gives values of θUrch = (â, b, ĉ, d) = (0.997797, 1.05107, 5.05733, 13.0312),

Confidence intervals

Relying on proposition 3.4, the asymptotic covariance matrices are evaluated at those maximum likelihood arguments : Essentially only â and b (resp ĉ and d) are correlated.

    var(â U rch , Y) var( bUrch , Y) var(ĉ U rch , Y) var( dUrch , Y)     =     0.0587 0.1020 1.6804 14.4793     Corr( θUrch , Y) =     1
To evaluate the actual coverage of confidence regions in the present sampling conditions (that may be far from asymptotics), 16000 simulations were launched, assuming the same number of strata and the same number of data points per stratum as the urchin catches (resp. sunflower starfishes) with θ as hypothetic true parameter, thus disregarding possible bias. As a practical working conclusions, Figures 8 and9 show how to correct theoretical asymptotical confidence intervals.

The results are quite different from one dataset to the other.

• On Urchins dataset, to get an actual 90% confidence region, we must expand as far as the asymptotic ellipse corresponding to a 99.964% normal approximation as shown in Figure 8.

• On Sunflower Starfish dataset, things work better and the 94% asymptotical confidence interval is quite a good surrogate for an actual 90% confidence region!

To understand Table is completely unappropriate for our model. The estimation is clearly biased with a shift to the right (verified on simulations not shown here) although we tried to correct bias as proposed in [START_REF] Hesterberg | Unbiasing the bootstrap-bootknife sampling vs. smoothing[END_REF]. The width of confidence intervals are underestimated for both species and does not even contain the θ-value. The hierarchical structure of the model may explain part of this bad behavior of bootstrap method but this would need further investigations not in the scope of this paper. The left column of Table 1 exhibits two different behaviors according to the species considered.

• The asymptotic variance of maximum likelihood parameters under-estimate strongly the true sampling characteristics in the Urchin case. This may be due to the large numbers of zero's for that species: consequently relatively less non zero data remain for the ρ ′ s (inverse of patch abundance) and the estimation of c and d that rule the between units variation of ρ's may become difficult.

• The Sunflower Starfishes case exhibits much better properties regarding the approximation of the covariance matrix. For this species, less zeros data occur and we guess that enough information is made available in the sample to get correct estimations.

Figures 10 and11 present the predictions for the random effects in each stratum.

Validation of the gamma assumption for random effects

We have assumed that the random effects µ and ρ were distributed according to gamma distributions. This choice was essentially made for technical convenience The pp-plot for µ suggests that the gamma distribution is appropriate (Fig 12 ); this is not true of the gamma pp-plot for ρ (Fig 13). First there are only 36 points estimates because 2 strata are empty and ρ's for these strata are not defined.

Second the probability plot does not adjust to a straight 45 degrees line. Looking more closely at four extreme points in the ρ pp-plot, we found that they come from strata with less than two non-zero data points. Excluding these 4 points produces the much more acceptable fit of Figure 14.

Simulations Studies

The previous section showed different behaviors depending on the species : the EM procedure provides rather reliable estimates for the starfish RLOL statistical features but not for the Urchin ones. The purpose of this section is to check the role of the sampling designs. Simulation studies are performed to explore the quality of the EM estimation procedure and to check the actual coverage of the asymptotic variance-covariance matrix approximation. For each simulation, the estimation procedure depicted in section 3 yields one point estimate and one estimation of the asymptotic covariance matrix. Assuming that the asymptotic approximation holds and using a normal approximation, confidence intervals can be given for the true value. As we work within a simulation context, the true value is known and one can compute the actual proportion of samples for which the asymptotic confidence interval covers the true value.

Simulation design

RLOL Results

The simulation study is achieved for two values of parameters θ corresponding to the two applications developped in section 4.1. We choose θ U rchin = (1, 1, 5, 13) and θ Sunstars = (1.9, 1.8, 1.9, 0.9) as true parameter references for the simulations.

We first present a study of the bias and then an investigation of the actual coverage of confidence intervals.

Bias study

We can study the bias by simulation according to the numbers of strata and the number of measure points within strata. Figures 15 and16 present the results for relative bias obtained with 1000 simulations in each configuration. As expected it decreases quickly with the number of strata and only marginal amelioration is obtained as soon as the number of data per stratum becomes reasonable.

Confidence intervals study

Using 1000 simulations in each cell, the empirical proportion of the asymptotic 90% confidence ellipsoids that cover the true value is given in Figures 17 and18.

With 1000 trials in a binomial distribution with probability p of success, a confidence interval for p = 0.90 is approximatively [88%, 92%] : cells from Figures 17 and18 that belongs to that interval have been colored in light grey. Results about confidence intervals strongly depend on the value of θ. The asymptotic approximation seems quite satisfying for θ Sunstars : the asymptotical conditions are quickly fulfilled and the design of the case study seems acceptable. For θ U rchin however, the present design should be strongly re-enforced (up to 40 points per stratum with 36 strata!) before yielding acceptable estimations, and confidence regions based on asymptotical theory are definitely too optimistic. These two sets of parameter recover two very different situations : the larger number of zeros in the Urchin case may render the estimation procedure more difficult than in the Starfish situation. However one should note that the difference is not markedly pronounced : 34% instead of 24%! Such a simulation study shows that the quality of variance covariance matrix estimation used to build an ellipsoid of confidence behaves has to be checked through this simulation approach by instance to verify whether the asymptotic conditions are fulfilled and that the analyst should beware of overconfidence.

Conclusion and Perspectives

The following conclusions have been reached: 3. Independence between the latent features ρ and µ has been a priori assumed for the random effects between units. This absence of prior correlation is quite a stringent hypothesis as we might expect ρ and µ to covary (e.g, low non-zero realized abundance could stem from either a small µ or a large ρ). Working with a gaussian copula for a joint bivariate distribution for the couple (µ, ρ) is a bad remedy, because we would have lost the conjugate properties and increased computational load. To keep partial conjugacy , a better idea is considering the natural extension of the gamma family, but such bivariate distributions are rather restrictive since they can only take into account positive correlation and need that the two marginals share the same shape parameter. However such a model would remain parsimonious with 4 parameters: one is gained to depict correlation and one is lost to depict the marginals'shape. The issue of correlation has been addressed in [START_REF] Ancelet | Exploiter l'approche hiérarchique bayésienne pour la modélisation statistique de structures spatiales[END_REF] who proved via simulation that the correlation between ρ and µ has little bearing on the property we are ultimately trying to predict in practice, i.e. the realized biomass in a tow. Finally, the correlation indicates that the latent variables ρ and µ are model concepts that should themselves not be overinterpreted; they don't actually characterize the true size and number of organism patches.

4. Stochastic EM inferential techniques (with importance sampling for the non explicit expectation steps) require a modest computational effort since the random effects are taken partially conjugate with the compound Poisson distributions. Auxiliary importance distributions can be proposed by careful inspection the structure of the joint distribution of the latent variables and integrating out as much as can analytically be done. Much advantage is taken from conditional independence, especially when computing the Fisher information matrix by re-sampling with the simulated missing data that have been previously generated to evaluate the maximum likelihood estimate.

However, the value of results given here depends on the errors involved with the use of maximum likelihood asymptotic formula on one hand and on the precision of Monte Carlo sampling algorithms on the other hand. Due to the multidimensional nature of the latent variables to be simulated , the variability between several trials of the importance sampling techniques when evaluating the information matrix (and its inverse) can be important enough, especially when few data makes a rather flat likelihood function.

5. Asymptotic errors bounds need to be checked and corrected if necessary. We relied on a simulation study to get a more reliable idea of their ranges. The simulated sets of zero-inflated data show that, in the Starfish case, one can readily trust the confidence intervals based on the information matrix while in the Urchin case, one should beware of being overconfident. The asymptotic conditions may not be encountered rapidly. For the Starfish case study, the design allowed a reasonable estimation of the RLOL model features. For the other species with a 10% higher probability of getting zero values, safisfying precision estimates with 40 strata need at least collecting 40 data points per stratum before the confidence coverage gets reasonably close to its theoretically recommended approximate value. Because 1600 stations represents generally unrealistically large sampling effort for a marine bottom-trawl survey in that Urchin example, statisticians need to inform practitioners (before launching the data collection) about possible underestimation of uncertainty.

6. Covariates for the fixed effect of environmental variable (depth, temperature and habitat type) could be added to the model, potentially enhancing ecological interpretation of the observed patterns in organism abundance and distribution. However, it may bring a lot of additional burden during the inferential computations since many of the conjugate properties would be lost.

For the same reasons, non exchangeable strata (with for instance an intrinsic CAR structure on the top of the hierarchy as described in (author?) [START_REF] Sudipto | Hierarchical Modeling and analysis for spatial data[END_REF]) have not been considered here. Simple (low dimensional) importance sampling should be replaced with brute force Hastings Metropolis techniques [START_REF] Hastings | Monte Carlo sampling methods using markov chains and their applications[END_REF]. In such a context, it may be worthwhile to work on encoding prior knowledge [START_REF] Kadane | Papers on elicitation with discussions[END_REF] into probability distributions and switch the problem into a Bayesian framework [START_REF] Berger | Statistical Decision Theory and Bayesian Analysis[END_REF], relying on ready-made tools such as WinBugs for inference [START_REF] Spiegelhalter | WinBUGS Version 1.3. User Manual[END_REF].

7. In the case study, the random effect models with compound Poisson distribution for the occurrence of zero-inflated data fit the data well and allow transfer of information between strata to help predict in data-poor units. Its hierarchical structure favors discussion between ecologists and statisticians, and helps query its interpretation in term of ecological situations with extra zeros.

This equation exhibits the infinite divisibility property of Y with regards to parameter µ, which offers a nice conceptual interpretation when returning to the marked Poisson process underneath this stochastic construction : the resulting quantity Y is obtained by collecting a random number of primarily (hidden) batches X i distributed at random with intensity µ. Such a conceptual latent process of aggregates would be intuitive for many ecologists. Conversely, one can easily check by writing the logarithm of their characteristic functions, that traditional models for zero-inflated data (think for instance of the delta-gamma model or the Zero-Inflated Poisson model such as [START_REF] Ridout | Models for count data with many zeros[END_REF]) lack of coherence for adapting to a change of the scale in the experiment.

Among the many choices for the probability distribution f of the random mark of the sum, this paper focuses, for parsimony and realism, on the exponential distribution for X (continuous case) that is :

f (x) = ρe -ρx so that f(ω) = ρ
ρ+iω and ĝ(ω) = e -µ( iω ρ+iω ) . For the discrete case, we suggest the corresponding geometric distribution :

f (x) = 1 x>0 × (1 -r) × r x e iwx leading to f (ω) = 1-r 1-re iω and ĝ(ω) = e -µ " r(1-e iω ) 1-re iω «
for the exponential compound Poisson count model.

B. Initialization of the Newton-Raphson algorithm

The main point on Newton-Raphson algorithm consists in choosing a good initial point. In this paper we use this algorithm to find the zero of ln(a)ψ(a) -C = 0 Note that function ψ verifies the following asymptotic series' expansion [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables[END_REF] :

ψ(x) ∼ x→∞ ln(x) - 1 2x - ∞ n=1 B 2n 2n x 2n ∼ x→∞ ln(x) - 1 2x - 1 12 x 2 + 1 120 x 4 + . . .
The convergence is very fast (see Figure 19) so that we choose to initiate Newton-Raphson algorithm with x 0 = 1 2C . Let Z be a random variable with gamma distribution, Z ∼ Γ(s, t). Using laplace transform it is easy to obtain the first moment of ln(Z) :

E e λ ln(Z) = E Z λ = t s Γ(s) +∞ 0 y λ y s-1 e -ty dy = Γ(s + λ) Γ(s)t λ .
Differentiating this equation with respect to λ, we have the expected value of ln(Z) (when λ = 0) and Z ln(Z) (when λ = 1):

∂E Z λ ∂λ λ=0 = E (ln(Z)) = ψ(s) -ln(t), (C.1) and ∂E Z λ ∂λ λ=1 = E (Z ln(Z)) = s t (ψ(s + 1) -ln(t)) . (C.2)
Taking the second order derivative, we show :

∂ 2 E Z λ ∂λ 2 λ=0 = E ln(Z) 2 = ψ ′ (s) + ψ(s) 2 -2 ln(t)ψ(s) + ln(t) 2 . (C.3)
Therefore the variance-covariance matrix between Z and ln(Z) is :

s t 2 1 t 1 t ψ ′ (s)
C.2. First and second moments for the sufficient statistics of the beta pdf Let S be a random variable with beta distribution S ∼ β(s, t).

E e λ ln(S) = Γ(s + t) Γ(s + t + λ) Γ(s + λ) Γ(s)
So that, by first and second differentiation, one gets, (the derivation is quite straightfully performed if working with ln E e λ ln(S) ) :

E (ln(S)) = ψ(s) -ψ(s + t) , E (ln(1 -S)) = ψ(t) -ψ(s + t) E ln(S) 2 = ψ ′ (s) -ψ ′ (s + t) + (ψ(s) -ψ(s + t)) 2
One can extend the properties of characteristic function by considering the function of the two arguments λ and µ

E e λ ln(S)+µ ln(1-S) = Γ(s + t) Γ(s + t + λ) Γ(s + λ) Γ(s) Γ(t + µ) Γ(t)
By cross-differentiation under regularity conditions (working with ln E S λ (1 -S) µ makes things easier here also) , the joint moment can be analytically obtained :

∂ 2 E S λ (1 -S) µ ∂λ∂µ λ=0,µ=0 = E (ln(S) ln(1 -S)) = -ψ ′ (s + t) + E (ln(S)) E (ln(1 -S))
Therefore the variance-covariance matrix between ln(S) and ln(1 -S) reads :

ψ ′ (s) -ψ ′ (s + t) -ψ ′ (s + t) -ψ ′ (s + t) ψ ′ (t) -ψ ′ (s + t)

D. EM algorithm principle

From a constructive point of view, one often writes

[x |θ ] = z ln[x, z |θ ] [z |θ ′ , x ] dz - z ln [z |θ, x] [z |θ ′ , x] dz = Q(θ, θ ′ ) -H(θ, θ ′ ) (D.2)
The maximum of θ → H(θ, θ ′ ) is achieved in θ = θ ′ [START_REF] Tanner | Tools for Statistical Inference : Observed Data and Data Augmentation Methods[END_REF].

So H(θ, θ ′ ) < H(θ ′ , θ ′ ).Let us consider D.2 for θ and θ

′ ln [x |θ ] -ln [x |θ ′ ] = (Q(θ, θ ′ ) -Q(θ ′ , θ ′ )) + (H(θ ′ , θ ′ ) -H(θ, θ ′ ))
EM algorithm is based upon an iterative procedure which exhibits θ such that

Q(θ, θ ′ ) > Q(θ ′ , θ ′ ) .
The best θ is obtained by

θ = argmax θ Q(θ, θ ′ )
During iteration we can monitor the value of the gradient for the log likelihood :

∂ ln [x |θ ] ∂θ = ∂ ln[x, z |θ ] ∂θ - ∂ ln [z |θ, x] ∂θ (D.3)
Integrating the right hand term with respect to conditional density [z |θ, x ] ,and keeping in mind that, for any sufficiently regular pdf f (z; θ) of variable z with parameter θ one can write:

z ∂ ln f (z,θ) ∂θ f (z; θ)dz = ∂ ∂θ z ∂ ln f (z,θ) ∂θ f (z; θ)dz = 0, we have ∂ ln [x |θ ] ∂θ = z ∂ ln[x, z |θ ] ∂θ [z |θ, x] dz - z ∂ ln [z |θ, x ] ∂θ [z |θ, x ] dz ∂ ln [x |θ ] ∂θ = z ∂ ln[x, z |θ ] ∂θ [z |θ, x] dz (D.4)
We may use this equality (computed by Monte Carlo method) to perform a gradient method to obtain the maximum likelihood or just to check along the iterations that the gradient is going to zero.

D.2. Score function

From now on, let's call Sc(θ, z, x) = ∂ ln[x,z|θ ] ∂θ the score, i.e the complete loglikelihood gradient and Sc(θ i , z, x) = ∂ ln [x,z|θ ] ∂θ i its i th component. ∇θ, equation D.4 proves that its conditional expectation (with respect to [z |θ, x]) is always equal to the likelihood gradient. Pushing the derivation game one step further leads to:

∂ ∂θ j ∂ ln [x |θ ] ∂θ i = z ∂Sc i ∂θ j [z |θ, x] + Sc i ∂ [z |θ, x] ∂θ j [z |θ, x] [z |θ, x] dz ∂ 2 ln [x |θ ] ∂θ i ∂θ j = z ∂ 2 ln [x, z |θ ] ∂θ i ∂θ j + Sc i Sc j - ∂ ln [x |θ ] ∂θ j [z |θ, x] dz (D.5)

D.3. Information matrix

To obtain the covariance matrix of the estimators at the maximum of likelihood, the empirical information matrix needs to be computed. The second order derivative is obtained by differentiating D.1:

∂ 2 ln [x |θ ] ∂θ i ∂θ j = ∂ 2 ln[x, z |θ ] ∂θ i ∂θ j - ∂ 2 ln [z |θ, x] ∂θ i ∂θ j (D.6)
At the maximum θ = θ, formula D.3 implies ∂ ln [x|θ ] ∂θ jj = 0 so that equation D.5 takes a more friendly aspect because the score term ∂ ln [x|θ ] ∂θ j in the right hand side vanishes at θ = θ . Equation D.6 becomes therefore much more handy because it only involves conditional expectations of first and second derivatives of the complete likelihood terms : Since we detail the computation for one particular s, we will omit to mention it in order to make the reading easier. We also note respectively y, D and N the vectors of data, catching efforts and corresponding number of clumps in one stratum. We define J as J(N, ρ, µ) = ρ, µ, N a, b, c, d, y, D .

∂ 2 ln( x θ ) ∂θ i ∂θ j = z   ∂ 2 ln x, z θ ∂θ i ∂θ j + ∂ ln x, z θ ∂θ i ∂ ln x, z θ ∂θ j   z θ, x dz
(E.1)

Then J satisfies the following set of equations :

∝ y, ρ, µ, N |a, b, c, d, D ∝ I i=1 [y i |N i , ρ] [N i |µ, D i ] [µ |a, b] [ρ |c, d ] ∝ I i=1 [y i |N i , ρ, µ] [N i |ρ, µ] µ a-1 e -µb ρ c-1 e -ρd
with the convention that [A|B] ∝ f (A, B) means that the coefficient of proportionality only depends on B. We note I ⋆ the number of zero value y and we reorder the vector y so that the I + = I -I * non zero y i are the first, so that J may be written as :

J(N , ρ, µ) ∝ I-I ⋆ i=1 y N i i e -ρy i ρ N i Γ(N i ) e -µD i (µ D i ) N i Γ(N i + 1) I i ⋆ =I-I ⋆ +1 δ(N i * )e -µD i ⋆ µ a-1 e -µb ρ c-1 e -ρd
Defining Y + = I i=1 y i , N + = I i=1 N i and D + = I i=1 D i , we obtain :

J(N , ρ, µ) ∝ I + i=1 y N i i Γ(N i )Γ(N i + 1) e -ρ(Y + +d) ρ N + +c-1 e -µ(D + +b) µ N + +a-1
Conditionally to the latent vector N , the random effects ρ and µ are independent.

Isolating the terms which depend on µ on one side and those depend on ρ on the other, we find that Then

E θ ′ µ s | y s = E θ ′ a ′ + N s+ b ′ + D s+ y = a ′ + E θ ′ N s+ |y s b ′ + D s+ .
If Z follows gamma distribution Γ(s, t), then E(ln(Z)) = ψ(s)ln(t) (see annex C), so that

E θ ′ ln(µ s ) | y s = E θ ′ ψ(a ′ + N s+ ) y s -ln(b ′ + D s+ ).
We have respectively for ρ s

E θ ′ ρ s | y s = c ′ + E θ ′ N s+ |y s d ′ + Y s+ , and
E θ ′ ln(ρ s ) | y s = E θ ′ ψ(c ′ + N s+ ) y s -ln(d ′ + Y s+ ).

E.2. Proof of proposition 3.2

Let us define J as the distribution of ρ, µ, N θ ′ , y, D in one particular stratum s. We will write J in a bottom-up perspective and consider the distribution of µ and ρ conditionned by N, because µ and ρ are conditionally independant.

J is given by :

J(N, ρ, µ) = ρ, µ, N θ ′ , y, D = ρ N , θ, y µ N, θ ′ , y, D N , θ, y, D
Using the independent conditional gamma distributions of µ and ρ and integrating according to µ and ρ given N, we can exhibit all the terms depending on

N . ρ µ J(N, ρ, µ)dµ dρ = N θ, y, D ∝ I + i=1 y N i i Γ(N i )Γ(N i + 1) I i ⋆ =I-I ⋆+1 δ(N i ⋆ ) (b ′ + D + ) N + (d + Y + ) N + Γ(a + N + )Γ(c + N + ) -1 E.3. Proof of proposition 3.4 606
In the following Z will stand for all the hidden variables i.e Z = (N, µ, ρ) , |M ij | is another notation for matrix M that details the content of the i th row and j th column, and ∂F (θ) ∂θ stands for the gradient of F written as a vector whose i th component is the scalar ∂F (θ) ∂θ i . The key equation involves rewriting equation D.7 as the expectation of the second order derivative of the complete log-likelihood and the variance of the score (its gradient) to be taken with regards to the conditional distribution Z x, θ (see annex D.3)

∂ 2 ln( x θ ) ∂θ i ∂θ j = E Z|x ∂ 2 ln( x, Z θ ) ∂θ i ∂θ j + Var Z|x ∂ ln [x, Z |θ ] ∂θ (E.
∂θ i ∂θ j = ∂ 2 ln( x, z θ ) ∂θ i ∂θ j = S     -ψ ′ (â) 1 b 0 0 1 b -â b2 0 0 0 0 -ψ ′ (ĉ) 1 d 0 0 1 d -ĉ d2    
As shown in Figure 20. , given Y s , Y s ′ and θ, the latent variables Z s and Z s ′ of two stratum s and s ′ are conditionnaly independent, therefore : To evaluate the variance of the score in stratum s, we will take advantage of successive conditioning due to the hierarchical structure depicted in Figure 2.

Var Z|x ∂ ln [x, Z |θ ] ∂θ = S s=1 Var Zs|x     ln(µ s ) -µ s ln(ρ s ) -ρ s    
Recalling that the latent variable Z s includes, in addition to (µ s , ρ s ), the vector N s , i-e the latent number of clumps for each record, the variance conditional decomposition formula gives:

Var Zs|x     ln(µ s ) -µ s ln(ρ s ) -ρ s     = E Ns|x     Var     ln(µ s ) -µ s ln(ρ s ) -ρ s |N s         +Var Ns|x     E     ln(µ s ) -µ s ln(ρ s ) -ρ s |N s        
So that we have

I e ( θ, x) = S       -ψ ′ (â) 1 b 0 0 1 b -â b2 0 0 0 0 -ψ ′ (ĉ) 1 d 0 0 1 d -ĉ d2       + S s=1 (A s + B s ) with A s = E N|x     Var     ln(µ s ) -µ s ln(ρ s ) -ρ s |N         and B s = Var Ns|x     E     ln(µ s ) -µ s ln(ρ s ) -ρ s |N s         .
Given

′, c s ′, d s ′) = (a + N s+ , b + D s+ , c + N s+ , d + Y s+ ) such that : E     ln(µ s ) -µ s ln(ρ s ) -ρ s |N     =      ψ(a ′ s ) -ln(b ′ s ) -a ′ Var     ln(µ s ) -µ s ln(ρ s ) -ρ s |N     =     -ψ ′ (a ′ ) 1 b ′ 0 0 1 b ′ -a ′ b ′2 0 0 0 0 -ψ ′ (c ′ ) 1 d ′ 0 0 1 d ′ -c ′ d ′2    
and the expression for A s follows easily. In the first derivative, the latent variables µ and ρ appear not surprisingly only through their arithmetic or geometric means (sufficient statistics for the gamma pdf ). Using standard notation μ for the arithmetic mean 1 S S s=1 µ s , we have :

∂ ln [x, z |θ ] ∂a = S ln(µ) + ln b -ψ(a) ∂ ln [x, z |θ ] ∂c = S ln(ρ) + ln d -ψ(c) ∂ ln [x, z |θ ] ∂b = S a b -µ ∂ ln [x, z |θ ] ∂d = S c d -ρ
The gradient of the complete log-likelihood (so-called the "score") may be split into two parts : the first one ∆ θ does not depend on the latent variable z while the other one ∆ z gathers terms depending on z (and possibly of θ), i.e :

  ∂ ln x, z θ ∂θ   = ∆ θ + ∆ z with ∆ θ = S     ln b -ψ(a) a b ln d -ψ(c) c d     ∆ z = S     ln(µ) -µ ln(ρ) -ρ    
In addition here, ∆ z does not contain terms with θ, consequently the second order derivatives are easy to obtain and don't involve the latent variable :

∂ 2 ln [x, z |θ ] ∂a∂a = -Sψ ′ (a) ∂ 2 ln [x, z |θ ] ∂c∂c = -Sψ ′ (c) ∂ 2 ln [x, z |θ ] ∂a∂b = S b ∂ 2 ln [x, z |θ ] ∂c∂d = S d ∂ 2 ln [x, z |θ ] ∂b∂b = - Sa b 2 ∂ 2 ln [x, z |θ ] ∂d∂d = -Sc d 2
G. Second derivative of the complete log-likelihood with discrete data

The complete log likelihood of the model, in the discrete case, reads : In addition here, ∆ z does not contain terms with θ, consequently the second order derivatives are easy to obtain and don't involve the latent variable; with Z standing for all the hidden variables i.e Z = (N, µ, p): and then the matrix B s is obtained by taking the covariance of this vector.Given N s additional advantage is taken from the conditional independence of p s and µ s (as shown on Figure 21 for the continuous case).

ln [x, z |θ ] = C -θ + (a -
Var     ln(µ s ) -µ s ln(ρ s ) -ρ s |N     =     -ψ ′ (a ′ ) 1 b ′ 0 0 1 b ′ -a ′ b ′2 0 0 0 0 ψ ′ (c ′ s ) -ψ ′ (c ′ s + d ′ s ) -ψ ′ (c ′ s + d ′ s ) 0 0 -ψ ′ (c ′ s + d ′ s ) ψ ′ (d ′ s ) -ψ ′ (c ′ s + d ′ s )    
and the expectation to obtain A s is performed via importance sampling.

A s =      E νs (ψ ′ (a ′ s )) -1 b ′ s 0 0 -1 b ′ s E N s+ |Y, θ (a ′ s ) b ′ s 2 0 0 0 0 E νs (ψ ′ (c ′ s ) -ψ ′ (c ′ s + d ′ s )) -E νs (ψ ′ (c ′ s + d ′ s )) 0 0 -E νs (ψ ′ (c ′ s + d ′ s )) E νs (ψ ′ (d ′ s ) -ψ ′ (c ′ s + d ′ s ))

Figure 1 :

 1 Figure 1: Realization of a marked Poisson process on a region of R 2 , the sample is conducted over a region D. Here the total catch is y = 7.7, the effective number of collected clumps is 8.

Figure 2 :

 2 Figure 2: DAG of the RLOL model

3 .Proposition 3 . 1 .

 331 Proofs of propositions are given in appendix E Assuming Y ∼ RLOL(θ ′ ) with θ ′ = (a ′ , b ′ , c ′ , d ′ ), S strata and I s records in stratum s as in 2.3 , then the complete conditional distributions of µ s and ρ s in one particular stratum s are given by

)

  where in stratum s, N s+ = Is i=1 N si denotes the total number of clumps caught, Y s+ = Is i=1 Y si is the entire quantity harvested and D s+ = Is i=1 D si is the whole catching effort.

Proposition 3 . 4 .

 34 Assuming Y ∼ RLOL(a, b, c, d) with S strata, and I s records in stratum s as in 2.3. Let us denote I e (θ) the empirical information matrix defined by

4. 1 .

 1 Real dataset -Gulf of St.Lawrence surveyA multi-species bottom-trawl survey of the southern Gulf of St.Lawrence (NW Atlantic) has been conducted each September since 1971. The purpose of this survey is to estimate the abundance and characterize the geographic distribution of marine biota. The survey follows a stratified random design, with 38 strata defined largely as homogeneous habitats using depth, temperature and sediments properties. The target fishing procedure at each fishing station is a 30-min straightline tow at a speed of 3.5 knots (i.e., 3.21km trawled distance). However the actual distance trawled can vary due to winds, currents and the avoidance of damaging rough bottoms; sampling effort is therefore variable among trawl tows, but this source of additional variability is easily accommodated in the models presented here ( the D s,k in eq 2.3). For our case study, we use data on the abundance of sea urchins and Sunflower starfishes collected during three survey years(1999)(2000)(2001), in a total of 540 bottom-trawl sets. The time period was chosen to minimize inter-annual changes in abundance while ensuring a sufficient sample size. The species were selected because inter-annual changes in their geographic distribution resulting from movements of individuals at the scale of survey sampling can be assumed to be approximately nil.The histograms of urchin and starfish catches in kg per survey tow clearly reflect zero-inflated distributions (Fig3 and 4). A large number of tows capture no urchin (nor starfish) and catches in non-zero tows tend to follow a skewed distribution. At the scale of the survey, sea urchins are distributed in patches of localized variable abundance, interspersed by numerous and relatively large areas where the species is absent (Fig5). Such patchy distributions of organisms are prevalent in ecological science. Data in two strata are always zero, thus rendering estimation impossible if we were to fit one model per stratum or to consider ρ s as fixed effects. Because

Figure 3 :Figure 4 :Figure 5 :

 345 Figure 3: Histogram of urchin biomass (kg/tow) from individual tows in the southern Gulf of St. Lawrence, bottom-trawl surveys: 1999-2000-2001

Figure 6 :Figure 7 :

 67 Figure 6: Comparisons between urchins dataset and averaged histogram (1000 simulations of datasets at θUrch )

Figure 8 :

 8 Figure 8: The ligthest ellipse corresponds to 90% confidence ellipsoid and the darkest one is 99.96% and contains 90% of the simulated values.

Figure 9 :

 9 Figure 9: The ligthest ellipse corresponds to 90% confidence ellipsoid and the darkest one is 94% and contains 90% of the simulated values.

Figure 10 :

 10 Figure 10: Predictions of the random effects µ s in each stratum correspond to the expected number of clumps collected during a measurement with standardcatching effort.

Figure 11 :

 11 Figure 11: Predictions of the inverse of ρ s in each stratum. These quantities give the expected biomass to be collected within a clump.

Figure 12 :

 12 Figure 12: pp-Plot with estimates of µ s versus a fitted gamma distribution.

Figure 13 :

 13 Figure 13: pp-Plot with estimates of ρ s versus a fitted gamma distribution. The extremal points correspond to strat with at least 75% of zeros

21

 21

For a given set

  of parameters θ = (a, b, c, d), we draw 1000 samples according to RLOL model given in eq 2.3 with a number S of strata and M measured points per each stratum. S has been chosen varying as k 2 with k = 3, 4, 5, 6, 8, 10, 12, 15 and M = 5, 10, 15, 20, 25, 30, 40.

Figure 15 :

 15 Figure 15: Urchins : Average relative bias in log scale depending on the number of strata and the number of measure points.

Figure 16 :

 16 Figure 16: Starfish : Average relative bias in log scale depending on the number of strata and the number os measure points.

Figure 17 :Figure 18 :

 1718 Figure 17: Urchin-like case. Effective proportion of 90% confidence intervals that cover the true value. Shading in particular cells reflects the degree of overlap: M-S combination that produces confidence intervals that are too liberal are in black whereas the lightest grey shade reflects confidence intervals that properly characterize parameter uncertainty

1 .

 1 Compound Poisson distributions can conveniently represent the presence of a large number of zeros and a skewed distribution of non-zero values. To deal the occurrence of zero-inflated data, very parsimonious models can be designed (with two parameters only) : a Poisson random sum of independent geometric random variables in the discrete case and with exponential random variables in the continuous one. They offer an alternative to the traditionnal delta gamma models and behave coherently when changing the scale of the catch effort, thanks to the Poisson process underpinning the model. 2. Compound Poisson distributions can be interpreted using a hierarchical framework. They describe the data collection involved in sampling individuals gathered in (latent) patches drawn from the homogeneous Poisson process with abundance tuned by the distributional parameter of the random components of the Poisson sum. The introduction of a random effect structure at the top of the hierarchy is straightforward and accommodates non homogeneity among strata that are themselves considered as homogeneous units. Such designs with random effects and data with extra zeros are commonly encountered in ecological analyzes, but gamma random effects are yet rarely advocated : variation between strata is typically modeled using a normal (or lognormal) distribution because its sufficient statistics match the commonsense interpretation of mean and variance. However, gamma random effects allow for partial conjugate properties with the compound Poisson model for zero-inflated data. Beyond this theoretical convenience, the parameters of the gamma distribution are well estimated in the Starfish like simulation examples and they can describe the entire range of variability between units for the real case study.

Figure 19 :

 19 Figure 19: Difference between log(x)ψ(x) and 1/2x

  [x, z |θ ] = [x |θ, z ] × [z |θ ] , but using Bayes rule, we may write the reverse logarithmic form : ln [x |θ ] = ln[x, z |θ ]ln [z |θ, x] (D.1) Let us remark that relation D.1 is valid whatever z represents. D.1. Recall about EM algorithm and control of the gradient Under regularity conditions for the joint distribution [x, z |θ ] and the conditional one [z |θ, x] , integrating relation D.1 with respect to the probability density [z |θ ′ , x] : ln

  x dz = ∂ ln[x| θ ] ∂θ j = 0, the second term in the right hand side of eq D.7 can be considered as the conditional variance of the gradient of the complete log-likelihood ln x, z θ . This expectation can be numerically computed with the same techniques to which recourse was made for the EM algorithm. E. Detailed proofs of propositions E.1. Proof of proposition 3.1

  µ|N, θ, y, D ∼ Γ(a + N + , b + D + ) ρ|N, θ, y ∼ Γ(c + N + , d + Y + ) For the expectation step we only need to compute E θ µ s | Y s , E θ ln(µ s ) | Y s and the same sufficient statistics concerning ρ. Since µ s |N s , θ, y follows a gamma distribution Γ(a + N s+ , b + D s+ ), the conditional expected value µ s given N s and θ = (a, b, c, d) is (a + N s+ )/(b + D s+ ).

2 )

 2 Computing the first term of the right hand side of equation E.2 is easy, since [x |z, θ ] = [x |z ] (consequently the complete log-likelihood ln([x, z |θ ]) can be separated as ln([x |z ]) + ln([z, |θ ])) and the gamma random effects [z|θ] belong to an exponential family. As a consequence, annex F shows that E ∂ 2 ln( x, Z θ )

Figure 20 :

 20 Figure 20: The random effects in each stratum are conditionally independent given the data and the set of parameters

F.

  Second derivative of the complete log-likelihood Let us first recall the complete log likelihood of the model : ln [x, z |θ ] = C -θ + (a -1) S s=1 ln µ s + Sa ln bb S s=1 µ s -S ln Γ(a) + (c -1) S s=1 ln ρ s + Sc ln dd S s=1 ρ s -S ln Γ(c)

  1) S s=1 ln µ s + Sa ln bb S s=1 µ s -S ln Γ(a) S ln Γ(c + d) Γ(c)Γ(d) + (c -1) S s=1 ln p s + (d -1) S s=1 ln(1p s )In the first derivative, the latent variables µ and p appear only through their arithmetic or geometric means (sufficient statistics for the gamma and beta pdf ). Using standard notation μ for the arithmetic mean 1 S S s=1 µ s , we have :∂ ln [x, z |θ ] ∂a = S ln(µ) + ln bψ(a) ∂ ln [x, z |θ ] ∂c = S ln(p) + ψ(c + d)ψ(c) ∂ ln [x, z |θ ] ∂b = S a b µ ∂ ln [x, z |θ ] ∂d = S ln(1p) + ψ(c + d)ψ(d)The gradient of the complete log-likelihood (so-called the "score") may be split into two parts : the first one ∆ θ does not depend on the latent variable z while the other one ∆ z gathers terms depending on z (and possibly of θ), i.e : + d)ψ(c) ψ(c + d)ψ(d)

∂ 2

 2 ln( x, z θ )∂θ i ∂θ j = S ′ (ĉ) + ψ ′ (ĉ + d) ψ ′ (ĉ + d) 0 0 ψ ′ (ĉ + d) -ψ ′ ( d) + ψ ′ (ĉ + d)As shown in Figure20for the continuous case , given Y s , Y s ′ and θ, the latent variables Z s and Z s ′ of two strata s and s ′ are conditionnaly independent, therefore :Var Z|x ∂ ln [x, Z |θ ] ∂θ = s ) -µ s ln(p s ) ln(1p s )To evaluate the variance of the score in stratum s, we will take advantage from successive conditioning due to the hierarchical structure depicted in Figure2still true for the discrete case. The variance conditional decomposition formula gives:′ (ĉ) + ψ ′ (ĉ + d) ψ ′ (ĉ + d) 0 0 ψ ′ (ĉ + d) -ψ ′ ( d) + ψ ′ (ĉ + d) Given N s ,µ s and ρ s are independent. Moreover the pdf [ρ s |N s , Y s , a, b, c, d] and [p s |N s , Y s , a, b, c, d] are gamma and beta so that analytic expressions are available for the expectation and variance of the gamma sufficient statistics, as detailed in equations C.1 to C.3. The key functions of N s+ are (a s ′, b s ′, c s ′, d s ′) = (a+ N s+ , b+ D s+ , c + N s+ , d + Y s+ -N s+ ) such that : s )ψ(c ′ s + d ′ s ) ψ(d ′ s )ψ(c ′ s + d ′ s )

616

  To sum it upI e (θ) = -∂ 2 ln [Y|θ] ∂θ i ∂θ j (G.1)At the maximum likelihood estimator θ, the following equality occurs :I e ( θ, Y) = S ′ (ĉ) + ψ ′ (ĉ + d) ψ ′ (ĉ + d) 0 0 ψ ′ (ĉ + d) -ψ ′ ( d) + ψ ′ (ĉ + d)

3 . 4 .

 34 s )ψ(c ′ s + d ′ s ) ψ(d ′ s )ψ(c ′ s + d ′ s ) s = â + N s+ , b ′ s = b + D s+ , c ′ = ĉ + N s+ and d ′ s = d + Y s+ -N s+ ( b ′ sis the only term that is not a function of N s+ , thus behaving like a constant with regards to the Var N s+ | θ,x operator)H. The discrete algorithmIf we adapt bluntly from the continuous version, the algoritm would write 1. Generate N (g) i = 0 wherever y i=0 fori = I -I + + 1, . . . , I.2. Generate a value of N + according toN + ∝ Γ(a ′ + N + )Γ(c ′ + N + )Γ(d ′ + Y + -N + )D N + + (b ′ + D + ) a ′ +N + I + j=1 Γ N + Y j D j (Y D) + Generate each N i for i = 1, .. . , I + , so that the vector N is distributed according to a multivariate hypegeometric Fisher distribution[START_REF] Mccullagh | Generalized Linear Models[END_REF] given by [N |N + ] = g(N; N + , Y , D/D + ) K N + with g(N; N + , Y , D) = I i=I * +1 Y j N j (D j /D + ) N j , K N + = y∈S g(y; N + , Y , D), and S = N ∈ Z I + + | I i=I * +1N i = N + . Associate to the vector the weightw (g) = K N (g)

Table 1 :

 1 1, we suggest to consider the median column as the refer-Comparison of the asymptotic 90% confidence interval with the one obtained by simulation for each parameter component for both species

	ence confidence interval (based on simulation/ EM re-estimation). The right col-
	umn gives bootstrap+ EM re-estimation. We notice that the Bootstrap approach

  N s , µ s and ρ s are independent. Moreover the pdf [ρ s |N s , Y s , a, b, c, d] and [µ s |N s , Y s , a, b, c, d] are gamma and analytic expressions are available for the expectation and variance of the gamma sufficient statistics, as detailed in equations C.1 to C.3. The key functions of N s+ are (a s ′, b s

APPENDICES

A. Compound Poisson process characteristic function

When X is real valued, we denote by f the Fourrier transform 1 of f (i.e the characteristic function of X) : f (ω) = E(e iωX ) From equation 2.1, the compound Poisson distribution g is such that :

1 For non negative integer valued random variables X the probability generating function

Pr(X = n)z n is the corresponding machinery for handling discrete distributions :

the same results can be found in this case by setting the change of variables z = e iω Importance Sampling relying this time on the multivariate hypergeometric distribution seems to stand naturally as the core of the algorithm to evaluate (3.23). But during our first trials, the above adaptation of the continuous version performed very badly, leading to a large variance of the importance weights, i.e. a degeneracy phenomenon that would put the all weight onto a very few contributing particles. In order to put more weight onto particles that have a good chance to efficiently attain the target distribution, a mixture was chosen as the importance distribution for a modified algorithm. The idea is similar in spirit to the auxiliary particle filtering of [START_REF] Pitt | Filtering via simulation : auxiliary particle filters[END_REF]. More precisely, the first step consists of determining an approximate mean of N s+ in stratum s, denoted N (ref )

s+ . One draws a L-sample of N s+ according to

The g distribution corresponds to the conditional distribution of N + given the sum of the data collected in stratum s but ignoring the individual records.

is given by the mean over a sample that is

and provides a good estimation of the location of N s+ . As previously we omit the index s to make the reading easier. Subsequently, the following algorithm relies on independent but non identically distributed simulations :

1. Generate N (g) i = 0 wherever y i = 0 for i = I -I + + 1, . . . , I. 3. Given µ (g) and ρ (g) , draw N

where K i denotes the normalizing constant.

4. Compute the weight of each particle g using