Catherine Larédo

Olivier David

Aurélie Garnier'

INFERENCE FOR PARTIALLY OBSERVED MULTITYPE BRANCHING PROCESSES AND ECOLOGICAL APPLICATIONS

Keywords: AMS 2000 subject classifications: Primary 62M05, 62M09, 62P12; secondary 92D40, 60J85 Parametric inference, Incomplete data, Multitype branching processes, ecology. 1

Multitype branching processes with immigration in one type are used to model the dynamics of stage-structured plant populations. Parametric inference is first carried out when count data of all types are observed. Statistical identifiability is proved together with derivation of consistent and asymptotically Gaussian estimators for all the parameters ruling the population dynamics model. However, for many ecological data, some stages (i.e. types) cannot be observed in practice. We study which mechanisms can still be estimated given the model and the data available in this context. Parametric inference is investigated in the case of Poisson distributions. We prove that identifiability holds for only a subset of the parameter set depending on the number of generations observed, together with consistent and asymptotic properties of estimators. Finally, simulations are performed to study the behaviour of the estimators when the model is no longer Poisson. Quite good results are obtained for a large class of models with distributions having mean and variance within the same order of magnitude, leading to some stability results with respect to the Poisson assumption.

1. Introduction. Understanding population dynamics requires models that admit the complexity of natural populations and the data ecologists can get from them. Thus analyzing ecological data raises questions ranging from modeling purposes to statistical inference. Among various methods, Leslie matrices or demographic matrix models are widely used for studying the dynamics of age or stage-structured populations (e.g. [START_REF] Caswell | Matrix Population Models[END_REF]. These models are deterministic with noise added to introduce some variability. In many cases however and especially for small populations, the demographic stochasticity has to be taken into account; these models are too simple [START_REF] Melbourne | Extinction risk depends stongly on factors contributiong to stochasticity[END_REF] and can no longer be used, even adding stochasticity into the dynamics with random effects and covariates [START_REF] Royle | N-mixture models for estimating population size from spatially replicated counts[END_REF][START_REF] Barry | The analysis of ring-recovery data using random effects[END_REF] or Bayesian approaches [START_REF] Raftery | Inference for a deterministic population dynamics model for bowhead whales[END_REF][START_REF] Gross | Bayesian estimation of a demographic matrix model from stage-frequency data[END_REF][START_REF] Clark | Population time series: Process variability, observation errors, missing values, lags, and hidden states[END_REF]. For these reasons, we use here stochastic models to study small populations dynamics.

The starting point of this work is a three-year field survey of feral populations (i.e. populations escaped from crops) of an annual crop species (oilseed rape) that was carrried out in the center of France (Selommes, Loir-et-Cher; [START_REF] Garnier | Dynamique et dispersion d'une espèce cultivée échappée des champs: le cas du colza[END_REF]. Unlike cultivated oilseed rape, very few facts are known about the dynamics of feral oilseed rape populations. In this study, the dynamics is modeled by a multitype branching process (five types including vegetative and reproductive plant stages along with seeds in the soil seedbank) with immigration in one type (seeds). Data consisted in populations counts in each type, except the seeds that could not be observed. Three main difficulties occur when studying this demographic dataset. (1) A large number of populations (K = 300) have been observed over a short period of time (n = 2, 3); (2) only count data have been collected; (3) some types could not be observed by ecologists (here seeds). These characteristics are clearly not specific to this survey and are frequently met in data coming from Population Genetics and Ecology (see e.g. [START_REF] De Valpine | Monte carlo state-space likelihoods by weighted posterior kernel density estimation[END_REF] and the references therein). These data could be studied as longitudinal data, but for concerns about the dynamics, better insights can be obtained by means of mechanistic models describing it.

Branching processes have largely been studied (see Athreya and Ney, 1972 for a general presentation; [START_REF] Mode | Multitype Branching Processes[END_REF] for Multitype Branching Processes and [START_REF] Haccou | Branching processes: Variation, growth, and extinction of populations[END_REF][START_REF] Kimmel | Branching Processes in Biology[END_REF][START_REF] Mode | Branching Processes in Epidemiology[END_REF] for applications in biology). Statistical inference has also been largely investigated [START_REF] Hall | Martingale Limit Theory and Its Applications[END_REF][START_REF] Guttorp | Statistical Inference for Branching Processes[END_REF] for general branching processes; [START_REF] Wei | Estimation of the means in the branching process with immigration[END_REF][START_REF] Winnicki | Estimation of the variances in the branching process with immigration[END_REF] for branching processes with immigration; [START_REF] Bhat | Maximum likelihood estimation for branching processes with immigration[END_REF][START_REF] Maaouia | Identification of multitype branching processes[END_REF][START_REF] Gonzalez | Non-parametric bayesian estimation for multitype branching processes through simulation-based methods[END_REF] for multitype branching processes). However, the precise multitype branching process with immigration used here is a combination of the previous ones and moreover statistical inference for multitype branching processes is usually based on the following different observations [START_REF] Maaouia | Identification of multitype branching processes[END_REF][START_REF] Gonzalez | Non-parametric bayesian estimation for multitype branching processes through simulation-based methods[END_REF]: the number of descendants of type j coming from all the type i individuals. We just observed the successive counts of "individuals" of each type. This is more realistic assumption since this situation frequently occurs with datasets from field studies, and inference is studied here in this framework.

We are interested in the estimation of the parameters involved in the population dynamics from the incomplete observations of count data collected simultaneously in several populations. This is an "Incomplete Data model", or "State Space model" (as defined for instance in [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]. It is also an inverse problem and a central theme in Ecology arising from its study is that parameters might not be identifiable knowing only the population dynamics [START_REF] Wood | Inverse problems and structured populations dynamics[END_REF]. In practice, the inference based on such data is performed using various E.M. algorithms eventually coupled with Monte Carlo methods [START_REF] Dempster | Maximum likelihood fron incomplete data via the em algorithm[END_REF][START_REF] Kuhn | Coupling a stochastic approximation version of EM with an MCMC procedure[END_REF][START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF][START_REF] Sung | Monte carlo likelihood inference for missing data models[END_REF][START_REF] Olsson | Sequential monte carlo smoothing to parameter estimation in non linear state space models[END_REF], and Bayesian or Hierarchical Bayesian methods [START_REF] Clark | Population time series: Process variability, observation errors, missing values, lags, and hidden states[END_REF][START_REF] Buckland | State-space models for the dynamics of wild animal populations[END_REF][START_REF] Thomas | A unified framework for modelling wildlife population dynamics[END_REF]. All these methods circumvent but cannot address the identifiability problem. However, identifiability is a prerequisite of statistical inference, and understanding the dynamics mechanisms strongly relies on how parameters are linked in the identifiability question. We propose here an integrated framework in order to analyze as accurately as possible the whole data set of the field survey. Introduction of covariates and a priori knowledge, errors coming from non exhaustive population samplings within some populations, use of various algorithms rely on this work and are studied in two companion papers [START_REF] David | Estimation of plant demographic parameters from stage-censuses[END_REF][START_REF] Garnier | Are seed immigration and local seed production efficient to favour the persistence of feral oilseed rape populations submitted to weed control? in preparation[END_REF].

The paper is organized as follows. Section 2 contains the description of the population dynamics and preliminary results (Proposition 2.1). The statistical inference for complete observations is studied in Section 3. We first prove identifiability for all the parameters and derive consistent and asymptotically Gaussian estimators (Proposition 3.1). Since seeds are not observed in practice, the problem of unobserved types is addressed in Section 4. This is a non linear non Gaussian state space model. The associated three-dimensional stochastic process is no longer Markovian. The model with Poisson distributions provides a useful example with explicit computations. We obtain a closed form of the dependence of present observations on the whole past for the three-dimensional process (Theorem 4.1). A question concerns the statistical model identifiability : it is studied according to the number of observed generations. We characterize the parameter subset where identifiability holds (Theorem 4.2) and study the parametric inference (Proposition 4.2). Section 5 presents simulation results to study how the estimation performs with respect to deviations to the Poisson model. Detailed proofs are given in the Appendix.

Model and preliminary results.

2.1. Dynamics of annual plants. We consider annual plants with the following life cycle. Seeds are released at the end of summer; they can either enter in a seed bank if buried or germinate in autumn.The emerged rosettes vernalize during winter, then bolt in spring and finally produce mature plants that shed seeds in summer and then die. Five developmental stages are considered: rosettes before winter R, rosettes after winter (vernalized rosettes) V , mature plants F , seeds located in the soil seed bank ("old seeds") S, and seeds located on the soil surface ("new seeds") T . New seeds and old seeds are separated because they have different demographic parameters. Within each cycle, new seeds can enter these populations at the end of summer (immigration). There exist two sources of seed immigration, seeds from adjacent mature crops and seeds from spillage during seed transport [START_REF] Crawley | Seed limitation and the dynamics of feral oilseed rape on the m25 motorway[END_REF]Claessen et al., 2005a;[START_REF] Garnier | Stochastic modelling of feral plant populations with seed immigration and road verge management[END_REF][START_REF] Pivard | Where do the feral oilseed rape populations come from? A large-scale study of their possible origin in a farmland area[END_REF]. This model is quite general and for dynamical purposes only, it could be simplified considering just seeds and mature plants. However, our concern is different since we aim at estimating as many parameters as possible given both the model and the data available. Keeping a five-type model allows using all the data collected in the field survey and thus leads to the best description of the plant population dynamics for inference.

2.2. Notations and assumptions. Consider first one population. From now on, the term year corresponds to one life cycle. It starts with the birth of the new seeds and ends just before the birth of the new seeds of the next generation. All the variables are integer random variables indexed by i ∈ N. For year i, denote by S i the number of "old seeds", T i the number of "new seeds", R i the number of rosettes before winter, V i the number of rosettes after winter and F i the number of mature plants. Six parameters describe these transitions:

(2.1) (c, d, a, b, a ′ , b ′ ) ∈ (0, 1) 6 with 0 < a + b ≤ 1, 0 < a ′ + b ′ ≤ 1. (2.2) P (seed in S i → S i+1 ) = a ; P (seed in S i → R i ) = b (2.3) P (seed in T i → S i+1 ) = a ′ ; P (seed in T i → R i ) = b ′ (2.4) P (rosette in R i → V i ) = c ; P (rosette in V i → F i ) = d.
Mature plants in F i produce "new seeds" T ′ i+1 according to the offspring distribution G(.).

A number I i of "new seeds" immigrate into the population at the beginning of year i; it is assumed to follow the distribution µ(.). Seeds in S i come from two sources, the "old seeds" S i-1 and the "new seeds" T i-1 . Denote by S ′ i and S ′′ i these two quantities. Similarly, rosettes before winter R i come from "old seeds" in S i and "new seeds" in T i . Denote by R ′ i and R ′′ i the rosettes coming from S i and T i . These variables satisfy :

(2.5)

S i = S ′ i + S ′′ i , R i = R ′ i + R ′′ i ,
(2.6)

T i = T ′ i + I i .
Note that the probabilities of dying for stages S, T, R, V are respectively

(1 -a -b), (1 -a ′ -b ′ ), (1 -c), (1 -d)
and that the probability of no offspring for a mature plant is G(0).

Let us now detail the framework and assumptions used in the sequel. The field survey consisted of a large number of feral oilseed rape populations (around 300 in [START_REF] Garnier | Dynamique et dispersion d'une espèce cultivée échappée des champs: le cas du colza[END_REF] observed over a short period of time (n = 2 or 3). These populations were isolated, so we assume here independence for these populations and the plant density was low in the surveyed populations so that density-dependence in plant survival and reproduction could be neglected [START_REF] Pivard | Where do the feral oilseed rape populations come from? A large-scale study of their possible origin in a farmland area[END_REF][START_REF] Garnier | Dynamique et dispersion d'une espèce cultivée échappée des champs: le cas du colza[END_REF]. Moreover, we assume Assumption 1. Offspring distribution G(.) All mature plants reproduce independently according to the same offspring distribution G(.) with expectation m and variance δ 2 . Assumption 2. There is no competition in survival and germination between seeds in the seed bank or "old seeds" S i and seeds on the soil or "new seeds" T i .

Assumption 3. Immigration distribution µ(.) (i) Immigration I i is independent of seed bank seeds S i , offspring seeds T ′ i and of previous years. (ii) The random variables (I i , i ∈ N ) are independent and identically distributed according to µ(.) with expectation u and variance ρ 2 .

The above assumptions are twofold: "independence" and "identically distributed". They do not have the same status : while releasing the "independence " assumption is quite difficult, the "identically distributed " assumption is done here for sake of clarity. The independence assumption in (A1) is justified because of the low plant density. There is no biological background for considering competition between the evolution of "old seeds" S i and "new seeds" T ′ i leading to (A2). In the field survey of feral oilseed rape populations, seed immigration mainly occured from spillage during seed transport and from adjacent mature crops, leading to the independence assumptions in (A3) [START_REF] Pivard | Where do the feral oilseed rape populations come from? A large-scale study of their possible origin in a farmland area[END_REF][START_REF] Garnier | Dynamique et dispersion d'une espèce cultivée échappée des champs: le cas du colza[END_REF]. Adding covariates or a priori information can easily be introduced within this framework, which amounts to remove the "identically distributed" assumption. This indeed has been done for the statistical analysis on the whole data set: using that many populations had been observed, covariates, a priori knowledge, and a dependence with respect to i, k of the parameters defined in (2.2)-(2.4) have been added [START_REF] Garnier | Are seed immigration and local seed production efficient to favour the persistence of feral oilseed rape populations submitted to weed control? in preparation[END_REF]. The framework detailed here is therefore quite general.

2.3. Preliminary results. From now on, time i is associated with a complete cycle of a plant. Let us still consider one population. The model is a discrete time stochastic process (X i ) with state space N 5 . Set (2.7)

X i = (S i , T i , R i , V i, F i ) and F i = σ(X 0 , X 1 , . . . , , X i ).
For x = (s, t, r, v, f ), x ′ = (s ′ , t ′ , r ′ , v ′ , f ′ ) ∈ N 5 , denote π 0 (x) the initial distribution and p(x, x ′ ) the conditional distribution of X 1 given X 0 :

(2.8)

π 0 (x) = P(X 0 = x) ; p(x, x ′ ) = P(X 1 = x ′ /X 0 = x).
For the Binomial distribution B(N ; p), we write P(Y = k) = B(N ; p)(k); Multinomial distributions on N l , M(N ; p 1 , . . . , p l ) are simplified omitting the last component, which leads for l = 3 and 0 < p

1 + p 2 < 1, (2.9) M(N ; p 1 , p 2 )(n 1 , n 2 ) = M(N ; p 1 , p 2 , 1 -p 1 -p 2 )(n 1 , n 2 , N -n 1 -n 2 )
Let ⋆ denote the convolution product of two distributions. Using these notations and (2.1), define the distributions (2.10) ν(s, t) = P(S 0 = s, T 0 = t),

(2.11)

p 1 (s ′ /s, t, r) = (M(s; a, b) ⋆ M(t; a ′ , b ′ ))(s ′ , r) (B(s; b) ⋆ B(t; b ′ ))(r) ,
(2.12)

p 2 (t ′ /f ) = (G ⋆f ⋆ µ)(t ′ ); p 3 (r/s, t) = (B(s; b) ⋆ B(t; b ′ ))(r), (2.13) p 4 (v/r) = B(r; c)(v) , p 5 (f /v) = B(v; d)(f ).
Using now definitions (2.7), (2.8) and notations (2.9)-(2.13), the following holds.

Proposition 2.1. Under (A1), (A2), (A3), (X i ) i≥0 is a time homogeneous Markov chain on N 5 with initial distribution π 0 (x) and transition probabilities p(x, x ′ ) satisfying

π 0 (x) = ν(s, t) p 3 (r/s, t) p 4 (v/r) p 5 (f /v) (2.14) p(x, x ′ ) = p 1 (s ′ /s, t, r) p 2 (t ′ /f ) p 3 (r ′ /s ′ , t ′ ) p 4 (v ′ /r ′ ) p 5 (f ′ /v ′ ). (2.15)
The process (X i ) i≥0 is also a multitype branching process with immigration in one type.

The last statement of Proposition 2.1 is immediate since, considering each stage of the plant as a type, each plant reproduces independently from the others according to the same offspring distribution with values in N 5 . However, the two types R i and V i have no offspring in the next generation, leading to a non positively regular multitype process as defined in [START_REF] Mode | Multitype Branching Processes[END_REF] or [START_REF] Athreya | Branching Processes[END_REF]. The process (S i , T i , F i ) is positively regular process but, for the reasons stated in 2.1, we prefer keeping here the five-dimensional process (X i ).

The proof of Proposition 2.1 is given in Appendix A.1.

3. Likelihood and inference for complete observations. We first study the case when all types are observed.

3.1. Notations and statistical framework. We assume in the sequel that the initial distribution of (S 0 , R 0 ), the offspring distribution G(.) and the immigration µ(.) belong to the parametric families : -distribution of (S 0 , T 0 ) : ν(θ 1 ; s, t), θ 1 ∈ Θ 1 ; -offspring G(.): G(θ 2 ; .) , θ 2 ∈ Θ 2 with mean m and variance δ 2 ; -immigration µ(.) : µ(θ 3 ; .), θ 3 ∈ Θ 3 ) with mean u and variance ρ 2 .

Let us denote by θ = (c, d, a, b, a ′ , b ′ , θ 1 , θ 2 , θ 3 ) (resp. θ 0 ) an arbitrary value (resp. the true value) of the parameter and by Θ the parameter set. We assume :

Assumption 4. Θ compact set of R l and θ 0 ∈ • Θ.
For the kth population, X k = (X k i , i ∈ N) is the Markov chain describing the population dynamics and

x k i = (s k i , t k i , r k i , v k i , f k i )
are the observations at time i. In order to simplify the expressions for the likelihood, we consider here that S i , T i are observed up to time n + 1. This has no consequence since seeds are not observed in practice. Hence we denote,

X k 0:n = X k 0 , . . . , X k n , S k n+1 , T k n+1 , (3.1) O k 0:n = (x k 0 , x k 1 , . . . , x k n , s k n+1 , t k n+1 ). (3.2)
The processes X k 0:n are repetitions of X 0:n = (X 0 , . . . , X n , S n+1 , T n+1 ). Joining all the populations, define :

(3.3) X 0:n (K) = (X 1 0:n , . . . , X K 0:n ), and O 0:n (K) = (O 1 0:n , . . . , O K 0:n ). Let π 0 (θ; x) (resp. p(θ; x, x ′
)) be the initial distribution (resp. the transition probabilities) associated with parameter θ and P θ the probability distribution of (X i ) on the canonical space and E θ the expectation w.r.t. P θ .

Likelihood.

Computing the likelihood of a Markov process with transition probabilities p(θ; x, x ′ ) is classical : for population k, it has for expression, using (2.14),

(3.4) L(θ; O k 0:n ) = π 0 (θ; x k 0 ) n i=1 p(θ; x k i-1 , x k i ) p 1 (s k n+1 /x k n )p 2 (t k n+1 /x k n ).
Joining the observations from the K independent populations, we obtain using Proposition 2.1, (2.10) -(2.13)

(3.5) l(θ; O 0:n (K)) = l K (θ) = K k=1 log L(θ; O k 0:n ).
Reordering the terms of (3.4) and (3.5) according to the parameters yields,

(3.6) l(θ; O 0:n (K)) = 5 j=0 l j K (θ).
The first term deals with the initial distribution of (S 0 , T 0 )

(3.7) l 0 K (θ) = K k=1 log ν(θ 1 ; s k 0 , t k 0 ) = l 0 K (θ 1 ).
Gather in the second term the transition from seeds S, T to rosettes R.

(3.8)

l 1 K (θ) = K k=1 n i=0 log(B(s k i ; b) ⋆ B(t k i ; b ′ ))(r k i ) = l 1 K (b, b ′ ).
The next two terms contain the transitions from R to V and V to F , they write

(3.9) l 2 K (θ) = K k=1 n i=0 v k i log c + (r k i -v k i ) log(1-c)+ C 2 (O 0:n (K)) = l 2 K (c), (3.10) l 3 K (θ) = K k=1 n i=0 f k i log d+(v k i -f k i ) log(1-d)+C 3 (O 0:n (K)) = l 3 K (d),
where the two terms C 2 (O 0:n (K)) and C 3 (O 0:n (K)) only depend on the observations. Set in the next term the transition from F to T ,

(3.11) l 4 K (θ) = K k=1 n i=0 log G(θ 2 ; .) ⋆f k i ⋆ µ(θ 3 ; .) (t k i+1 ) = l 4 K (θ 2 , θ 3 )
The last term concerns the seeds S :

(3.12)

l 5 K (θ) = K k=1 n i=0 log( (M(s k i ; a, b) ⋆ M(t k i ; a ′ , b ′ ))(s k i+1 , r k i ) (B(s k i ; b) ⋆ B(t k i ; b ′ ) (r k i ) ) Joining the two terms containing a, a ′ , b, b ′ yields (3.13) l 6 K (θ) = l 1 K (θ) + l 5 K (θ) = l 6 K (a, b, a ′ , b ′ ).
The terms l 0 k (θ), l 2 K (θ), l 3 K (θ), l 4 K (θ) and l 6 K (θ) depend on disjoint sets for the parameters. Hence, maximizing the loglikelihood can be performed maximizing separately these five terms.

Remark 3.1. Usually, statistical inference for Stochastic Processes is investigated in the asymptotic framework small K (mostly K = 1) and large n (leading to asymptotics results n → +∞). Here, we have that n is small and K large (e.g. magnitude 300). This situation often occurs in Ecology [START_REF] De Valpine | Monte carlo state-space likelihoods by weighted posterior kernel density estimation[END_REF].

3.3.

Study of maximum likelihood and other estimators. This is a K sample of i.i.d. random variables, each variable being a part of a branching process path. We have to use simultaneously the repetitions and the Markov structure to estimate the parameters. So, deriving the properties of the statistical model is not standard. The various terms l i K (θ) of the loglikelihood are associated with different parametric inference problems.

The first term l 0 K (θ) deals with the estimation of θ 1 based on a sample of K i.i.d. random variables with distribution ν(θ 1 ; .) on N 2 , which is standard. The terms l 2 K (θ) = l 2 K (c) and l 3 K (θ) = l 3 K (d) are related to the estimation of parameters c, d. Denote by ĉK , dK the maximum likelihood estimators (MLE) obtained maximizing l 2 K (c) and l 3 K (d). They depend on the successive observations

(r k i , v k i ) (resp. (v k i , f k i ))
for {i = 0, . . . , n; k = 1, . . . , K} and are explicit (see Appendix A.2).

Parameters (a, b, a ′ , b ′ ) are only present in l 6 K (θ) defined in (3.13). Maximum likelihood estimators for (a, b, a ′ , b ′ ) can be defined maximizing l 6 K (θ). To prove identifiability and consistency, we have rather consider here conditional least squares (CLS) estimators. Conditionally on (S i , T i ), the marginal distribution of S i+1 (resp. R i ) is the sum of the two independent distributions, B(S i , a) and B(T i , a ′ ) (resp. B(S i , b) and B(T i , b ′ )). Therefore, we can define the CLS estimators (â K , â′ K ) and ( bK , b′ K ) minimizing the Conditional Least Squares:

J 1 K (a, a ′ ) = K k=1 n i=0 (s k i+1 -as k i -a ′ t k i ) 2 ; (3.14) J 2 K (b, b ′ ) = K k=1 n i=0 (r k i -bs k i -b ′ t k i ) 2 . (3.15)
The remaining term is l 4 K (θ 2 , θ 3 ) = l 4 K (m, u) since we are concerned by the estimation of m and u. This is the only part of the likelihood associated with the branching mechanism. This likelihood containing the convolution product G ⋆f ⋆ µ is untractable, and methods based on conditional least squares or weighted conditional least squares are used for branching processes [START_REF] Hall | Martingale Limit Theory and Its Applications[END_REF][START_REF] Guttorp | Statistical Inference for Branching Processes[END_REF][START_REF] Wei | Estimation of the means in the branching process with immigration[END_REF]) leading to moment estimations of G and µ. Noting that E(T i+1 /F i ) = mF i + u, we can just consider for the estimation of m and u, the CLS process,

(3.16) J 4 K (m, u) = K k=1 n i=0 (t k i+1 -mf k i -u) 2 .
Let ( mK , ûK ) be the CLS estimators minimizing (3.16). All the above estimators are explicitly defined in Appendix A.2 and we can state :

Proposition 3.1. Assume (A1), ( A2), ( A3) and (A4). Then, under P θ 0 , all the parameters c, d, a, b, a ′ , b ′ , m, u are identifiable and, as K → ∞, (i) (ĉ K , dK , âK , bK , â′ K , b′ K , mK , ûK ) are consistent and asymptotically Gaussian at rate √ K; (ii) ĉK , dK , (â K , bK , â′ K , b′ K ), ( mK , ûK ) are asymptotically independent, with explicit covariance matrix given in Appendix A.2.

Let us stress that, before studying in detail this inference problem, it was difficult to assert the classical properties stated in Proposition 3.1. Adding immigration could lead to non identifiability or estimating problems for m and u. Moreover, maximum likelihood estimators, for multitype branching processes, are based on the observations of G i,j (k), i.e. offspring of type j from type i parents (see [START_REF] Guttorp | Statistical Inference for Branching Processes[END_REF][START_REF] Gonzalez | Non-parametric bayesian estimation for multitype branching processes through simulation-based methods[END_REF][START_REF] Maaouia | Identification of multitype branching processes[END_REF]. We did not require this information for the inference and just used the counts of individuals of each type in successive generations. Hence, getting identifiability and consistency for the parameters is the only difficulty here. Other properties are classical but requires using the exact structure of the data provided the regularity of the statistical model.

Remark 3.2. Estimating additional moments of G(.) and µ(.) can be performed similarly using other functionals than Conditional Least Squares (see [START_REF] Winnicki | Estimation of the variances in the branching process with immigration[END_REF] for the variance estimation of G(.) and µ(.) ).

The proof is given in Appendix A.2.

4. Incomplete model study in the Poisson case. The set-up is now different: only R i , V i , F i are observed while S i and T i are unobserved. Clearly, algorithms simulating the missing data given the model and the parameters at each step can be used to get estimation. This approach is complementary to our concern that aims at understanding which mechanisms can be estimated. For this, we have to study the process (R i , V i , F i ) for i = 0, . . . , n. It is a discrete time stochastic process, which is no longer Markov : the distribution of (R i+1 , V i+1 , F i+1 ) given the past now depends on the whole past and not only on (R i , V i , F i ). This appears explicitely later on. This is similar to problems encountered when studying Hidden Markov Models [START_REF] Genon-Catalot | Conditional likelihood estimators for hidden markov models and stochastic volatility models[END_REF][START_REF] Genon-Catalot | Leroux's method for general hidden markov models[END_REF][START_REF] Cappé | Inference in Hidden Markov Models[END_REF]. For a first approach, we restrict our attention to a very informative example, the case of Poisson distributions, which leads to explicit computations.

Probabilistic properties in the Poisson case.

Let us specify all the distributions appearing in the populations dynamics.

Assumption 5. The offspring distribution G(.) is Poisson P(m), the immigration distribution µ(.) is Poisson P(u). Assumption 6. The variables S 0 and T 0 are independent and distributed according to Poisson distributions: S 0 ∼ P(σ) and T 0 ∼ P(τ ).

For Poisson distributions, we denote P(λ)(k) = P(X = k). Recall a property of Multinomial and Poisson distributions.

Lemma 4.1. Assume that N is a random variable distributed according to a Poisson distribution P(λ) and that X = (X 1 , X 2 , . . . , X l ) is a l-dimensional random variable such that the conditional distribution of X given N is a Multinomial distribution M(N ; a 1 , . . . , a l ) with l i=1 a i = 1. Then, the random variables {X i , i = 1, . . . , l} are independent and verify X i ∼ P(a i λ).

Y i = (R i , V i , F i ) and G i = σ(Y j , j = 0, . . . , i).
Clearly, G i is the information available up to time i. To state the main result of this section, define the three sequences of G i-1 measurable random variables:

(4.2) Γ 0 = σ ; Γ ′ 0 = τ ; (4.3) for i ≥ 1, Γ i = aΓ i-1 + a ′ Γ ′ i-1 ; Γ ′ i = mF i-1 + u. (4.4) for i ≥ 0, Λ i = bΓ i + b ′ Γ ′ i .
Then, the following holds :

Theorem 4.1. Under Assumptions (A1)-( A6), the initial distribution π0 (y) of (Y i ), and the conditional distribution L(Y i+1 /G i ) satisfy, using (4.1) and definitions (2.13),(4.4), for y = (r, v, f

) ∈ N 3 P(Y 0 = (r 0 , v 0 , f 0 )) = π0 (y 0 ) = P(Λ 0 )(r 0 )p 4 (v 0 /r 0 )p 5 (f 0 /v 0 ), (4.5) P(Y i+1 = (r, v, f )/G i ) = P(Λ i+1 )(r)p 4 (v/r)p 5 (f /v). (4.6)
The explicit dependence of R i on the whole past (Y 0 , . . . , Y i-1 ) appears more simply with the following expression for Λ i ,

Λ 0 = bσ + b ′ τ = c 0 (θ) ; (4.7) Λ 1 = b ′ mF 0 + c 1 (θ) with c 1 (θ) = abσ + a ′ bτ + b ′ u, (4.8) Λ i = b ′ mF i-1 + a ′ bm(F i-2 + aF i-3 + • • • + a i-2 F 0 ) + c i (θ), with (4.9) c i (θ) = a i bσ + a i-1 a ′ bτ + a ′ bu 1 -a i-1 1 -a + b ′ u for i ≥ 2. (4.10)
The result of Theorem 4.1 is a consequence of the proposition stated below.

Proposition 4.1. Under Assumptions (A1)-(A6), the random variables S i+1 , T i+1 are conditionally independent given G i , and their conditional distributions satisfy, using the random variables Γ i and

Γ ′ i defined in (4.2), (4.3), (4.11) L(S i+1 /G i ) ∼ P(Γ i+1 ) and L(T i+1 /G i ) ∼ P(Γ ′ i+1 ).
Let us prove Theorem 4.1 assuming Proposition 4.1. We just have to check the expression of the conditional distribution of R i+1 . By (2.5), we have

R i+1 = R ′ i+1 + R ′′ i+1 .
Using Proposition 2.1 and (2.12), the distribution of R i+1 conditionnally on (S i+1 , T i+1 ) is equal to p 3 (r/S i+1 , T i+1 ). Applying Proposition 4.1, S i+1 and T i+1 are conditionally independent given G i and distributed according to two independent Poisson distributions. Hence, another application of Lemma 4.1 yields that the conditional distribution of

R i+1 given G i is P(bΓ i+1 + b ′ Γ ′ i+1 ) = P(Λ i+1
), which is (4.4). The proof of Proposition 4.1 is given in Appendix A.3. 4.2. Likelihood of the incomplete observations. The inference is now based on the observations of (Y i ) recorded up to time n for K independent populations. Denote by

Y k i = (R k i , V k i , F k i )
the process describing its dynamics in population k and

y k i = (r k i , v k i , f k i ) the observations at time i. We set, (4.12) Y k 0:n = (Y k 0 , . . . , Y k n ), and Y 0:n (K) = (Y 1 0:n , . . . , Y K 0:n ).
Observations up to time n are denoted (4.13) Õk 0:n = (y k 0 , y k 1 , . . . , y k n ) and Õ0:n (K) = ( Õ1 0:n , . . . , ÕK 0:n ).

Let us first compute the likelihood for one population, population k. Successive conditionnings yield

(4.14) L(θ; Õk 0:n ) = P θ (Y k 0 = y k 0 ) n i=1 P θ (Y k i = y k i /Y k 0:i-1 = y k 0:i-1 ).
Contrary to the previous section, each term of this product depends on i and on the observations up to time i -1. Theorem 4.1 gives the expression of these conditional probabilities. Since the random variables Λ i now depend on θ and on the past, we define Λ i (θ) = Λ i (θ; Y 0:i-1 ), and for population k,

(4.15) Λ k i (θ) = Λ i (θ; Y k 0:i-1 ) , λ k i (θ) = Λ i (θ; y k 0:i-1 ). (4.16) P θ (Y k i = y k i /y k 0:i-1 ) = P(λ k i (θ))(r k i ) p 3 (c; v k i /r k i ) p 4 (d; f k i /v k i ).
Joining the observations in the K populations, the likelihood writes, using notations (4.12), (4.13), (4.16), (4.17)

L(θ; Õ0:n (K)) = K k=1 L(θ; Õk 0:n ).
The log-likelihood splits into three terms, (4.18) l(θ, Õ0:n (K)) = log L(θ; Õ0:n (K)) =

3 i=1 l i K (θ, Õ0:n (K)),
where, using Assumptions (A5), (A6) and Theorem 4.1,

(4.19) l 1 K (θ) = l 1 K (θ, Õ0:n (K)) = K k=1 n i=0 log P(λ k i )(r k i ), (4.20) l 2 K (θ, Õ0:n (K)) = l 2 K (c, O 0:n ) ; l 3 K (θ, Õ0:n (K)) = l 3 K (d, O 0:n ).
Hence, estimating parameters c, d is exactly the same as in the previous section; their inference is omitted in the sequel.

Parametric inference.

It remains to study the estimation of the parameter θ = (σ, τ, a, b, a ′ , b ′ , u, m). As before, let θ 0 be true value of the parameter.

We first have to investigate which parameters are identifiable when only these incomplete observations are available. By identifiability, we mean here identifiability of a statistical model M = (P θ , θ ∈ Θ) :

∀θ, θ ′ ∈ Θ, {P θ = P θ ′ } ⇒ {θ = θ ′ }.
Recall that n denotes the time index of a plant lifecycle and that we consider populations recorded up to time n (n = 0 corresponding here to the observation of one complete cycle). Using now the definitions of the terms c i (θ) given in (4.7), the following holds : Theorem 4.2. Assume (A1), (A2), (A3), (A4), (A5), (A6). Then,

1. if n = 0, only c 0 (θ) = bσ + b ′ τ is identifiable; 2. if n = 1, (b ′ m, c 0 (θ), c 1 (θ)) is identifiable; 3. if n = 2, ( a ′ b b ′ , b ′ m, c 0 (θ), c 1 (θ), c 2 (θ)), is identifiable; 4. if n ≥ 3, and a = a ′ b b ′ , then φ = (a, a ′ b b ′ , b ′ m, b ′ u, bσ, b ′ τ ) is identifiable; if n ≥ 3 and a = a ′ b b ′ , then (a = a ′ b b ′ , b ′ m, b ′ u, bσ + b ′ τ ).
Note that identifiability of additional parameters cannot be gained increasing n beyond 3. Larger values of n result in improving the asymptotic variance for the estimation of φ.

Remark 4.1. Stating the above theorem for the first values of n is unusual. However, in the field survey of feral oilseed rape populations, observations unfortunately had been collected up to n = 2, leading to the unability of estimating a, the annual survival rate in the seed bank, which is a parameter of much concern in Ecology. Using (4.19), define the maximum likelihood estimator φK as a solution of (4.21)

l 1 K ( φK ) = sup{ l 1 K (φ) ; φ ∈ Φ}.
Under Assumption (A7), the function φ

→ n 0 (Λ i (φ; Y 0:i-1 )-log Λ i (φ, Y 0:i-1
)) is a.s. P φ twice differentiable on Φ, and we can define, for φ = (φ 1 , . . . , φ 6 ), the 6 × 6 matrix,

I(φ) p,q = n i=0 E θ 0 ( 1 Λ i (φ) ∂Λ i (φ) ∂φ p ∂Λ i (φ) ∂φ q
) with 1 ≤ p, q ≤ 6.

Proposition 4.2. Assume (A1)-( A3), ( A5), ( A6) and (A7). Then φK is strongly consistent. If moreover the matrix I(φ 0 ) is invertible, then

√ K( φK -φ 0 ) D → N (0, I(φ 0 ) -1 ) under P φ 0 as K → +∞.
Remark 4.2. The matrix I(φ 0 ) that appears in the asymptotic variance of φK can be estimated, using the explicit expressions for the derivatives of Λ i (φ), by the empirical estimate for 1 ≤ p, q ≤ 6,

Îp,q = 1 K K k=1 n i=0 1 λ k i ( φK ) ∂λ k i ( φK ) ∂φ p ∂λ k i ( φK ) ∂φ q .
The proofs of Theorem 4.2 and Proposition 4.2 are given in Appendix A.4.A.5.

5. Simulation study. In the case of incomplete observations, the estimators strongly rely on the assumption that both the offspring and the immigration distributions are Poisson distributions. We investigate in this section how the estimation behaves when these distributions are no longer Poisson. 5.1. Methods. We considered in these simulations that the offspring distribution G(.) or the immigration distribution µ(.) are more dispersed than Poisson distributions. Two series of simulations were performed: the first one concerns deviations to a Poisson distribution of the offspring G(.), and the other one for the immigration µ(.). For this, we used Negative Binomial distributions for G(.) (resp. µ(.)) with mean m (resp. u) and increasing values of the variance variance δ 2 (resp. ρ 2 ); the (variance /mean) ratio ranged from 2 to 1000. For each given set of parameters, we performed M = 100 repetitions including K = 300 populations. These populations dynamics were run during n = 4 years to get the observations for rosettes R 0 to R 4 and mature plants (F 0 to F 3 ), using biologically plausible values for demographic parameters. Indeed, for a, b, a ′ , b ′ , we used the parameters given in Claessen et al., 2005a,b : (5.1) a = 0.15 , a ′ = 0.006 , b = 0.5 , b ′ = 0.5.

We also had to fix some values for c, d, m = E(G), u = E(µ). We used values estimated in [START_REF] Garnier | Are seed immigration and local seed production efficient to favour the persistence of feral oilseed rape populations submitted to weed control? in preparation[END_REF] (5.2) c = 0.21 , d = 0.01 , m = 13 , u = 80.

The value m = 13 corresponds to the mean fecundity of plants mown twice. [START_REF] Colbach | Genesys : a model of the influence of cropping system on gene escape from herbicide tolerant rapeseed crops to rape volunteers-i. temporal evolution of a population of rapeseed volunteers in a field[END_REF]. The value u = 80 corresponds to average immigration when there is no cultivated field in the neighbourhood [START_REF] Garnier | Are seed immigration and local seed production efficient to favour the persistence of feral oilseed rape populations submitted to weed control? in preparation[END_REF]. Finally, for the initial distributions of S 0 , T 0 , we assumed (5.3) S 0 ∼ P(σ), T 0 ∼ P(τ ) with σ = τ = 50.

For each value of (variance/mean) ratio, the inference of the identifiable parameters (see Theorem 4.2) in the case n ≥ 3. was computed on each of the M = 100 repetitions of the K = 300 populations trajectories. Mean and standard deviation of these estimates were then empirically estimated from the M = 100 values obtained. Simulations and statistical analyses were performed using R-8-0 (R Development Core Team, 2005).

Results

. Some technical difficulties occured when we tried to estimate jointly the six identifiable parameters (a,

a ′ b b ′ , b ′ m, b ′ u, bσ, b ′ τ ):
the algorithm we used often did not converge because of the non-linearity in the statistical model. Since we were mainly interested in the estimations for G(.) and µ(.), we assumed that the quantities a and a ′ b b ′ were known. With this simplification, we just had to deal with a linear statistical model; we restricted our attention to the estimation of the parameters where conditionally on (S i , T i ), (S ′ i+1 , R ′ i ) and (S ′′ i+1 , R ′′ i ) are distributed according to two independent Multinomial distributions M(S i ; a, b) and M(T i ; a ′ , b ′ ). Hence, the marginal distribution of S i+1 is,

(b ′ m, b ′ u, bσ, b ′ τ ). b ′ m = 6.5 b ′ u = 40 b ′ σ = 25 b ′ τ = 25 δ 2 /m
P(S i+1 = s ′ /F i ) = (M(S i ; a, b) ⋆ M(T i ; a ′ , b ′ ))(s ′ , r) (B(S i , b) ⋆ B(T i , b ′ ))(r) = p 1 (s ′ /s, t, r).
Joining these results yields that (X i ) is a time homogeneous Markov chain with state space N 5 .

A.2. Proof of Proposition 3.1. We use in the sequel the Kullback-Leibler divergence K(P, Q) of distribution Q w.r.t. P . Recall its definition,

K(P, Q) = - log dQ dP dP = -E P ( dQ dP ) if Q ≪ P ; (A.1) K(P, Q) = +∞ otherwise. (A.2)
This quantity is non-negative and equal to 0 if and only if Q = P P a.s.

Let us first consider l 2 K (c) and l 3 K (d) defined in (3.9), (3.10). The maximum likelihood estimators are

(A.3) ĉK = K k=1 n i=0 v k i K k=1 n i=0 r k i ; dK = K k=1 n i=0 f k i K k=1 n i=0 v k i .
Since the K populations are independent, applying the strong law of large numbers yields the strong consistency of ĉK and dK . The random variables

Z k i = n i=0 (V k i -c 0 R k i ) are i.i.d. centered with variance c 0 (1-c 0 )E θ 0 ( n i=0 R i ), so that the Central Limit Theorem yields that, (A.4) as K → +∞ , √ K(ĉ K -c 0 ) D → N (0, c 0 (1 -c 0 ) E θ 0 ( n i=0 R i ) ) under P θ 0 .
The proof concerning dk is similar : dK → d 0 P θ 0 a.s. and

(A.5) √ K( dK -d 0 ) D → N (0, d 0 (1 -d 0 ) E θ 0 ( n i=0 V i )
).

Consider now the estimation of (a, b, a ′ , b ′ ). Let us first check identifiability. Applying the strong law of large numbers to l 6 K (θ; O 0:n ) defined in (3.13), we get that, under P θ 0 , as K → +∞,

1 K l 6 K (θ; X 0:n (K)) → E θ 0 n i=0 log(M(S i ; a, b) ⋆ M(T i ; a ′ , b ′ ))(S i+1 , R i ).
We can express the limit above using the Kullback-Leibler divergence defined in (A.1),

C(θ 0 ) -E θ 0 n i=0 K(M(S i ; a 0 , b 0 ) ⋆ M(T i ; a ′ 0 , b ′ 0 ), M(S i ; a, b) ⋆ M(T i ; a ′ , b ′ )),
where C(θ 0 ) is a constant depending only on θ 0 , and K(P, Q) is the Kullback-Leibler divergence of the two random probability distributions. Each term of the sum above is non positive and equal to 0 if and only if the two distributions are identical a.s. under P θ 0 . It is easy to check, using the first and second moments of these two distributions, that this implies (a, b, a

′ , b ′ ) = (a 0 , b 0 , a ′ 0 , b ′ 0 ). Hence,(a 0 , b 0 , a ′ 0 , b ′ 0 )
is a strict maximum of the limit above, which leads to the identifiability of these parameters.

Consider now the CLS estimators (â K , â′ K ) which minimize J 1 K (a, a ′ ) defined in (3.14). By the strong law of large numbers, under P θ 0 , as K → ∞,

(A.6) 1 K J 1 K (a, a ′ ) → n i=0 E θ 0 ({(a 0 -a)S i + (a ′ 0 -a ′ )T i } 2 ) + A(θ 0 ), where A(θ 0 ) = a 0 (1 -a 0 )E θ 0 ( n i=0 S i ) + a ′ 0 (1 -a ′ 0 )E θ 0 ( n i=0 E θ 0 T i )
only depends on θ 0 . Since S i and T i are non negative random variables, this limit possesses a strict minimum at (a, a ′ ) = (a 0 , a ′ 0 ). Denote by τ M the transposition of a matrix M and set Z the (n × 2) matrix with rows equal to (S i , T i ), τ S the vector (S 0 , . . . , S n ), τ S the vector (S 1 , . . . , S n+1 ), and Z k , S k , Sk , T k their values for population k. Then J 1 K (a, a ′ ) writes,

1 K J 1 K (a, a ′ ) = 1 K i=1,...,K τ Sk -Z k a a ′ Sk -Z k a a ′ (A.7) so that âK â′ K = K k=1 τ Z k Z k -1 K k=1 τ Z k Sk .
As functions of (a, a ′ ), 1 K J 1 K (a, a ′ ) and its limit defined in (A.6) are twice continuously differentiable a.s.. The parameter set Θ is compact, so we just have to control the continuity modulus of the process 1

K J 1 K (a, a ′ ). It is de- fined, for η > 0, by w(K, η) = sup{ 1 K |J 1 K (a 1 , a ′ 1 ) -J 1 K (a 2 , a ′ 2 )| ; (a 1 , a ′ 1 ) -(a 2 , a ′ 2 ) < η}.
By the Cauchy-Schwarz inequality, w(K, η) is bounded by

{ 1 K k,i (2S k i+1 -(a 1 +a 2 )S k i -(a ′ 1 +a ′ 2 )T k i ) 2 } 1/2 { 1 K k,i ((a 1 -a 2 )S k i +(a ′ 1 -a ′ 2 )T k i ) 2 } 1/2 .
The first term is a random variable converging P θ 0 -a.s. to a deterministic positive limit, and is thus bounded in probability. The second term is bounded by the r.v. 2η( 1

K K k=1 n i=0 (S k i + T k i )) 1/2 → 2η (E θ 0 n i=0 (S i + T i )) 1/2 P θ 0 a.s. Hence, as K → +∞, limsup w(K, η) ≤ φ(η)
, where φ(η) → 0 as η → 0. This ensures the consistency of (â K , â′ K ) defined in (A.7) [START_REF] Dacunha-Castelle | Probabilités et Statistiques. 2.Problèmes à temps mobile[END_REF][START_REF] Van Der Vaart | Asymptotic Statistics[END_REF].

Consider now the asymptotic normality. The function (a, a

′ ) → J 1 K (a, a ′ ) being C 2 a.s., the gradient of J 1 K (a, a ′ ) is DJ 1 K (a, a ′ ) = - K k=1 τ Z k Sk -Z k a a ′ .
The 2x2 matrix containing the second derivatives of

J 1 K (a, a ′ ) is ∇J 1 K (a, a ′ ) = 2 K k=1 τ Z k Z k . Using now that (â K , â′ K ) is consistent, a Taylor expansion of DJ 1 K at (a 0 , a ′ 0 ) yields (A.8) 0 = 1 2 √ K DJ 1 K (a 0 , a ′ 0 ) + 1 2 √ K ∇J 1 K a -a 0 a ′ -a ′ 0 + o P (1),
where o P (1) denotes a remainder term that goes to 0 in P θ 0 -probability.

The r.v. vectors S -Z k a 0 a ′ 0 are i.i.d. centered, hence 1 2 √ K DJ 1 K (a 0 , a ′ 0 ) converges in distribution under P θ 0 to a centered Normal distribution with covariance matrix E θ 0 ( τ ZV (θ 0 )Z) where V (θ 0 ) is the (n+1)×(n+1) diagonal matrix with diagonal elements a 0 (1a 0 )S i + a ′ 0 (1a ′ 0 )T i . Moreover, by the strong law of large numbers, 1 2K ∇J 1 K converges a.s. to E θ 0 ( τ ZZ). By the Cauchy-Schwarz inequality, this matrix is invertible and using (A.8) yields

√ K âK -a 0 â′ K -a ′ 0 = E θ 0 ( τ ZZ) -1 1 √ K K k=1 τ Z k (S k -Z k a 0 a ′ 0 ) + o P (1). Therefore, setting Σ 1 (θ 0 ) = (E θ 0 ( τ ZZ)) -1 E θ 0 ( τ ZV (θ 0 )Z)(E θ 0 ( τ ZZ)) -1 , (A.9) √ K âK -a 0 â′ K -a ′ 0 → D N (0, Σ 1 (θ 0 )).
Similarly, setting τ R the vector (R 0 , . . . , R n ), define (A.10) bK b′

K = K k=1 τ Z k Z k -1 K k=1 τ Z k R k .
The proof concerning the estimation of (b, b ′ ) is similar; so we get, setting W (θ 0 ) diagonal matrix with diagonal elements b 0 (1

-b 0 )S i + b ′ 0 (1 -b ′ 0 )T i and Σ 2 (θ 0 ) = (E θ 0 ( τ ZZ)) -1 E θ 0 ( τ ZW (θ 0 )Z)(E θ 0 ( τ ZZ)) -1 , (A.11) √ K bK -b 0 b′ K -b ′ 0 → D N (0, Σ 2 (θ 0 )).
Finally, let us study the estimation of (m, u) based on the CLS process J 4 K (m, u) defined in (3.16). Under the assumptions (A1), (A3), G(.) and µ(.) have finite variances δ 2 and ρ 2 . Denote by δ 2 0 and ρ 2 0 the variances associated with θ 0 . Then, as K → ∞, under P θ 0 ,

1 K J 4 K (m, u) → E θ 0 n i=0 (δ 2 0 F i + ρ 2 0 ) + E θ 0 n i=0 [(m 0 -m)F i + (u -u 0 )] 2 .
Clearly, the above functional has a strict minimum at (m 0 , u 0 ), leading to the identifiability of (m, u). Let τ F denote the vector (F 0 , . . . , F n ), G the (n + 1) × 2 matrix with rows equal to (F i , 1), τ T the vector (T 1 , . . . , T n+1 ), F k , G k , T k their values in population k. Then, the conditional least square estimator is

(A.12) mK ûK = K k=1 τ G k G k -1 K k=1 τ G k T k .
Consistency and asymptotic normality are obtained with a proof similar to the one detailed above. Define W ′ (θ 0 ) the diagonal matrix with diagonal elements δ 2 0 F i + ρ 2 0 , and

Σ 3 (θ 0 ) = (E θ 0 ( τ GG)) -1 E θ 0 ( τ GW ′ G)(E θ 0 ( τ GG)) -1 , then (A.13) √ K mK -m 0 ûK -u 0 → D N (0, Σ 3 (θ 0 )) in distribution under P θ 0 .
Using now that the likelihood splits into five terms that can be maximized separately leads to the asymptotic independence of the estimators stated in Proposition 3.1.

A.3. Proof of Proposition 4.1. Let us prove by induction Proposition 4.1. By Assumption A6, the property holds for i = 0 : S 0 and T 0 are independent and S 0 ∼ P(σ) = P(Γ 0 ), T 0 ∼ P(τ ) = P(Γ ′ 0 ). Assume that the property holds for i ≥ 1, i.e. using definitions (4.1), (4.3), S i and T i are independent conditionnally on G i-1 and L(S i /G i-1

) = P(Γ i ), L(T i /G i-1 ) = P(Γ ′ i ). Using (2.5), (S i+1 , R i ) = (S ′ i+1 , R ′ i ) + (S ′′ i+1 , R ′′ i ), where the conditional dis- tribution of (S ′ i+1 , R ′ i ) (resp.(S ′′ i+1 , R ′′ i )) given (S i , T i ) is the Multinomial distribution M(S i ; a, b) (resp. M(T i ; a ′ , b ′ )).
Using Assumption (A2), these two distributions are independent conditionally on S i , T i and applying now Lemma 4.1 conditionally on G i-1 to {S i ∼ P(Γ i ), M(S i ; a, b)} and to

{T i ∼ P(Γ ′ i ), M(T i ; a ′ , b ′ )} yields that the four variables S ′ i+1 , R ′ i , S ′′ i+1 , R ′ i are in- dependent conditionally on G i-1 and that S ′ i+1 ∼ P(aΓ i ), S ′′ i+1 ∼ P(a ′ Γ ′ i-1 ), R ′ i ∼ P(bΓ i ) and R ′′ i ∼ P(b ′ Γ ′ i ).
Hence, S i+1 and R i are independent conditionally on G i-1 and S i+1 ∼ P(aΓ i + a ′ Γ ′ i ), and R i ∼ P(bΓ i + b ′ Γ ′ i ). Let us now prove that E(S i+1 /G i ) = E(S i+1 /G i-1 ). Let φ 1 and φ 2 two measurable functions of (Y 0:i-1 , S i+1 ) and (Y 0:i-1 , R i , V i , F i ) and compute,

E G i-1 (φ 1 (S i+1 )φ 2 (R i , V i , F i )) = E G i-1 ( φ 1 (s i+1 )φ 2 (r i , v i , f i )ds i+1 dr i dv i df i ).
Using (2.13) and (2.14) , set ψ 2 (r) = φ 2 (r, v, f )p 4 (v/r)p 5 (f /v)dv, then

E G i-1 (φ 1 (S i+1 )φ 2 (R i , V i , F i )) = E G i-1 (φ 1 (S i+1 )ψ 2 (R i )) = E G i-1 (φ 1 (S i+1 ))E G i-1 (ψ 2 (R i )) = E G i-1 (φ 1 (S i+1 ))E G i-1 (φ 2 (R i , V i , F i )),
since S i+1 and R i are independent conditionally on G i-1 . Hence, S i+1 and (R i , V i , F i ) are independent conditionally on G i-1 and,

E(S i+1 /G i ) = E(S i+1 /G i-1 ) = P(aΓ i + a ′ Γ ′ i+1 ).
Consider two measurable functions φ 3 , φ 4 of (Y 0:i , S i+1 ) and (Y 0:i , T i+1 ), then using (2.12), (2.14), Assumption (A2) and the conditional independence given

G i-1 of S i+1 and (R i , V i , F i ) yield, (A.14) E G i (φ 3 (S i+1 )φ 4 (T i+1 )) = E G i ( (φ 3 (s i+1 )φ 4 (t ′ )p 2 (t ′ /F i )ds i+1 dt ′ ) = E G i (φ 3 (S i+1 ))( φ 4 (t ′ )p 2 (t/F i )dt ′ ) = E G i (φ 3 (S i+1 ))E G i (φ 4 (T i+1 )).
Therefore, conditionally on G i S i+1 and T i+1 are independent and using now Assumption (A5), L(T i+1 /G i ) = P(mF i + u). The property holds for i + 1 with Γ i+1 = aΓ i + a ′ Γ ′ i and Γ ′ i+1 = mF i + u, which are the definitions given in (4.3).

A.4. Proof of Theorem 4.2. The likelihood l 1 K (θ) defined in (4.19) sums up the available information associated with the incomplete model. An application of the strong law of large numbers yields,

1 K l 1 K (θ) → L(θ 0 , θ) = E θ 0 n i=0 (-Λ i (θ, Y 0:i-1 ) + R i log Λ i (θ, Y 0:i-1 )).
According to Theorem 4.1, conditionally on G i-1 , the random variables R i are Poisson distributions with parameter Λ i (θ 0 , Y 0:i-1 ) under P θ 0 . Hence,

L(θ 0 , θ) = E θ 0 n i=0 (-Λ i (θ, Y 0:i-1 ) + Λ i (θ 0 , Y 0:i-1 ) log Λ i (θ, Y 0:i-1 )).
Using the explicit form of the Kullback-Leibler divergence between two Poisson distributions, L(θ 0 , θ) writes,

L(θ 0 , θ) = - n i=0 E θ 0 {K(P(Λ i (θ 0 , Y 0:i-1 )), P(Λ i (θ, Y 0:i-1 )))} + C(θ 0 ), with C(θ 0 ) = E θ 0 n i=0 (-Λ i (θ 0 , Y 0:i-1 )+Λ i (θ 0 , Y 0:i-1 ) log Λ i (θ 0 , Y 0:i-1
)) only depends on the observations. The identifiability condition for θ 0 is therefore equivalent to {L(θ 0 , θ) = 0} ⇒ {θ = θ 0 }.

The Kullback-Leibler divergence of two Poisson distributions K(P(µ 0 ), P(µ)) is non negative and equal to 0 if and only if µ = µ 0 . Hence, the limit L(θ 0 , θ) presents a strict maximum at θ 0 if and only if Λ i (θ, Y 0:i-1 ) = Λ i (θ 0 , Y 0:i-1 ) P θ 0a. s. for i = 0, . . . , n.

Since the Λ i (θ, Y 0:i-1 ) depend on (F i-1 , . . . , F 0 ), the above condition can hold only if the coefficients associated with the random variables F i in the above expression are equal. The proof below is just elementary algebra based on this property. We obtain, using for the c i (θ) the definitions given in (4.7) and (4.15) :

If n = 0, Λ 0 (θ) = c 0 (θ) is deterministic. Only c 0 (θ) is identifiable: it is the first condition stated in Theorem 4.2. If n = 1, we have, Λ 1 (θ, Y 0:i-1 ) = b ′ mF 0 + c 1 (θ). Since F 0 is random, the two random variables Λ 1 (θ, Y 0:i-1 ) and Λ 1 (θ 0 , Y 0:i-1 ) are P θ 0 -a.s. equal iff b ′ m = b ′ 0 m 0 and c 1 (θ) = c 1 (θ 0 ), which leads to the identifiability of (b ′ m, c 0 (θ), c 1 (θ)). If n = 2, Λ 2 (θ) = b ′ m F 1 + a ′ b b ′ b ′ mF 0 + c 2 (θ),
We thus get two additional conditions which lead to the identifiability of (

a ′ b b ′ , b ′ m, c 0 (θ), c 1 (θ), c 2 (θ)). If n = 3, Λ 3 (θ) = b ′ mF 2 + a ′ b b ′ b ′ m(aF 1 + F 0 ) + c 3 (θ)
. Hence, F 0 , F 1 , F 2 being random variables, we get that (a, a ′ b b ′ , b ′ m) are identifiable. Now, identifying (bσ, b ′ τ, b ′ u) consists in solving a linear system using the conditions on c 0 (θ), c 1 (θ), c 2 (θ)). We obtain :

-

if {a = a ′ b b ′ }, (a, a ′ b b ′ , b ′ m, b ′ u, bσ, b ′ τ ) is identifiable. -if a = a ′ b b ′ , the identifiable parameters are (a, b ′ m, b ′ u, bσ + b ′ τ ). Noting that, for i ≥ 1, c i+1 (θ) -ac i (θ) = (1 -a + a ′ b b ′ )b ′ u and that only (a, b ′ m, a ′ b b ′ ) enters in the F i 's coefficients.
it can be checked that observing more generations does not lead to the idenfiability of parameters.

A.5. Proof of Proposition 4.2. Let φ = (a, a ′ b b ′ , b ′ m, b ′ u, bσ, b ′ τ ) with a = a ′ b b ′ . Under Asssumption (A7), we can define K 1 , K 2 , a 1 , a 2 such that Φ ⊂ [K 1 , K 2 ] 5 × [a 1 , a 2 ], with 0 < K 1 < K 2 < +∞ and 0 < a 1 < a 2 < 1.
The normalized loglikelihood process can be written as

(A.15) 1 K l 1 K (φ, Y 0:n (K)) = 1 K K k=1 J(φ, Y k 0:n ) with J(φ, Y k 0:n ) = n i=0 R k i log Λ i (φ, Y k 0:i-1 ) -Λ i (φ, Y k 0:i-1 ).
These random variables are i.i.d. and we have to study the behaviour of their empirical distribution with respect to φ. For getting the consistency of the associated maximum likelihood estimator, we have to prove that the parametric class {P φ , φ ∈ Φ} is Glivenko-Cantelli (see e.g. [START_REF] Van Der Vaart | Asymptotic Statistics[END_REF]. There is no close form for the density of these variables, so that no generic argument can here be applied: we thus propose a direct proof. The functional {φ → J(φ, Y 0:n )} is a.s twice continuously differentiable on Φ. Let DΛ i (φ) denote the gradient in R 6 of Λ i (φ, Y 0:i-1 ). Then,

| J(φ ′ ; Y 0:n )-J(φ ′′ ; Y 0:n ) |≤ n i=0 R i sup φ∈Φ ( 1 Λ i (φ) +1)DΛ i (φ) φ ′ -φ ′′ .
Under (A7), we have that, for all i, Λ i (φ) ≥ K 1 > 0. Let us compute DΛ i (φ) Using (4.4) for the definitions of c i (θ) and Λ i (θ), let us set

Λ i (φ) = b ′ mF i-1 + a ′ b b ′ (b ′ m)(F i-2 + aF i-3 + • • • + a i-2 F 0 ) + c i (φ). For i = 0, we have ∂Λ 0 (φ) ∂a = ∂Λ 0 (φ) ∂(a ′ b/b ′ ) = ∂Λ 0 (φ) ∂(b ′ m) = ∂Λ 0 (φ) ∂(b ′ u) = 0; ∂Λ 0 (φ) ∂(bσ) = ∂Λ 0 (φ) ∂(b ′ τ ) = 1.
Using the convention that non properly defined terms are set to 0 (e.g. the sum (F i-3 + 2aF i-4 + . . . ) is set to 0 for i ≤ 2, we get, for i ≥ 1,

∂Λ i (φ) ∂a = a ′ b b ′ (b ′ m)(F i-3 + 2aF i-4 + • • • + (i -2)a i-3 F 0 ) + ∂c i (φ) ∂a , with
∂c i (φ) ∂a = (bσ)ia i-1 + a ′ b b ′ (b ′ u) d da ( 1 -a i-1 1 -a ) + a ′ b b ′ (b ′ τ )(i -1)a i-2 ; ∂Λ i (φ) ∂(a ′ b/b ′ ) = b ′ m(F i-2 + aF i-3 + • • • + a i-2 F 0 ) + (b ′ τ )a i-1 + b ′ u 1 -a i-1 1 -a ; ∂Λ i (φ) ∂(b ′ m) = F i-1 + a ′ b b ′ (F i-2 + aF i-1 + • • • + a i-2 F 0 ) ; ∂Λ i (φ) ∂(b ′ u) = 1 + a ′ b b ′ 1 -a i-1 1 -a ; ∂Λ i (φ) ∂(bσ) = a i ; ∂Λ i (φ) ∂(b ′ τ ) = a ′ b b ′ a i-1 .
All these partial derivatives are positive for i ≥ 1 and, except the first one

∂Λ i (φ)
∂a , they are bounded from above by M 1 (1 + i-1 i=0 F j ) where M 1 is a constant determined by Φ (since on Φ, 0 < a 1 ≤ a ≤ a 2 < 1). Noting now that the application {i → ia i-1 } satisfies ∀i ≥ 1, α 1 < ia i-1 < α 2 on [a 1 , a 2 ] ⊂ (0, 1) with 0 < α 1 < α 2 < +∞, we can also bound ∂Λ (1 + i-1 j=0 F j )R i .

Using now that E φ 0 (R i i-1 j=0 F j ) = E φ 0 (Λ i (φ 0 ) i-1 j=0 F j ), Z n satisfies,

E φ 0 Z n ≤ nE φ 0 n i=0 (1 + F j ) 2 ≤ n n i=0 1 + E φ 0 F 2 i 2 .
This is finite since G(.) and µ(.) have finite variance and n is prescribed. Let us define the r.v.

Z k n = n i=0 (1 + i-1 j=0 F k j )R k i .
Then the continuity modulus of l1 K (θ) verifies :

w(K, η, 1 K l1 K ) = sup φ ′ -φ ′′ ≤η 1 K K k=1 | J(φ ′ ; Y k 0:n ) -J(φ ′′ ; Y k 0:n ) | ≤ η 1 K K k=1 Z k n .
Using now that E φ 0 Z n < ∞, we can apply the strong law of large numbers to w(K, η, 1 K l1 K ), which is a sufficient condition for ensuring the consistency of φK .

It remains to study the asymptotic normality of the estimators. It is easy to check that the random variables τ W k = (W k 1 , . . . , W k 6 ) with

W k p = n i=0 ( R k i Λ k i (φ 0 ) -1)
∂Λ k i (φ 0 ) ∂φ p are i.i.d. centered under P φ 0 , with covariance matrix Σ is for 1 ≤ p, q ≤ 6 Σ p,q = Σ p,q (φ 0 ) = E φ 0 n i=0

1 Λ i (φ 0 ) ∂Λ i ∂φ p (φ 0 ) ∂Λ i ∂φ q (φ 0 ).
This matrix is well-defined and finite is since, for all i, Λ i (φ) ≥ K 1 > 0 and | ∂Λ i (φ) ∂φp |≤ M 3 (1 + i-1 0 F j ) for p = 1, . . . 6. Now, φK is a zero for the score function, so that a Taylor expansion at φ 0 yields, using the consistency of the vector φK ,

0 = 1 √ K D l1 K ( φK ) = 1 √ K D l1 K (φ 0 )+ 1 K ∇ l1 K (φ 0 )+R K (φ 0 , φK √ K( φK -φ 0 )
The term ∇ l1 K (φ 0 ) contains the second derivatives of l1 K (φ) w.r.t. φ p , φ q and RK (φ 0 , φK ) is the remainder term of the Taylor expansion. So, the strong law of large numbers yields :

1 K ∇ l1 K (φ 0 ) → -n i=0 E φ 0 ( 1 Λ i (φ 0 ) ∂Λ i (φ 0 ) ∂φp ∂Λ i (φ 0 )
∂φq ). The remainder term R K (φ 0 , φK ) is bounded uniformly on Φ by φKφ 0 Z K , with Z K = sup{ 1 K ∇ l1 K (φ), φ ∈ Φ}. Using that Z K is bounded uniformly on Φ by n 2 M 6 (1 + n i=0 F i F j ) 2 , we get that R K (φ 0 , φ) goes to 0 under P φ 0 , which leads to the result stated in Proposition 4.2 provided that Σ(φ 0 ) is invertible.

Aknowledgments

This study was funded by the European project "sustainable Introduction of Genetically Modified Crops into European Agriculture" and by the project "Flux des (trans-)gènes et impact sur la biodiversité" of the Agence nationale de la Recherche.

Aurélie Garnier Université Paris Sud, Laboratoire Ecologie, Sytématique & Evolution, UMR8079, F-91405 Orsay, France;CNRS, 91405 0rsay, France; AgroParisTech, 75231 Paris, France. E-mail: au.garnier@gmail.com

Fig 1 .

 1 Fig 1. Schematic Dynamics of feral oilseed rape populations.

First

  consider one population and omit the index k in what follows. Set (4.1)

  Assume now that n ≥ 3, a = a ′ b b ′ and let us study the inference of the identifiable parameters. Let us denote by φ = (a, a ′ b b ′ , b ′ m, b ′ u, bσ, b ′ τ ) (resp. φ 0 ) an arbitrary (resp. the true value) of the paramete, by Φ the parameter set and asssume Assumption 7. Φ is a compact set of R * + 5 × (0, 1); and φ 0 ∈ • Φ.

  ) and standard deviation (sd) of estimators when the offspring G(.) is a binomial distribution with mean m and variance δ 2 instead of a P(m) distribution. Immigration µ ∼ P(u); a = 0.16, a ′ = 0.006, b = b ′ = 0.5, m = 13, u = 80, σ = τ = 50. Value b ′ m = 6.5 b ′ u = 40 b ′ σ = 25 b ′ τ = 25 ρ 2 /

  ) and standard deviation (sd) of estimators when the immigration µ(.) is a negative binomial distribution with mean u and variance ρ 2 instead of a P(u) ditribution.Offspring G(.) ∼ P(m); a = 0.16, a ′ = 0.006, b = b ′ = 0.5, m = 13, u = 80, σ = τ = 50.

  i

F

  j ) , | J(φ; Y 0:n ) -J(φ ′ ; Y 0:n ) |≤ ηM 3 Z n with Z n = n i=0

The results are given in the two tables below. Concerning deviations of G(.) from the distribution P(m), we obtained that the estimation procedure performed very well, even for large deviations from the Poisson case : the bias remained less than 5% for all the four estimated quantities with a variance/ratio up to 1000. (see Table 5.2) When immigration was assumed to follow a Negative Binomial distribution, the estimation procedure performed quite well for values of (variance/mean) ratios up to 50 : the biases remained less than 10% for all four identifiable quantities. For larger values of variance/mean ratio, the bias could 17% of parameter value approximately.

APPENDIX A: APPENDIX SECTION

A.1. Proof of Proposition 2.1 . Consider first the initial distribution π 0 (x). Successive conditionings yield

The first distribution is ν(s, t). Using definitions (2.4), (2.13), the last two conditional distributions are equal to p 4 (v/r) and

, leading to Binomial distributions B(S 0 ; b) for the marginal distribution of R ′ 0 (resp. B(T 0 ; b ′ ) for R ′′ 0 ). By (A2), these two distributions are independent conditionally on (S 0 , T 0 ) which yields (2.12) and the expression of π 0 (x). Let us now study two successive generations. Using notations (2.7), the conditional distribution of X i+1 given F i can be expressed for x = (s, t, r, v, f ),

The last three conditional distributions have already been computed for getting π 0 (x). They are respectively equal to p 3 (r ′ /s ′ , t ′ ), p 4 (v ′ /r ′ ) and p 5 (f ′ /v ′ ) defined in (2.12), (2.13). Let us compute P(T i+1 = t ′ /s ′ ; F i ). Seeds on the ground at cycle (i + 1) come from two sources, offspring of mature plants F i and seed immigration during cycle i + 1, T i+1 = T ′ i+1 + I i+1 . Using (A1), the distribution of T ′ i+1 is G ⋆F i . By (A3), I i+1 is independent of S i+1 , T ′ i+1 and F i , so that P(T i+1 = t ′ /s ′ ; F i ) = (G ⋆f ⋆ µ)(t ′ ) = p 2 (t ′ /f ). The last distribution is P(S i+1 = s ′ /F i ). Using (A2) and (2.5) yields that S i+1 = S ′ i+1 + S ′′ i+1