Yuting Chen

Samis Trevezas
email: samis.trevezas@ecp.fr

Paul-Henry Cournède

A regularized particle filter EM algorithm based on Gaussian randomization with an application to plant growth modeling

Keywords: Gaussian randomization, stochastic EM algorithm, regularized particle filter, plant growth model, LNAS model, state space model

 can be used in this frequentist-based approach to perform the Monte Carlo E-step at each iteration of the stochastic EM algorithm. This regularized version is particularly adapted to the framework of Gaussian randomization since the last iterations of the EM algorithm are characterized by low variance in the parameter distributions. A toy example with available analytic solutions, a synthetic example and a real data application with scarce observations to the LNAS (Log-Normal Allocation and Senescence) model of sugar beet growth are presented to highlight some theoretical and practical aspects of the proposed methodology.

Introduction

The current paper is motivated by the need to have a simpler method to perform parameter estimation in plant growth models which are usually characterized by complex parameterization and scarce experimental data. Nevertheless, our proposed method is generic and could be used in any statistical model under certain conditions that we state in the sequel. The basic idea consists in interpreting some of the parameters of the model as the means of unobserved states following normal distributions with unknown variances.

The resulting data and parameter augmentation method enable us to estimate more easily the parameters of the initial model with the help of an iterative state estimation technique for the extended model. The term Gaussian randomization is introduced for such type of modifications. The proposed method should not be confused with the Bayesian approach. The variance parameters involved in the Gaussian randomization are assumed to be unknown in contrast to the Bayesian approach in which normal priors are selected for the parameters of the model.

The estimation is performed in the classical framework of incomplete data models via an appropriate stochastic variant of an EM-type algorithm [START_REF] Dempster | Maximum Likelihood from incomplete data via the EM algorithm[END_REF]. In contrast to classical stochastic variants ( [START_REF] Celeux | The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem[END_REF], [START_REF] Wei | A Monte Carlo implementation of the EM algorithm and the poor man's data augmentation algorithms[END_REF], [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF]) the proposed regularized particle-filter EM algorithm is particularly adapted to the extended model resulting from the Gaussian randomization, since the last iterations of the algorithm are performed under very low variance scenarios. Some recent efforts to obtain appropriate stochastic EM-type algorithms for complex models arising in plant growth model applications can be found in [START_REF] Trevezas | A sequential Monte Carlo approach for MLE in a plant growth model[END_REF] and [START_REF] Trevezas | Simulation techniques for parameter estimation via a stochastic ECM algorithm with application to plant growth modeling[END_REF]. These algorithms are designed for a non-explicit M-step and a numerical maximization algorithm is involved as well in one of the conditional maximization steps. In this paper we present a simple way of turning a non-explicit M-step into an explicit one and opt for the postregularized particle filter (post-RPF, e.g., [START_REF] Musso | Regularization schemes for branching particle systems as a numerical solving method of the nonlinear filtering problem[END_REF], [START_REF] Oudjane | Multiple model particle filter[END_REF] and [START_REF] Gland | Stability and uniform approximation of nonlinear filters using the Hilbert metric and application to particle filters[END_REF]) to perform the E-step. We mention that the use of this type of particle filter was restricted until now to Bayesian type estimation.

The method is proved to be well-adapted for mathematical models of plant growth and more generally in life sciences, for which model parameterization is generally a difficult process. Such type of models are characterized by a large number of interacting processes with a large number of model parameters, nonlinear dynamics and scarce data for parameter estimation, resulting from costly experimental data acquisition. Bayesian inference methods offer interesting perspectives to cope with such characteristics, mostly due to their robustness when confronted to scarce data, and have been used for this pur-pose in plant growth modeling ( [START_REF] Wyckoff | Predicting tree mortality from diameter growth: a comparison of maximum likelihood and bayesian approaches[END_REF], [START_REF] Gaucherel | Parameterization of a process-based tree growth model: Comparison of optimization, mcmc and particle filtering algorithms[END_REF]). However, in such situations the influence of a priori distributions is prevalent ( [START_REF] Chen | Assessment of parameter uncertainty in plant growth model identification[END_REF]) while it is generally not easy to assess precisely.

Our approach can be seen as an effort to combine the robustness of Bayesian methods while preserving the prevalence of observation data and the interest of point estimation for mechanistic plant models, in which the biophysical meaning and value of parameters are crucial.

The paper is organized as follows. In Section 2, we introduce the notion of Gaussian randomization and present a theoretical framework for transforming the parameter estimation problem in the initial model into a problem that can be solved with the help of an appropriate EM-type algorithm in the resulting incomplete data model. In Section 3, we illustrate this idea in a toy Gaussian model, where analytic solutions are available. In Section 4, we first describe a parametric model which arises in the description of the sugar beet growth, then give a state space model representation and finally apply the idea of Gaussian randomization for this model. In Section 5, we introduce a regularized particle filter approximation of the E-step of the EM algorithm in order to make feasible parameter estimation and give a synthetic example. In order to illustrate the performance of the algorithm with real data, we present in Section 6 a case study with irregular observations (scarce data) from the sugar beet plant.

Gaussian randomization

Let us consider an arbitrary parametric statistical model m {(Ω, A, P µ ); µ ∈ M }, where M is a euclidean subset, and a random vector Y (representing the data vector) defined on this probability space and taking values in a measured space {(Y, B, ν)}, where ν is a reference measure. In this paper the vector Y admits a density p(y; µ) w.r.t. ν for each µ ∈ M and for a given observation Y = y, the likelihood function is denoted by L(µ).

The following assumption will be used in the sequel:

Assumption 1. i) The likelihood function L(µ) is continuous on M ,
ii) The model has a unique MLE µ * , i.e., for all µ ∈ M,

0 ≤ L(µ) ≤ L(µ * ),
and the second inequality is strict for µ = µ * .

Let us suppose that maximization of L(µ) or log L(µ) w.r.t. some of the parameters of µ is very difficult to perform. In particular, let us decompose the unknown parameter µ = (µ 1 , µ 2 ), where µ 1 could represent a part of the parameter vector that could not be updated explicitly in an iterative conditional maximization procedure by fixing µ 2 . The general idea consists in considering an enlarged model, by adding parameters and hidden variables to the existing model, in such a way that maximization w.r.t. to the augmented parameter vector in the enlarged model is equivalent to the initial maximization problem.

In this direction, we give a first definition. Let d x denote the dimension of a vector x. We treat the case where µ 1 ∈ R dµ

p(y|θ; µ 2 ) = L(θ, µ 2 ).
Let us take σ 2 = (σ 2 i ) 1≤i≤dµ 1 ∈ (R * + ) dµ 1 to be a minimal vector representation of the corresponding covariance matrix Σ. The parameterization corresponding to m(µ 1 ) is given by:

φ = (µ, σ 2 ) = (µ 1 , µ 2 , σ 2 ) ∈ Φ = M × (R * + ) dµ 1 ⊂ R d φ . (1)
Let us now assume that we allow some variance parameters of Θ to be null. This is equivalent to say that the corresponding parameters are not randomized. Denote by µ 11 the subvector of µ 1 with associated strictly positive variances. According to Definition 1-i) this model is no longer a Gaussian randomization of m w.r.t. µ 1 but rather a Gaussian randomization of m w.r.t. µ 11 . In order to treat in a common framework the class of all Gaussian randomizations of m w.r.t. to subvectors of µ 1 we allow null variances in the parameterization. Most importantly, if σ 2 is assigned to the null vector 0 dµ 1 , then the parameter φ = (µ 1 , µ 2 , 0 dµ 1 ) could be identified with the parameter µ = (µ 1 , µ 2 ) of the initial model m. In this way m could be understood as a submodel of m(µ 1 ). In the sequel we justify theoretically this intuition. Let us first state the following lemma. Proof: First, we compute the likelihood function L(φ) of the extended model m(µ 1 ). We have:

L(φ) = p(y; φ) = θ∈M 1 p(θ, y; φ)dθ = θ∈M 1 p(θ; µ 1 , σ 2 ) p(y|θ; µ 2 )dθ = E µ 1 ,σ 2 p(y|Θ; µ 2 ) , (2)
where Θ ∼ N dµ 1 (µ 1 , Σ) according to Definition 1-i). By using Definition 1-ii) and Assumption 1-ii) we get successively that for any φ ∈ Φ,

L(φ) = E µ 1 ,σ 2 L(Θ, µ 2 ) ≤ L(µ * ), (3)
and this completes the proof.

In order to justify the existence of the MLE in this extended model, we need to extend the likelihood function of m(µ 1 ) to parameter values which include null variances as explained at the beginning of this section and give the correct interpretation in this type of boundary values. We will need the following lemma.

Lemma 2. Let x ∈ R dx , I ⊂ {1, . . . , d x } and J = {1, . . . , d x } \ I. Let also σ 2 I and σ 2
J be the subvectors of σ 2 ∈ R dx with components indexed by I and J respectively. If we consider the family of non-singular random vectors {X σ 2 ; X σ 2 ∼ N dx (x, Σ)}, where

Σ = diag{σ 2 i } 1≤i≤dx , then we have X σ 2 weakly -----→ ||σ 2 J ||→0 X σ 2 I ∼ N dx (x, Σ I ),
where Σ I is a singular covariance matrix which results from Σ by putting σ 2 J = 0 J . In particular, if I = ∅, then

X σ 2 weakly ----→ ||σ 2 ||→0
x.

Proof: We denote by φ σ 2 , φ σ 2 I and φ 0 the characteristic functions of X σ 2 , X σ 2 I and the constant vector x respectively. It suffices to show that as σ 2 J → 0 J , φ σ 2 converges pointwise to φ σ 2 I . Indeed, let t = (t 1 , . . . , t dx ) ⊤ ∈ R dx . We have

φ σ 2 (t) = e it ⊤ x-1 2 t ⊤ Σt -----→ ||σ 2 J ||→0 e it ⊤ x-1 2 t ⊤ Σ I t = φ σ 2 I (t). (4)
In particular, if I = ∅, then Σ I is the null matrix and by (4), φ σ 2 converges pointwise to

φ 0 .
Proposition 1. Let I ⊂ {1, . . . , d µ 1 } and J = {1, . . . , d µ 1 }\I. Let also µ I be the subvector of µ 1 with components indexed by I, and σ 2 I , σ 2 J and Σ I as given in Lemma 2. We have

L(φ) = L(µ, σ 2 I , σ 2 J ) -----→ ||σ 2 J ||→0 L(µ, σ 2 I ),
where L(µ, σ 2 I ) is the likelihood function associated with the m(µ I ) Gaussian randomization of m. In particular, if I = ∅, then

L(φ) = L(µ, σ 2 ) ----→ ||σ 2 ||→0 L(µ),
where L(µ) is the likelihood function associated with the initial model m.

Proof: By (3) the likelihood L(φ) is expressed as the expectation of the random variable L(Θ, µ 2 ), where Θ ∼ N dµ 1 (µ 1 , Σ). For clarity, let us denote {Θ σ 2 } the family of random vectors Θ indexed by σ 2 . Then, we can rewrite (3) as

L(φ) = E µ 1 ,σ 2 L(Θ σ 2 , µ 2 ) . (5)
If we apply Lemma 2 for the family {Θ σ 2 }, we have that for any (µ 1 , σ 2 I ),

Θ σ 2 weakly -----→ ||σ 2 J ||→0 Θ σ 2 I ∼ N dµ 1 (µ 1 , Σ I ),
or equivalently, for any continuous and bounded function f : R

dµ 1 → R, E µ 1 ,σ 2 f (Θ σ 2 ) -----→ ||σ 2 J ||→0 E µ 1 ,σ 2 I ,0 J f (Θ σ 2 I ) , (6)
where we have used the characterization of weak convergence of a family of random vectors (indexed by a continuous parameter) to a random vector in a specific limit point.

By Assumption 1, the likelihood function L(µ) is continuous on M and bounded by L(µ * ).

Consequently, for any µ 2 , the function g : R dµ 1 → R, where µ 1 → L(µ 1 , µ 2 ) is continuous on M 1 and bounded. Using this and (6), we have

E µ 1 ,σ 2 L(Θ σ 2 , µ 2 ) = E µ 1 ,σ 2 g(Θ σ 2 ) -----→ ||σ 2 J ||→0 E µ 1 ,σ 2 I ,0 J g(Θ σ 2 I ) = E µ I ,σ 2 I L(Θ I,σ 2 I , µ J , µ 2 ) , (7)
where Θ I,σ 2 I is the subvector of Θ σ 2 indexed by I. In the last expectation σ 2 i > 0 only for i ∈ I and consequently the last term corresponds to the analogue of ( 5) for an equivalent representation of the likelihood function associated to the m(µ I ) Gaussian randomization of m. If I = ∅, then obviously the last term coincides with L(µ) and the proof is complete.

By Proposition 1 we can consider in a unified framework the class { m(µ I )} I of 2 dµ 1 -1 Gaussian randomizations of m w.r.t. to subvectors of µ 1 together with the initial model m, which corresponds to the choice I = ∅. For this reason we extend the domain of L(φ), from Φ given by (1) to Φ = M × (R + ) dµ 1 by allowing null variances as follows:

L(µ, σ 2 I , 0 J ) = L(µ, σ 2 I ), (8)
L(µ, 0 dµ 1 ) = L(µ).

With similar arguments as in the proof of Proposition 1 we can deduce from the continuity of L(µ) on M that L(φ) is continuous on Φ (as (µ, σ 2 ) → (µ 0 , σ 2 0 ) ∈ Φ show that (Θ µ 1 ,σ 2 , µ 2 ) converges weakly to (Θ µ 1,0 ,σ 2 0 , µ 2,0 )). Consequently, both equations ( 8) and (9) determine the continuous extension of L from Φ to Φ by using Proposition 1.

φ * = (µ * , 0 dµ 1 ), (10)
L(φ * ) = L(µ * ). ( 11
)
Conversely, if Assumption 1-i) holds and the MLE φ * associated with the model m µ 1 exists and is unique, then it determines a unique MLE µ * associated with the model m and satisfies [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF] and [START_REF] Dempster | Maximum Likelihood from incomplete data via the EM algorithm[END_REF].

Proof: ⇒) The restriction of L(φ) on Φ = M × (R * + ) dµ 1 is upper bounded by L(µ * ) by Lemma 1. Now, if φ ∈ Φ \ Φ, take the nonempty set J = {j ∈ {1, . . . , d µ 1 } : σ 2 j = 0} and I the complement of J. Then, L(φ) can be expressed as the left-hand member of equations [START_REF] Cournède | Development and evaluation of plant growth models: Methodology and implementation in the PYGMALION platform[END_REF] or [START_REF] Moral | Non-linear filtering: interacting particle solution[END_REF]. By the same equations we obtain that L(φ) equals the likelihood function of the extended model m(µ I ) or of the initial model m respectively. Therefore, by applying Lemma 1 for µ I or by Assumption 1 respectively, we infer that L(φ) is upper bounded by L(µ * ) also on Φ \ Φ, and consequently on Φ. By [START_REF] Moral | Non-linear filtering: interacting particle solution[END_REF] we have that L(µ * , 0 dµ 1 ) = L(µ * ) and consequently the MLE exists and if it is unique it satisfies [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF] and [START_REF] Dempster | Maximum Likelihood from incomplete data via the EM algorithm[END_REF]. We have that (µ * , 0 dµ 1 ) ∈ M × {0} dµ 1 and by Assumption 1 and ( 9), it is unique on this set. Let us now assume that there exists φ = (µ,

σ 2 ) ∈ M × R dµ 1 + \ {0} dµ 1 such that L(µ, σ 2 ) = L(µ * ). Take the nonempty set I = {i ∈ {1, . . . , d µ 1 } : σ 2
i > 0} and J the complement of I. By ( 8), (3), the remark following [START_REF] Chen | Assessment of parameter uncertainty in plant growth model identification[END_REF] and our assumption, we have

L(µ, σ 2 ) = E µ I ,σ 2 I L(Θ I , µ ′ 2 ) = L(µ * ), (12)
where µ ′ 2 concatenates the components of µ 1 indexed by J and µ 2 . Since by Assumption 1, the function L(µ) is upper bounded by L(µ * ), we obtain that the random variable

L(Θ I , µ ′ 2 ) has the same upper bound, that is, L(µ * ) -L(Θ I , µ ′ 2
) is a positive random variable. But, by ( 12)

E µ I ,σ 2 I L(µ * ) -L(Θ I , µ ′ 2 ) = 0, so we get L(Θ I , µ ′ 2 ) = L(µ * ), P µ I ,σ 2 I -a.s. , or 1 = P µ I ,σ 2 I L(Θ I , µ ′ 2 ) = L(µ * ) = 1 {µ ′ 2 =µ ′ * 2 } P µ I ,σ 2 I (Θ I = µ * I ) = 0,
where the second equality follows by Assumption 1 and the last equality follows from the fact that Θ I follows a non-singular multivariate normal distribution, so it induces on (R dµ I , B(R dµ I )) a measure which is absolutely continuous w.r.t. the Lebesgue measure which assigns measure 0 to the singleton {µ * I }. We conclude that for any φ ∈ M × R dµ 1 + \ {0} dµ 1 , we have L(µ, σ 2 ) < L(µ * ) and consequently the MLE φ * is unique. ⇐) If an MLE exists on Φ, then L is upper bounded on Φ by L(φ * ), and consequently is upper bounded on M × {0} dµ 1 by the same value. By [START_REF] Moral | Non-linear filtering: interacting particle solution[END_REF] this implies that L is upper bounded on M by L(φ * ). If φ * ∈ M × {0} dµ 1 , then [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF] and [START_REF] Dempster | Maximum Likelihood from incomplete data via the EM algorithm[END_REF] are satisfied for some µ * ∈ M . Since the upper bound of L is attained on µ * the MLE exists. It is also unique, since if µ * 1 and µ * 2 are two distinct MLE, then (µ * 1 , 0 dµ 1 ) and (µ * 2 , 0 dµ 1 ) are two distinct MLE associated with m µ 1 , and this contradicts the initial assumption. Now we will show that

φ * ∈ M × {0} dµ 1 is a necessary condition. Let us assume that φ * ∈ M × R dµ 1 + \ {0} dµ 1 . Take the nonempty set I = {i ∈ {1, . . . , d µ 1 } : (σ 2 i ) * > 0}
and J the complement of I. By applying [START_REF] Doucet | Sequential Monte Carlo methods in practice[END_REF] for φ * we have

L(φ * ) = E µ * I ,(σ 2 I ) * L(Θ I , (µ ′ 2 ) * ) . ( 13
) If L(µ I , (µ ′ 2 ) * ) < L(φ *
) for all µ I ∈ R dµ I then this would contradict [START_REF] Gaucherel | Parameterization of a process-based tree growth model: Comparison of optimization, mcmc and particle filtering algorithms[END_REF]. So, we conclude that there exists µ * * I such that L(µ * * I , (µ

′ 2 ) * ) = L(φ * ), or φ * * = (µ * * I , (µ ′ 2
) * , 0 dµ 1 ) = φ * is an MLE, which contradicts the assumption of uniqueness. Therefore, φ * ∈ M × {0} dµ 1 and the proof is complete.

A consequence of Theorem 1 is that we can transfer the maximization problem corresponding to the initial model m to a maximization problem corresponding to the model m µ 1 . But since the latter is formulated as an incomplete data problem, we could design an appropriate version of the EM-algorithm to solve the maximization problem or more often a stochastic variant of the EM-algorithm (Expectation-Maximization, see [START_REF] Dempster | Maximum Likelihood from incomplete data via the EM algorithm[END_REF]) to handle the usually non-explicit state estimation problem. It is noteworthy that if all parameters belonging to µ are randomized then the M-step becomes explicit. Otherwise, a GEM (G:Generalized) or stochastic GEM could be used to find the MLE. An example of this kind could be an ECM (C:Conditional), see [START_REF] Meng | Maximum likelihood estimation via the ecm algorithm: A general framework[END_REF], where parameters are updated in a cyclic fashion. In the following Proposition we give the form of the Q-function, and the update formulas for the randomized parameter µ 1 .

Proposition 2. Let φ ′ = (µ 1 , µ 2 , σ 2 ) ′ ∈ M × (R +
* ) dµ 1 be the current parameter update associated with the model m µ 1 . The update equations for the parameter (µ 1 , σ 2 ) are given by:

µ 1 = E φ ′ Θ | Y = y (14)
σ 2 = ( σ 2 i ) 1≤i≤dµ 1 = Var φ ′ (Θ i | Y = y) 1≤i≤dµ 1 ( 15
)
and µ 2 can be updated independently by maximizing

Q 2 (µ 2 ; φ ′ ) = E φ ′ log p(y|Θ; µ 2 ) | Y = y . ( 16
)
Proof: Since φ ′ ∈ M × (R * + ) dµ 1 , this parameter corresponds to the Gaussian randomization m(µ 1 ) and by the definition of the Q-function and Definition 1 we have

Q(φ; φ ′ ) = E φ ′ log p(Θ, y; µ 1 , µ 2 , σ 2 ) | Y = y (17) = E φ ′ log p(Θ; µ 1 , σ 2 ) | Y = y + E φ ′ log p(y|Θ; µ 2 ) | Y = y . ( 18
)
The last term in the right-hand member of equation ( 18) depends only on µ 2 and coincides with ( 16), while the first term denoted by Q 1 (µ 1 , σ 2 ; φ ′ ) does not depend on µ 2 .

Consequently, the last statement of this theorem is true. By using the independence of the components of Θ and the density of the normal distribution we have that Q 1 can be maximized equivalently as follows:

Q 1 (µ 1 , σ 2 ; φ ′ ) max ∼ - 1 2 dµ 1 i=1 log σ 2 i - 1 2 dµ 1 i=1 E φ ′ {(Θ i -µ 1,i ) 2 | Y = y} σ 2 i . (19)
Consequently, it is elementary to show by relation ( 19) that ( 14) and ( 15) hold.

Remark 1.

A question which arises naturally concerns the adaptation of this method in the case that some of the components of the parameter µ have range different than R.

Then a Gaussian randomization as explained in Definition 1 cannot be applied directly and a modification is needed by suitable reparameterizations whenever possible. Examples of this kind are given in subsection 4.3.

A toy example of Gaussian randomization

In order to illustrate the proposed method, we start this section by a simple example where the MLE is known. In the next section we consider an application to a model which arises from plant growth modeling and the MLE has no explicit form.

Let us consider a simple Gaussian model. In particular, take a vector Y = y of n i.i.d.

observations {Y i = y i } n i=1 from a univariate normal distribution N (µ y , σ 2 y ), where µ y ∈ R is an unknown parameter, σ 2 y ∈ R * + .
We study the effect of Gaussian randomization when i) σ 2 y is known and ii) σ 2 y should be estimated. Let us first treat case i). Since the initial model m is Gaussian it is well-known that the MLE exists, is unique and is given by the sample mean, that is, µ * y = ȳ. Now assume that a Gaussian randomization is introduced for the parameter µ y . In particular, set µ = µ y and consider the Gaussian randomization m(µ). In this case,

Y ∼ N n Θ1 n , σ 2 y I n , Θ ∼ N (µ, σ 2 ), ( 20
)
where 1 n is the column-vector of ones and I n is the identity matrix of rank n. By writing,

Y i = µ + σξ + σ y ǫ i , i = 1, . . . , n,
where {ξ, ǫ i ; i = 1, . . . , n} are mutually independent standard Gaussian random variables, we get easily that

Y ∼ N n µ1 n , Σ σ 2 , ( 21
)
where the elements of the covariance matrix Σ σ 2 are given by

Σ σ 2 (i, j) = σ 2 + σ 2 y if i = j, σ 2 if i = j. (22)
Notice that in the extended model m(µ) the r.v. {Y i } n i=1 are no longer independent. Now, we take φ = (µ, σ 2 ), Φ = R × R * + and we write the equivalent log-likelihood in the sense of maximization

L(φ) max ∼ - 1 2 log det Σ σ 2 - 1 2 (y -µ1 n ) ⊤ Σ -1 σ 2 (y -µ1 n ), φ ∈ Φ. ( 23
)
By Theorem 1, we should have that a continuous extension of L(φ) from Φ to Φ = R × R + is feasible, the MLE exists in the extended model, is unique and is given by φ * = (µ * , 0).

As far as the extension is concerned, we verify by ( 23) that lim σ 2 →0 L(φ) = L(µ), where L(µ) is the likelihood function of the initial model m, and consequently a continuous extension is possible. In order to verify the aforementioned properties of the MLE, it is advantageous to use the simple form of the covariance matrix Σ σ 2 given by [START_REF] Musso | Improving regularized particle filters[END_REF]. We can

easily get that det Σ σ 2 = (σ 2 y ) n-1 (nσ 2 + σ 2 y ), n ≥ 1, (24)
and also that

Σ -1 σ 2 (i, j) =        1 σ 2 y 1 - σ 2 nσ 2 + σ 2 y if i = j, n ≥ 1, - σ 2 σ 2 y (nσ 2 + σ 2 y ) if i = j, n ≥ 2. ( 25
)
Now, by inserting ( 24) and ( 25) in ( 23), we get

L(φ) max ∼ - 1 2 log(nσ 2 + σ 2 y ) - 1 2σ 2 y n i=1 (y i -µ) 2 + σ 2 2σ 2 y (nσ 2 + σ 2 y ) n i=1 (y i -µ) 2 . (26)
By taking the partial derivative of (26) w.r.t. µ we get

∂ L(φ)/∂µ = 1 nσ 2 + σ 2 y n i=1 (y i -µ). ( 27
)
By setting [START_REF] Wei | A Monte Carlo implementation of the EM algorithm and the poor man's data augmentation algorithms[END_REF] equal to zero, we get indeed as unique solution µ * the sample mean ȳ. By inserting this solution into (26) the third term of the right-hand member of (26) gets null.

So, the maximization w.r.t. σ 2 is equivalent to the maximization of the first term, which is a strictly decreasing function of σ 2 and consequently attains its maximum at σ 2 = 0.

This shows that φ * = (µ * , 0) as expected. Conversely, and this is the most interesting direction in applications, since the MLE has the form (µ * , 0) in the extended model, we conclude by Theorem 1 that µ * is the MLE in the initial model.

It is noteworthy that in interesting problems that this idea can be applied, the direct maximization problem in the extended model will in general be no easier than the maximization in the initial model. But, its formulation as a hidden variable model gives directly a recursive procedure via an appropriate version of the EM algorithm to estimate the parameters of the model. In this example the M-step is explicit since the conditional distribution of Θ given Y is normal. Let φ ′ denote the current update of φ. We obtain easily that

[Θ|Y = y] ∼ N n(σ 2 ) ′ ȳ + σ 2 y µ ′ n(σ 2 ) ′ + σ 2 y , (σ 2 ) ′ σ 2 y n(σ 2 ) ′ + σ 2 y . (28)
Consequently, by Proposition 2, the update equations for µ and σ 2 are given by the mean and the variance of this normal distribution respectively. Note also that the likelihood of the extended model given by ( 23) has a unique stationary point which coincides with the MLE, so the sequence of EM iterates given by ( 28) is guaranteed to converge to the MLE.

Let us now treat case ii). The previous analysis shows that the MLE of (µ y , σ 2 ) does not involve the parameter σ 2 y and for this reason the MLE of µ y and σ 2 are identical with those of case i). After maximizing L(φ) given by [START_REF] Oudjane | Multiple model particle filter[END_REF] with respect to (µ y , σ 2 ), the MLE

An application to a plant growth model

Let us consider a nonlinear dynamic stochastic model of individual plant growth, the LNAS (Log-Normal Allocation and Senescence) model [START_REF] Cournède | Development and evaluation of plant growth models: Methodology and implementation in the PYGMALION platform[END_REF], whose aim is to describe the growth of the sugar beet plant in terms of the evolution of its biomass (mass for living organisms) with a daily time step. This model can therefore be regarded as a simplified version of the GreenLab model [START_REF] Lemaire | A morphogenetic crop model for sugar-beet (beta vulgaris l.)[END_REF] which is more adapted to experimental trials when only compartment data are available.

Description of the model

In this example, we present the case where the structure of the sugar beet can be simplified as the set of two types of organ compartments, the leaves and the root, from now on labeled as compartment 1 and 2 respectively. The development of sugar beet can be described via three basic mechanisms, biomass production, biomass allocation to these compartments and leaf senescence.

First, we describe the procedure of biomass production. The initial mass of the plant, corresponds to the mass of the seed and is denoted by q 0 . The subsequent biomasses q n , n ≥ 1, result from photosynthesis and under no modeling errors are assumed to be given by an adaptation of the empirical Beer-Lambert law [START_REF] Guo | Parameter optimization and field validation of the functional-structural model greenlab for maize[END_REF],

q n = µ a • PAR n • 1 -exp{-λ q g n } , (33)
where µ a is a strictly positive biological parameter related to the radiation use efficiency, q g n is the total mass of green leaves at day n (leaves that are still (photosynthetically active), PAR n is the photosynthetically active radiation at the n-th day and λ an empirical parameter. In [START_REF] Lemaire | A morphogenetic crop model for sugar-beet (beta vulgaris l.)[END_REF], λ is decomposed in the product of the Beer-Lambert extinction coefficient k by the specific leaf area divided by S p , a characteristic surface related to inter-individual competition.

In the sequel, for notational simplicity, we denote by F n the function given by (33) which relates q n with q g n by omitting the dependence on the other variables. In order to give some flexibility to this model and to explain some observed variability around the theoretical values given by (33), we consider a modified model of the form:

Q n = F n (Q g n ; µ a , λ) e ξ q n , n ≥ 1, ( 34
)
where Q n and Q g n are the stochastic counterparts of q n and q g n respectively, related via the function F n given by (33), and {ξ q n } n≥1 is an independent sequence with ξ q n ∼ N (0, σ 2 q ). Now, we describe the procedure of biomass allocation. Since we assume two organ compartments, the allocation of each produced biomass Q n is divided into two parts, one for the leaves compartment and one for the root. Let G(x; µ, σ 2 ) denote the distribution function of the log-normal distribution with location and scale parameters given by µ and σ respectively, and G c (x; µ, σ 2 ) its corresponding survival function, that is, G c (x) = 1 -G(x). Let also τ be the cumulative daily temperature (thermal time in ❽.days) and define the function

γ(τ ; γ 0 , γ f , µ γ , σ 2 γ ) = γ 0 + (γ f -γ 0 ) G(τ ; µ γ , σ 2 γ ), τ > 0, ( 35
)
where γ 0 , γ f are parameters such that 0 < γ 0 < γ f < 1. The allocated percentage to compartment 1 (leaves) on day n associated with thermal time τ n , denoted by Γ n , is assumed to be a logit-normal distribution with mean the logit transform of γ(τ n ; γ 0 , γ f , µ γ , σ 2 γ )) and variance σ 2 γγ , that is,

log Γ n 1 -Γ n = log γ(τ n ; µ γ , σ 2 γ ) 1 -γ(τ n ; µ γ , σ 2 γ ) + ξ γ n , (36)
where ξ γ n ∼ N (0, σ 2 γγ ). In this way, the total leaf mass and the total root mass can be computed recurrently by:

Q f n+1 = Q f n + Γ n Q n , (37)
Q r n+1 = Q r n + (1 -Γ n )Q n , (38)
where the initial conditions Q f 0 = Q r 0 = 0 hold, Q n and Γ n satisfy (34) and (36) respectively.

The last important process to consider is leaf senescence. We model with the following equation the senescent leaf mass Q s n until day n, which corresponds to the part of the total leaf mass at day n that is no longer photosynthetically active (it may have already fallen from the plant or is still on the plant but yellow). We assume that

Q s n = G(τ n ; µ s , σ 2 s ) Q f n . (39)
By (39) we deduce directly that the total green leaf mass is given by:

Q g n = Q f n -Q s n = G c (τ n ; µ s , σ 2 s ) Q f n . ( 40
)
The potentially available data consist in the compartment mass measurements {Q g n , Q r n } N n=1 , where N is a given time horizon, observed under measurement errors. Let us denote Z n,1 and Z n,2 the measurements of Q g n and Q r n respectively corrupted by error. In particular, we assume that for n = 1, . . . , N ,

Y n = (log Z n,j ) j = log Q g n + ǫ g n log Q r n + ǫ r n = log G c (τ n ; µ s , σ 2 s )Q f n + ǫ g n log Q r n + ǫ r n , (41)
where {ǫ g n , ǫ r n } n≥1 are mutually independent and ǫ g n ∼ N (0, σ 2 g ), ǫ r n ∼ N (0, σ 2 r ).

We are now in the position to represent the model given by equations ( 34)-(41) as a state space model.

State space model representation

A possible choice of the state vector of the model is

X n = (Q f n , Q r n ).
In order to have the state and the observation equation of the model we need the following representation:

X n+1 = Φ n (X n , ξ n ; µ), (42)
Y n = R n (X n , ǫ n ; µ), (43)
where

X n = (X n,1 , X n,2 ) = (Q f n , Q r n
) and {ξ n } n≥1 and {ǫ n } n≥1 are appropriately chosen sequences of random vectors. The parameter vector of our model is denoted by µ and includes all structural parameters {µ a , λ, µ γ , σ 2 γ , γ 0 , γ f } eventually together with some of error variances {σ 2 γγ , σ 2 q , σ 2 g , σ 2 r }. In the following lemma we explicit the state equation and the observation equation corresponding to our model. 34)-(41) can be described as a state space model. The state equation is given by the following system of equations:

Lemma 3. The model given by equations (

X n+1,1 = Φ n,1 (X n , ξ n ; µ) = X n,1 + F n G c (τ n ; µ s , σ 2 s ) X n,1 ; µ a , λ e ξ n,1 1 + γ -1 (τ n ; µ γ , σ 2 γ ) -1 e -ξ n,2 (44)
X n+1,2 = Φ n,2 (X n , ξ n ; µ) = X n,2 + F n G c (τ n ; µ s , σ 2 s ) X n,1 ; µ a , λ e ξ n,1 1 + γ -1 (τ n ; µ γ , σ 2 γ ) -1 -1 e ξ n,2 , ( 45
)
where {ξ n } n≥1 is an independent from X 0 sequence of mutually independent bi-dimensional centered Gaussian random vectors with diagonal covariance matrix and variances σ 2 1 and σ 2 2 for X n,1 and X n,2 respectively. The observation equation is given by

Y n = R n (X n , ǫ n ) = log G c (τ n ; µ s , σ 2 s ) + log X n,1 log X n,2 + ǫ n,1 ǫ n,2 , ( 46
)
where {ǫ n } n≥1 is an i.i.d. sequence of 2-dimensional centered Gaussian random vectors with diagonal covariance matrix, variances σ 2 i for each ǫ n,i , i = 1, 2, and assumed to be independent from {ξ n } n≥1 .

Proof: Due to the choice of the state vector X n = (Q f n , Q r n ) and equations ( 37) and (38), we only have to show that Γ n Q n and (1 -Γ n )Q n coincide with the second summand of the right-hand member of ( 44) and (45) respectively. In both terms the numerator is common and equals Q n . Indeed, this can be deduced by (34) by substituting Q g n with its expression in (40). Let us now denote by φ : (0, 1) × R → R such that φ(x, y) = log(x/1x) + y. Note by (36) that the left-hand member of this equation corresponds to the logit transformation of Γ n , which is the inverse of the standard logistic sigmoid function. We deduce that

Γ n = 1 1 + e -φ(γ(τn),ξ n,2 ) = 1 1 + γ -1 (τ n ) -1 e -ξ n,2 , (47)
and

1 -Γ n = 1 1 + e φ(γ(τn),ξ n,2 ) = 1 1 + γ -1 (τ n ) -1 -1 e ξ n,2 . ( 48
)
By multiplying both members of (47) and (48) with Q n given by (34) and substituting these expressions in the state equations we obtain (44) and (45). The observation equation (46) of the model can be deduced directly by (41).

Gaussian randomization of the plant growth model

In this part we apply the ideas developed in Section 2 for the state space model described in Lemma 3. Assume for simplicity of presentation that all error variances {σ 2 γγ , σ 2 q , σ 2 g , σ 2 r } are known or a priori fixed to some values. The interest will be focused on the estimation of the structural parameters and let µ = (µ a , λ, µ γ , σ γ , γ 0 , γ f ). Note also that this state space model can equivalently be described as a non-homogeneous hidden Markov model [START_REF] Cappé | Inference in Hidden Markov Models[END_REF], and usual estimation techniques are based on some variant of the EM algorithm, by combining an approximation of the E-step and a usually explicit M-step. Nevertheless, we can easily see that even in the complete log-likelihood of this model there is no explicit solution for µ, therefore leading to a non-explicit M-step. Otherwise, conditional maximization approaches should be invoked and this often leads to complicated generalized EM-type algorithms (see, eg., [START_REF] Trevezas | A sequential Monte Carlo approach for MLE in a plant growth model[END_REF]), where a numerical maximization procedure should be implemented as well. It is also well known (see, eg., [START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF]) that the EM-type sequences do not necessarily converge to the true MLE (if it exists). In typical plant growth models with a large number of parameters, it is very difficult to ensure the existence of the MLE and the subsequent convergence of the EM iterates to the true MLE. The best that we can hope is to have a robust algorithm and the solutions to the initial maximization problem to be biologically relevant. For this reason, we take Assumption 1 for granted.

Let us first consider a suitable reparameterization of the state space model m in such a way that a Gaussian randomization can be applied. Let µ = (µ a , λ, µ γ , σ γ , γ 0 , γ f ). By Remark 1 and since µ a , σ 2 γ ∈ (0, +∞) we deduce that a possible reparameterization can be given by η = (log µ a , λ, µ γ , log σ 2 γ , γ 0 , γ f ), where the range of η is R 6 . The functional parameters µ s and σ s which are related to the senescent process can be estimated only based on the senescent leaf mass and thus are fixed in our study case. The results of Section The parameter of the extended model is given by φ = (η, σ 2 ) ∈ R 6 ×(R * + ) 6 . By Proposition 2, the problem of parameter estimation of φ is transformed in an iterative state estimation problem. Starting with a fixed value φ (0) , the EM-update equations given by ( 14) and [START_REF] Le Gland | An analysis of regularized interacting particle methods for nonlinear filtering[END_REF] produce a sequence φ (n) which converges to a stationary point of the initial likelihood.

Unfortunately, the conditional moments involved in these equations cannot be computed explicitly for nonlinear models of this type and for this reason we need to implement a stochastic variant of the EM-algorithm to approximate them. In order to tackle this problem, we describe in the next section an appropriate particle filtering method.

Particle filtering approximation of the E-step

At each iteration of the EM-algorithm, given the current parameter update φ ′ , our objective is to approximate

E φ ′ Θ | Y 0:n = y 0:n and Var φ ′ (Θ i | Y 0:n = y 0:n ) 1≤i≤dµ 1
where n corresponds to the observation length, since the update equations for the parameters of the model given in Proposition 2 generally lead to non-explicit solutions. This is also the case for the specific application that we present in the previous section. One could think of several alternative algorithms to make this approximation, like Markov chain Monte-Carlo or particle filters. In this paper we opt for the latter type of algorithms.

The Monte-Carlo particle filters (PFs) have been studied a lot in the literature. The most adapted to our application context can be cast into the framework of the generic sequential importance sampling (SIS) algorithm, or SISR, when a resampling step is introduced (see, e.g., [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-gaussian bayesian tracking[END_REF] or [START_REF] Doucet | Sequential Monte Carlo methods in practice[END_REF]). The first efficient implementation of Monte-Carlo particle filters of this type dates back to Gordon, Salmond and Smith (1993). The authors give a generic formulation of the bootstrap PF, using importance sampling ideas in the correction step and systematic use of resampling in order to avoid degeneracy of the particle filter. For this reason, the algorithm is also named SIR algorithm (sampling/importance resampling). The term interacting particle filter (IPF) is also used by [START_REF] Moral | Non-linear filtering: interacting particle solution[END_REF] for the same PF algorithm. For other names of the same algorithm see [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-gaussian bayesian tracking[END_REF] (page 178) and the references therein. Before describing an appropriate particle filter for our purposes let us explain how the filtering approach is feasible in our context.

The filtering approach

In our application context the initial model m is a state space model and Θ can be considered as a part of the hidden state vector. The most natural way of dealing with Θ dynamically in time is to construct copies Θ n which remain constant in time, that is, for all n, we have Θ n = Θ n+1 . Despite these constant dynamics, the version of the particle filter that we describe in the sequel enables us to modify sequentially in time the distribution of Θ given the increasingly data vector. The corresponding augmented state vector will be denoted by X a n = (X n , Θ n ). In this context, all we need to approximate is

E φ ′ Θ n | Y 0:n = y 0:n and Var φ ′ (Θ n,i | Y 0:n = y 0:n ) 1≤i≤dµ 1
. If we have a way to approximate the intractable filtering density p(x a n |y 0:n ) by a suitable p(x a n |y 0:n ), then since x a n = (x n , θ n ), we can obtain an approximation of the first two conditional moments of Θ n given the complete data vector with the help of the marginal p(θ n |y 0:n ). In this direction, we describe in the sequel an appropriate particle filtering technique to make this approximation.

The regularized particle filter

A particular feature of the Gaussian randomization is that the optimal variance parameters equal zero (see Theorem 1). When the variances are close to zero, and this will indeed happen after some iterations of the EM-algorithm, classical particle filters cannot perform well. This is due to the fact that the particles make a discrete approximation of the filtering distribution and consequently degenerate easier under a low variance scenario due to the resampling mechanism (this is known as sample impoverishment). A successful strategy to overcome this problem was devised in [START_REF] Musso | Regularization schemes for branching particle systems as a numerical solving method of the nonlinear filtering problem[END_REF] with the so-called postregularized particle filter (post-RPF). The regularization refers to the use of a kernel smoothing method to change the discrete approximation of the filtering density (induced by the weighted particles in the classical IPF) into an absolutely continuous approximation. The term post-RPF is used to indicate that the regularization step takes place after the correction step. A comprehensive exposition of RPFs and improvements can be found in [START_REF] Musso | Improving regularized particle filters[END_REF] and some theoretical results in [START_REF] Le Gland | An analysis of regularized interacting particle methods for nonlinear filtering[END_REF] and [START_REF] Gland | Stability and uniform approximation of nonlinear filters using the Hilbert metric and application to particle filters[END_REF]. In [START_REF] Musso | Regularization schemes for branching particle systems as a numerical solving method of the nonlinear filtering problem[END_REF] and [START_REF] Oudjane | Multiple model particle filter[END_REF], the authors compare the performance of the RPFs with the IPFs in some classical tracking problems. Let us now describe the adaptation of this method in our context. We denote by K the regularization kernel (usually Epanechnikov or Gaussian) associated to X a n and by K h the rescaled kernel given by:

K h (x a ) = h -d K(h -1 x a ),
where h > 0 is the bandwidth parameter and d is the dimension of the state vector. In our application, we used the Gaussian kernel and with M particles the optimal bandwidth parameter (when the underlying density is Gaussian) is given by (see, eg., [START_REF] Musso | Improving regularized particle filters[END_REF]):

h M = 4 d + 2 1 d+4 M -1 d+4 . (49)
Let us assume that at the end of the n-th step, the estimate of p(x a n |y 0:n ) is given by:

p(x a n |y 0:n ) = M i=1 w(i) n K h M (x a n -xa (i) n ),
where {(

w(i) n , xa (i) n )}
is the weighted sample before the regularization step. The basic filtering step can be described as follows:

❼ Sampling with regularization: For all i = 1, . . . , M :

-Generate I (i) ∈ {1, . . . , M }, with P(I (i) = j) = w(j) n . -Generate ǫ (i) ∼ K(x a ). -Compute xa (i) n = xa (I (i) ) n - + h M Σ1/2 n ǫ (i)
, where Σ1/2 n is the square root of the empirical covariance matrix (whitening is used)

❼ Prediction step -Generate x(i) n+1 ∼ p(x n+1 |x a (i) n ). -Set θ(i) n+1 = θ(i) n .
❼ Correction step and regularization:

-Set w(i) n+1 = p(y n+1 |x a (i) n+1 ), for all i = 1, . . . , M , -Regularize the weighted sample {( w(i) n+1 , xa (i) n+1 )} by taking p(x a n+1 |y 0:n+1 ) = M i=1 w(i) n+1 K h M (x a n+1 - xa (i) n+1 ).
As far as the regularization is concerned, it is noteworthy to mention the recent approach presented in [START_REF] Rossi | Nonlinear filtering in discrete time: A particle convolution approach[END_REF] (see also [START_REF] Campillo | Convolution Particle Filter for Parameter Estimation in General State-Space Models[END_REF]). The authors proposed a particle filter where a convolution kernel is used to regularize the likelihood of the observations as well. The new filter, called convolution particle filter, seems to be more appropriate to use when the likelihood cannot be computed explicitly or when signal-to-noise ratio is very low or very high (see, eg., [START_REF] Rossi | Nonlinear filtering in discrete time: A particle convolution approach[END_REF]). 1: Comparison of the estimations from two different initializations, both for the restricted (150000 particles) and the complete dataset (40000 particles). ♠ : estimated loglikelihood based on the restricted dataset (14 dates); ♣ : estimated log-likelihood based on the complete dataset (160 dates). of the restricted dataset, there is a difference in the estimates of the last three parameters (γ f , µ γ , σ γ ) and the effect of the initialization seems to be non negligible, especially in the case of σ γ . This can be explained by the sensitivity of model outputs Q g and Q r to parameters. A sensitivity analysis of the LNAS model based on Sobol's method was presented in [START_REF] Cournède | Development and evaluation of plant growth models: Methodology and implementation in the PYGMALION platform[END_REF] and showed that the last three parameters are the least influential ones among the six which were estimated. The increased variability of these parameters (see the corresponding standard deviations) from independent runs of the algorithm (20 repetitions) based on the same dataset is also an argument in this direction, since the great loss of information seems to affect their estimation quality. The results of loglikelihood evaluations also seem to suggest that despite that the estimated MLE could slightly differ from two different initializations, the most sensitive functional parameters can be estimated without difficulty even in the case of the restricted dataset, together with the observation and model noise parameters. By comparing the estimation results (with the same initialization) between the restricted and the complete dataset, we also conclude that the estimates of the most influential parameters are closer to the real values with the complete dataset than the restricted one, as expected, due to the great amount of missing information.

In Fig. 1, we illustrate the convergence of the estimations of the six functional parameters (for the restricted dataset) based on a single run of the algorithm. The uncertainty related to parameter estimation is evaluated by parametric bootstrap. As a better approximation of the MLE, we considered the means from the estimates of both initializations.

In Table 2 we present the results that we obtained with a bootstrap sample of size 100.

Except for the least influential parameters, the bootstrap means are very close to the MLE in the case of the complete dataset, but when a great amount of data is missing a small bias is introduced. By comparing the standard deviations of the estimates for the complete and the restricted dataset, it is clear that the additional amount of information reduces considerably the standard deviations of all the estimated parameters, except for the 3 least influential functional parameters. Moreover, the standard deviations of γ f and σ γ are prohibitively high to allow for reliable estimations. However, for all the other parameters, despite the increased standard deviations in the case of the restricted dataset, we finally get reliable estimates of the true values within the range of one standard In order to alleviate the burden of a large constant number of particles, we implemented and compared two different strategies of increasing the sample size during iterations, that is, a geometric and a quadratic increase. The geometric increase provided better results.

In Table 3, we compare the results of the geometric increase with those that we obtained with the constant number of particles. The increase of the particles' number takes place every time that the Monte Carlo error in the evaluation of the Q-function is considered to be significant, for example, by detecting a zig-zagging in parameter estimates. Without a strong effort to optimize the increasing strategy, the computational time was reduced by 23%. The results show that the geometric increase induces slightly higher variability in the estimates among different runs of the algorithm than in the case with a constant particle number. Consequently, the gain in computational time from the decrease of the total number of simulations is counterbalanced by this effect. In this paper, we do not focus on aspects concerning the most efficient use of Monte Carlo resources and further

A real-data application

In this section, we present the identification results from a real dataset with the LNAS model (described in Section 4), where parameter estimation is performed with the RPF-EM algorithm, as described in Section 5.

This dataset concerns limited experimental observations realized by the French institute for Sugar Beet research in 2010 (details of the experimental protocols are presented in [START_REF] Lemaire | Analysis of the density effects on the source-sink dynamics in sugar-beet growth[END_REF]). Dry matter of root (denoted by Q r ) and leaves (denoted by Q g ) were collected

on 50 plants at 14 different dates. The observation vector, which records at each available date the average mass per square meter, is given by Table 6. The details of the model are given in Section 4.3. The estimation results are presented in Table 5 together with estimated standard errors from 20 independent runs. The Monte Carlo sample size was increased geometrically and other implementation details can be found in Section 5.3.

Note that the influence of the initialization is negligible for all the parameters but σ γ , which is the least influential parameter. This is also suggested by the small difference of the log-likelihood evaluations at the corresponding solutions. The variability from independent runs of the algorithm was relatively high for the modeling noise parameters (σ q , σ γγ ), but by averaging a moderate number of independent runs a better precision can be obtained.

The uncertainty related to parameter estimation was assessed by parametric bootstrap, by averaging the estimates obtained from both initializations. In Table 6, we present the means and the standard deviations of a 100-bootstrap sample. Notice that the estimates of the 3 most important parameters (high sensitivity indices) are very satisfactory, with low standard deviations relatively to their mean values. Nevertheless, the estimations from the least influential parameters are less satisfactory with higher standard deviations as expected from the synthetic case in Section 5.3. A small bias in the mean estimated values from the bootstrap sample also confirms the difficulty in estimating these parameters.

Note also that the estimates of the noise parameters are satisfactory.

Conclusion

In the context of complex dynamic stochastic models, parameter estimation is a difficult and challenging problem. For a likelihood-based approach, many complex stochastic EM-type algorithms have been proposed to approach the maximum likelihood estimator.

In real data applications the M-step is never explicit and complex numerical algorithms should be combined to make parameter estimation feasible. The proposed Gaussian randomization process allows us to make the M-step explicit and to circumvent the numerical difficulties, under the assumption that the initial model has a unique maximum likelihood estimator. The proposed RPF-EM algorithm, which consists of an E-step performed by the post-regularized particle filter circumvents classical problems such as the divergence of the filters induced by weak noises and low variance scenarios during the last iterations of the algorithm. Additionally, noise related variance parameters can easily be estimated by conditional maximization steps. Synthetic examples and results based on a real dataset suggest that the proposed algorithm performs well when confronted to scarce observations. As a conclusion, the proposed RPF-EM algorithm based on Gaussian randomization appears as a promising alternative for robust estimation of plant growth models, and moreover in a general way for models in life sciences.

1 . 1 .

 11 Definition Let m be a statistical model which satisfies Assumption 1 and µ = (µ 1 , µ 2 ), where µ 1 ∈ R dµ 1 . The statistical model m(µ 1 ) will be called Gaussian randomization of m w.r.t. µ 1 if m(µ 1 ) is an incomplete data model, which consists in: i) a Gaussian hidden vector Θ, where Θ ∼ N dµ 1 (µ 1 , Σ), and Σ = diag{σ 2 i } 1≤i≤dµ 1 , where σ 2 i > 0, ii) an observed vector Y , where the conditional distribution of Y given Θ = θ depends only on the parameter µ 2 and satisfies:

Lemma 1 .

 1 The likelihood function L(φ) of the extended model m(µ 1 ) as given in Definition 1 is upper bounded by L(µ * ).

Theorem 1 .

 1 Let us consider the class m µ 1 of all Gaussian randomizations { m(µ I )} I of the initial model m, where I ⊂ {1, . . . , d µ 1 } (for I = ∅, m(µ ∅ ) = m) and µ I is defined in Proposition 1. If Assumption 1 holds, then the unique MLE µ * associated with the model m determines a unique MLE φ * associated with the model m µ 1 . In particular, we have

2

  are now directly applicable to η. In particular, let m(η) the Gaussian randomization of m w.r.t. η. By Definition 1 we have that m(η) is an incomplete data model which consists in: i) a Gaussian hidden vector Θ, where Θ ∼ N 6 (η, Σ), and Σ = diag{σ 2 i } 1≤i≤6 , where σ 2 i > 0, ii) an observed vector Y , where conditioned on Θ = θ is a state space model which satisfies (42) and (43), with the parameter µ replaced by θ.

Figure 1 :

 1 Figure 1: Evolution of the estimated values of six functional parameters: µ a , λ, γ 0 , γ f , µ γ and σ γ with the RPF-EM algorithm (150000 particles) for the restricted dataset.

Table

				Initialization 1
			Initial distribution	Restricted dataset ♠	Full dataset ♣
	Parameter	Real value	Mean	Std. Mean estimate Std. Mean estimate Std.
	µ a	3.56	3.40	0.15	3.664	0.049	3.564	0.005
	λ	56.6	60.0	5.0	59.73	0.22	55.68	0.12
	γ 0	0.625	0.58	0.10	0.679	0.026	0.628	0.002
	γ f	0.1035	0.12	0.025	0.088	0.013	0.047	0.006
	µ γ	550.0	500.0	50.00	492.91	2.44	624.55	9.07
	σ γ	950.0	880.0	50.00	1338.75	20.25	1633.29	98.66
	σ g	0.1	0.08	-	0.106	0.001	0.108	0.001
	σ r	0.1	0.08	-	0.083	0.001	0.091	0.001
	σ q	0.05	0.035	-	0.047	0.010	0.037	0.005
	σ γγ	0.05	0.035	-	0.059	0.006	0.061	0.005
	Log-likelihood	-136.875 ♠ -426.870 ♠	-	-135.514	0.181	-929.404	1.626
		-935.350 ♣ -17312.1 ♣
				Initialization 2
			Initial distribution	Restricted dataset ♠	Full dataset ♣
	Parameter	Real value	Mean	Std. Mean estimate Std. Mean estimate Std.
	µ a	3.56	3.60	0.15	3.668	0.058	3.566	0.003
	λ	56.6	52.0	5.0	59.94	0.50	55.73	0.07
	γ 0	0.625	0.7	0.10	0.686	0.036	0.628	0.002
	γ f	0.1035	0.09	0.025	0.070	0.024	0.048	0.005
	µ γ	550.0	580.0	50.00	510.24	4.60	623.14	6.69
	σ γ	950.0	1000.0	50.00	1610.93	46.88	1609.01	82.70
	σ g	0.1	0.08	-	0.106	0.002	0.108	0.001
	σ r	0.1	0.08	-	0.083	0.001	0.091	0.001
	σ q	0.05	0.035	-	0.048	0.012	0.040	0.004
	σ γγ	0.05	0.035	-	0.055	0.011	0.060	0.004
	Log-likelihood	-136.875 ♠ -388.465 ♠	-	-135.534	0.327	-929.261	1.264
		-935.350 ♣ -4106.015 ♣

Table 2 :

 2 Uncertainty assessment in the case of both restricted and complete simulated datasets by parametric bootstrap. The means and the standard deviations are given on the basis of a 100-bootstrap sample for the restricted (14 dates) and the complete (160 dates) datasets. investigations are needed in this direction.

				Restricted dataset ♠	Full dataset ♣
	Parameter Real value	MLE	Mean	Std.	MLE	Mean	Std.
	µ a	3.560		3.666	3.615	0.137	3.565	3.567	0.058
	λ	56.6		59.8	57.7		4.6	55.7	55.7	1.2
	γ 0	0.625		0.680	0.660	0.057	0.628	0.626	0.017
	γ f	0.1035		0.079	0.097	0.033	0.048	0.040	0.030
	µ γ	550.00		501.58 503.60 77.95 623.85 633.29 53.62
	σ γ	950.00		1474.84 1292.37 314.81 1621.15 1852.70 308.58
	σ g	0.1		0.106	0.101	0.016	0.108	0.106	0.005
	σ r	0.1		0.083	0.079	0.021	0.091	0.092	0.006
	σ q	0.05		0.047	0.049	0.012	0.038	0.038	0.004
	σ γγ	0.05		0.057	0.054	0.016	0.061	0.052	0.006
		Initial distribution		Constant →	Geometric increase ↑
	Parameter	Real value	Mean	Std.	Mean estimate Std. Mean estimate Std.
	µ a	3.56	3.40	0.15		3.664	0.049	3.662	0.082
	λ	56.6	60.0	5.0		59.73	0.22	59.67	0.54
	γ 0	0.625	0.58	0.10		0.679	0.026	0.677	0.063
	γ f	0.1035	0.12	0.025		0.088	0.001	0.089	0.001
	µ γ	550.0	500.0	50.00		492.91	2.44	495.27	7.71
	σ γ	950.0	880.0	50.00		1338.75	20.25	1372.40	32.92
	σ g	0.1	0.08	-		0.106	0.001	0.107	0.002
	σ r	0.1	0.08	-		0.083	0.001	0.085	0.001
	σ q	0.05	0.035	-		0.047	0.010	0.051	0.013
	σ γγ	0.05	0.035	-		0.059	0.006	0.061	0.009
	Log-likelihood -136.875 -426.870	-		-135.514	0.181	-135.648	0.173

Table 3 :

 3 Comparison of parameter estimates with two different configurations of particle numbers.

→: a constant number of 150000 particles. ↑: starting from 80000, the number of particles increases geometrically until a maximum number of 150000.

Table 4 :

 4 Experimental dataset provided by the French institute for Sugar beet research (ITB) based on an experiment in 2010. 14 dates of measurement are available (denoted by n). The observation vector contains the averaged mass (from 50 plants) per square meter of green leaf compartment (denoted by Y n,1 ) and root compartment (denoted by Y n,2 ) in g.

Table 5 :

 5 Comparison of the estimations resulting from two different initializations in the case of the 2010 experimental dataset.

			Initialization I			Initialization II
	Parameter	Initial distribution	Estimation	Initial distribution	Estimation
		Mean	Std.	Mean	Std.	Mean	Std.	Mean	Std.
	µ a	3.40	0.2	3.565	0.062	3.45	0.4	3.563	0.059
	λ	58.0	3.0	59.78	0.37	57.0	5.0	59.33	0.46
	γ 0	0.90	0.10	0.844	0.017	0.78	0.15	0.836	0.021
	γ f	0.18	0.05	0.190	0.002	0.20	0.05	0.198	0.001
	µ γ	650.0	40.0	642.64	5.96	630.00	50.0	641.02	6.46
	σ γ	350.0	30.00	319.59 13.62	300.0	40.0	297.78 14.80
	σ g	0.10	-	0.142	0.004	0.12	-	0.143	0.001
	σ r	0.10	-	0.164	0.004	0.12	-	0.167	0.001
	σ q	0.02	-	0.042	0.012	0.03	-	0.043	0.012
	σ γγ	0.05	-	0.064	0.012	0.03	-	0.064	0.012
	Log-likelihood -173.914	-	-165.586 0.116 -169.355	-	-165.591 0.132
			Parameter MLE Mean	Std.
			µ a	3.564 3.562 0.123
			λ	59.55 59.70	3.13
			γ 0	0.840 0.843 0.044
			γ f	0.194 0.182 0.053
			µ γ	642.33 654.89 62.40
			σ γ	308.69 336.62 109.95
			σ g	0.142 0.143 0.032
			σ r	0.165 0.158 0.024
			σ q	0.042 0.039 0.008
			σ γγ	0.064 0.073 0.015

Table 6 :

 6 Uncertainty assessment of estimates based on the 2010 experimental dataset performed by parametric bootstrap.

of σ 2 y can be obtained by maximizing the function

with respect to σ 2 y . This gives that (σ 2 y ) * = n -1 n i=1 (y iȳ) 2 , that is, the sample variance, which obviously coincides with the MLE of the initial model of i.i.d. normal observations. The point here is that these parameters that conditional on the others can be solved explicitly (as function of the others) in the initial likelihood formulation, they still keep this advantage when estimated with the EM algorithm. Indeed, in this example the second term of the Q-function given by ( 16) can be written as:

By maximizing the above quantity we get:

where µ y and σ 2 are the current updates of µ y and σ 2 given by [START_REF] Wyckoff | Predicting tree mortality from diameter growth: a comparison of maximum likelihood and bayesian approaches[END_REF], and since µ y → µ * and σ 2 → 0, we conclude by (30) that σ 2 y converges to (σ 2 y ) * .

Remark 2. It is noteworthy that the addition of these artificial parameters does not affect the MLE of the initial model, even if one could think at first sight that the extended model is not identifiable. For example, consider the extreme and trivial case where n = 1 and σ 2 y is known. We have one available observation and two parameters µ y and σ 2 to estimate.

But in fact here, the subsequent degeneration of σ 2 to zero, which we know that it is the true value, causes no problem. On the contrary, notice the difference if one wished to estimate in the initial model µ y and σ 2 y jointly (a truly non identifiable model). We have

We deduce that the likelihood function is unbounded and the MLE does not exist. Additionally, it does not admit any continuous extension to be able to approach from the larger model the obvious solution µ y = y when µ y is the only parameter to be estimated.

A simulation case study for a plant growth model

With the purpose of testing the performance of the proposed algorithm, we performed simulations based on the true experimental conditions of one experiment conducted by the French institute for Sugar beet research in 2010 (see [START_REF] Baey | Evaluation of the predictive capacity of five plant growth models for sugar beet[END_REF] for a detailed description of the experimental protocol). Here, we present the results of simulations in two different scenarios: (i) the complete dataset is observed, that is, 320 observations (from day 1 to day 160) and (ii) we have only access to 28 observations from this dataset, green leaf and root masses measured at the 14 dates included in O 2010 : The details of the model are given in Section 4.3. Different strategies were compared for the augmentation of the Monte Carlo sample size (geometric and quadratic) and a large constant sample size was also used to have a more ideal (but more expensive) EM algorithm (150000 particles for the restricted dataset and 40000 for the complete dataset).

Two different initializations were carried out in each case and the averaging technique ( [START_REF] Cappé | Inference in Hidden Markov Models[END_REF])

was used to smooth parameter estimates after a small burn-in period. As a stopping rule we used the one proposed in [START_REF] Booth | Maximizing generalized linear mixed model likelihoods with an automated monte carlo em algorithm[END_REF], which claims convergence when the relative change in the estimates from three successive iterations is reasonably small [START_REF] Booth | Maximizing generalized linear mixed model likelihoods with an automated monte carlo em algorithm[END_REF]. A conditional version of RPF-EM is implemented with the aim of providing reasonable estimates for the noise parameters (8 conditional steps for the estimations of the noise parameters).

The parameters' real values and the estimation results are presented in Table 1. The standard errors of the variability of parameter estimates from independent runs of the algorithm (20 repetitions based on the same dataset) are also given, which corresponds to the algorithmic uncertainty. The log-likelihood values of the estimates are estimated by the means of 10 independent evaluations (the standard errors were small, approximately 0.005).

In Table 1, we give an idea of the influence of the initial distribution. Note that for the complete dataset the estimates based on these two different initializations are very close, which suggests a negligible effect of the initial distributions. However, in the case


Y T Chen

M Rabilloud

J M Thiolet

B Coignard

M H Metzger

Benchmarking French regions according to their prevalence of healthcare-associated infections

Keywords: Hospital infection, Benchmarking, Cluster analysis, Hierarchical logistic regression modelling, Bayesian method

Objective: To propose an original method of benchmarking regions based on their prevalence of healthcare-associated infections and to identify regions with unusual results.

Design:

To study between-region variability with a three-level hierarchical logistic regression model and a Bayesian non-parametric method.

Setting: French 2006 national healthcare-associated infections point prevalence survey.

Participants: 336,858 patients from 2,289 healthcare facilities in 27 regions. Patients with an imported healthcare-associated infection (1% of the data, 20.7% of infected patients), facilities with less than 5 patients and patients who had at least one missing value for the variables taken into account were excluded (5.0% of patients).

Mean Outcome Measure: Binary outcome variable indicates whether a given patient was infected.

Results: Two clusters of regions were identified; one cluster of five regions had lower adjusted prevalence than the other one of 22 regions, while no region with unusually high prevalence could be identified. Nevertheless, the degree of heterogeneity of odds ratios between facilities for facility-specific effects of use of invasive devices was more important in some regions than in others.

Conclusions:

The adjusted regional prevalence of healthcare-associated infections can serve as an adequate benchmark to identify regions with concerning results. Although no outlier regions were identified, the proposed approach could be applied to the data of the 2012 national survey to benchmark regional healthcare policies. The estimation of facility-specific effects of use of invasive devices may orient future regional action plans.

INTRODUCTION

In France, control, prevention, and surveillance of healthcare-associated infections (HAIs) are organized geographically according to local, regional, and national administrative divisions.

Since the 1990s, all public and private healthcare facilities have been urged to constitute infection control committees and carry out IC programs. Currently, at the regional and interregional level, prevention and surveillance programs of health care facilities are developed with the support of a Centre de Coordination de la Lutte contre les infections nosocomiales (CClin) and at the national level by a steering committee coordinated by the French Ministry of Health. In addition, national HAI surveillance activities are coordinated by the Réseau d'Alerte, d'Investigation et de Surveillance des Infections Nosocomiales (RAISIN) [START_REF] The | RAISIN -a national programme for early warning, investigation and surveillance of healthcare-associated infection in France[END_REF] .

As part of HAI prevention strategy, a national point prevalence survey of HAIs has been organized every 5 years since 1996 [START_REF]The-French-Prevalence-Survey-Study-Group. Prevalence of nosocomial infections in France: results of the nationwide survey in 1996[END_REF][START_REF] Thiolet | Prévalence des infections nosocomiales, France[END_REF] . The most recent survey with available data was conducted in 2006; it involved 2,337 health care facilities and more than 350,000 patients [START_REF] Thiolet | Prévalence des infections nosocomiales, France[END_REF] .

Compared to the incidence surveys, which could represent long periods of surveillance, point prevalence survey data may not be appropriate for inter-facility comparison because they represent a "snapshot" with limited sample size and large confidence intervals. Nevertheless, many countries still use prevalence survey data to compare groups of facilities linked by geographic proximity or with similar patient case-mixes or cooperation programs 3 -5 .

In 2009, France created Regional Health Agencies (RHA) under the Hospital, Patients, Health and Territories Act. Two levels were defined to monitor health policy in France: a national program defines actions to reduce HAI infections and then at regional level, the Regional health agencies implement policies to ensure that health care provision meets the needs of the population. The objective is to improve coordination among ambulatory, hospitals and social care sector services while controlling national health expenditures [START_REF] Chevreul | France: health system review[END_REF] . Given this governance structure, comparison of HAI prevalence at a regional level would seem to be a relevant approach. Aggregated data, which avoid problems of limited sample size could permit regional level comparisons and identification of outlier regions, which could undergo reappraisal of their prevention strategies. However, before using HAI prevalence as a basis for inter-regional comparisons of the quality of infection control, the measurement should be adjusted for patients' characteristics and exposure to invasive devices [START_REF] Archibald | Hospital-acquired infections in the United States: the importance of inter hospital comparisons[END_REF][START_REF]The National Nosocomial Infections Surveillance (NNIS) System. Nosocomial infection rates for inter hospital comparison: limitations and possible solution[END_REF] and inter-facility as well as inter-regional prevalence heterogeneity resulting from sample size heterogeneity which may bias comparisons at the regional level. In this paper, a three-level hierarchical logistic regression model was implemented with a Bayesian non-parametric method to assess between-region heterogeneity.

The objective was to propose an original method of benchmarking regions according to their prevalence of healthcare-associated infections and to identify outlier regions.

METHODS

Data sources

Data were obtained from the 2006 national PPS database; they included information on 358,353 patients from 2,337 health care facilities. The participating facilities including regional hospitals, general hospitals, university hospitals, local hospitals and medicinesurgery and obstetrics clinics represented 83.3% of the HCFs and 93.6% of the inpatient beds in France. On a selected day in June, the 2006 national HAI point prevalence survey used a standardized protocol [START_REF] Thiolet | Prévalence des infections nosocomiales, France[END_REF] prepared by the French HAI surveillance network. Data on all patients hospitalized for more than 24 hours were entered into the database by each participating facility and aggregated and anonymized at regional and national levels [START_REF]RAISIN-Réseau-d'alerte-d'investigation-et-de-surveillance-des-infections-nosocomiales. Enquête nationale de prévalence des infections nosocomiales, juin 2006 -Protocole d'enquête[END_REF] . For each patient, the health care facility size, type and status were noted, along with the specialty of the receiving facility, the patient's characteristics, invasive device exposure and surgical procedure exposure. For those who were reported to have an infection, the location of the infection(s), the origin (imported or acquired in the current facility), the date of diagnosis and the microorganism(s) identified were documented.

Definition of Variables

An HAI was defined as an infection acquired during hospital care which was not present or incubating at the time of admission. The specific anatomical site definitions for HAI were standardized and adapted [START_REF]RAISIN-Réseau-d'alerte-d'investigation-et-de-surveillance-des-infections-nosocomiales. Enquête nationale de prévalence des infections nosocomiales, juin 2006 -Protocole d'enquête[END_REF][START_REF]et de la Solidarité. 100 Recommandations pour la surveillance et la prévention des infections nosocomiales[END_REF] from both the definitions of the Centers for Disease Control and Prevention [START_REF] Garner | CDC definitions for nosocomial infections[END_REF] (CDC) for acute care and the definitions of Mc Geer for long-term facilities [START_REF] Mcgeer | Definitions of infection for surveillance in longterm care facilities[END_REF] .

A binary outcome variable was indicated whether a given patient had an HAI on the day of the survey.

Explanatory variables included in the current analysis were classified into patient characteristics and invasive device exposure. Patient characteristics included sex, age, ward (medicine, surgery, gynecology, intensive care, rehabilitation, long-term care, cancer, and psychiatry), immunodeficiency, surgical procedure within the past 30 days, and McCabe score. The McCabe score is a three-group classification that defines the severity of a patient's underlying medical conditions as nonfatal (score 1), ultimately fatal (in about five years, score 2), or rapidly fatal (in about one year, score 3) [START_REF] Mccabe | Gram negative bacteremia. I. Etiology and ecology[END_REF] . It was specified in the model as an ordinal variable. Invasive devices included central venous or peripheral catheters (nominal variable with three levels: no catheter, peripheral catheter, and central venous catheter), urinary catheterization, and intubation (binary variables: presence vs. absence).

Exclusion criteria

Because facility level was accounted for patients who may have acquired an HAI from a referring facility were excluded (1% of the data, 20.7% of infected patients) to avoid the ambiguity concerning the accountable facility. Patients who had at least one missing value for the variables taken into account in our analysis (sex, age, ward, immunodeficiency, surgical procedure within the past 30 days and McCabe score) were also excluded (5.0% of patients).

Since our proposed method took into account facility effect, those facilities which had less than 5 patients were also excluded (five facilities with one patient, two facilities with two patients and one facility with three patients). Thus, the analyzed dataset consisted of 2,289 health care facilities and 336,858 patients.

Data analysis and statistics

A three-level hierarchical model was implemented, with patients (level 1) nested within facilities (level 2) and facilities nested within regions (level 3). To reduce the size of the dataset, the analysis was carried out in two steps (see details in the Appendices). In the first step, a single-level logistic regression model was used to predict the risk of infection per patient according to patients' characteristics and surgery procedures. These predictions were averaged at the facility-level to generate an offset variable for each facility, which was later introduced into the facility-level dataset (step 2).

The second step involved building the three-level hierarchical regression to model the HAI prevalence using aggregate data at facility-level. The three factors not taken into account in the previous model (central venous or peripheral catheterization, urinary catheterization, and intubation) were introduced treating the facility as a random effect in order to quantify the heterogeneity of these invasive device effects between facilities. The catheterization variable was designed to be nominal since those patients who had central venous catheter in most cases had also peripheral catheter. By reducing the number of the risk factors taken into account at the facility level, aggregation of the data was simplified.

We then focused on the facility-and region-level effects on prevalence. These effects were not assumed to be constant across facilities or regions and not assumed to be completely independent because of geographic proximity and cooperative programmes among facilities.

We assumed that facility effects were randomly drawn from a Student'st distribution which is equivalent to an assumption of exchangeability but is more flexible than the normal distribution.

At regional level, a non-parametric distribution was used, based on a Bayesian non-parametric method [START_REF] Ohlssen | Flexible random-effects models using Bayesian semi-parametric models: applications to institutional comparisons[END_REF] . We assumed that the real region effects were generated from a mixture of masspoints with each corresponding to a cluster of regions with similar prevalence. To generate this discrete distribution, we used a truncated Dirichlet process by placing a stick-breaking prior [START_REF] Ishwaran | Gibbs sampling methods for stick-breaking priors[END_REF] over the unknown distribution F and assuming a finite number of clusters (27 in our case). Given observed data, the posterior random-effect distribution over F was obtained using Gibbs sampling [START_REF] Gelman | Bayesian Data Analysis (2nd edn)[END_REF] . After a burn-in phase of 150 000 iterations, the posterior distribution of the random effects of each region was obtained from a large number of iterations (around 100 000). At each iteration, a rank was attributed to each region according to its randomeffect value. This gave the posterior distribution of the rank of each region.

RESULTS

Description of patients' and HCFs' characteristics by region

Among the 27 regions, 2,289 health care facilities, and 336,858 patients, the overall prevalence of patients with at least one HAI was 4.0%. The average number of patients per facility was 147 (min 5, max 2278). At regional level, HAI prevalence ranged from 2.1% (Corse) to 4.8% (Haute-Normandie) whereas at facility level, HAI prevalence ranged from 0 to 60% (3 HAI among 5 hospitalized patients) with a median of 4.1% (first and third quartiles: 1.4% and 6.8% respectively). The number of health care facilities per region varied from 2 in Guyane to 287 in Ile-de-France (Table I). Patients' characteristics also varied by region with the median of the regional average age being 63 years whereas the minimum was observed in region Territoires d'Outre-Mer (TOM) with an average of 40 years of age.

Effects of patients' characteristics and surgical procedure exposure

The results of the multivariate analysis at patient level illustrate the effect of patients' characteristics and exposure to surgical procedures on HAI (Table II). All the factors showed significant effects. Not surprisingly, there were large variations of the risk of HAI between patient wards with the odds of an HAI in intensive care unit about 7 times higher than those in the psychiatry unit.

Mean effect of invasive devices and heterogeneity of this effect between facilities

The use of an invasive device (catheterization or intubation) increased the odds of an HAI, but there was remarkably little inter-facility variation in the odds ratios for catherization or intubation (Table II). The medians of the OR distribution regarding the three studied factors were very similar between regions (Figure 1), yet in some regions, the heterogeneity of the ORs between facilities was more important than in other regions. For example, region Poitou-Charentes showed a median OR of 3.22 for the effect of urinary catheterization, which is really close to the overall average OR. Nevertheless, the third quartile of the distribution 3.57 was the highest of all the regions. We highlight the fact that for regions with few facilities, like Guyane (2 facilities) and TOM (5 facilities), the interpretation of ORs variability was not pertinent.

Inter-region comparison of HAI prevalence

The distribution of the regional prevalence of HAI identified clearly two clusters of regions.

The results of the model with the assumption of a mixture of two Gaussian distributions confirmed this hypothesis (Figure 2, a). One cluster included five regions had a lower adjusted prevalence than the other cluster of 22 regions (Figure 2,b). Hence, no regional outliers with unusually high adjusted prevalence were identified.

Furthermore, the ranking based on non-adjusted prevalence values appeared to be significantly different from the one obtained with adjusted values (Figure 2,b). For instance, based on the raw prevalence, Aquitaine was ranked sixth while after adjustment it became the first. We also noted that the adjusted ranks were given with very wide credible intervals (CI), which indicates a tremendous variation in rank estimations.

DISCUSSION

Within the context of a standardized national HAI surveillance system, inter-facility and interregion comparisons of HAI prevalence may help to improve coordination among health care services and to define strategies for infection control at regional level 3 -5 .

The simplest way to make such comparisons is a direct analysis of the distribution of HAI prevalence 17 . Some other approaches like control charts or the Standardized Infection Ratio are also widely used to monitor and compare quality between healthcare settings. However, most of these methods are either not risk-adjusted or performed poorly with small sample sizes and numerous risk factors [START_REF] Gustafson | Practical risk-adjusted quality control charts for infection control[END_REF] . The use of a multi-level regression model is considered more informative, as it permitted us to take into account inter-regional variability in the number of facilities among regions and the numbers of patients in facilities, as well as patientlevel factors. In this way, not only the inherent multilevel structure of the dataset are considered, but with a proper adjustment of the prevalence, the regional and facility effect can be identified seamlessly and the effect of invasive device which varies at facility-level can be studied more efficiently without too much misleading noises. When it comes to prediction, the uncertainty may therefore be more accurately assessed.

In modelling correlated data, random effects are routinely assumed to be normally distributed (e.g., the use of league tables and hierarchical models to compare education or healthcare providers [START_REF] Goldstein | League tables and their limitations: Statistical issues in comparisons of institutional performance[END_REF][START_REF] Normand | Statistical methods for profiling providers of medical care: issues and applications[END_REF] ). This assumption or normality presupposes that all the regions draw from the same population. However, the variability between regions may be wider than the variability allowed by a Gaussian model. Furthermore, assuming a normal distribution of the random effects tends to involve an excessive regression of the extreme values towards the overall mean [START_REF] Ohlssen | Flexible random-effects models using Bayesian semi-parametric models: applications to institutional comparisons[END_REF] . This can lead to spurious interpretations. For example, an improved ranking could simply reflect the shrinkage effect related to adjustment to reflect the normal distribution.

To avoid this, our analysis was based on a Bayesian non-parametric approach [START_REF] Ohlssen | Flexible random-effects models using Bayesian semi-parametric models: applications to institutional comparisons[END_REF] employing a truncated Dirichlet process to model the regional adjusted prevalence with a less restrictive random-effect distribution than the Gaussian one [START_REF] Walker | Bayesian nonparametric inference for random distributions and related functions (with discussion)[END_REF][START_REF] Burr | A Bayesian semi-parametric model for random effects meta-analysis[END_REF] . As with the Gaussian assumption, this approach protects against extreme estimates of HAI prevalence in regions with few facilities but also against excessive regression to the mean which allowed the identification of clusters of regions.

The heterogeneity of the facility effects was modelled using a Student's-t distribution known to be more tolerant to outliers than the normal one. Even so, the posterior distribution showed an obvious skewness, which supposes a skewed normal or Student's-t distribution [START_REF] Lee | Flexible parametric models for random-effects distributions[END_REF][START_REF] Lachos | Robust linear mixed models with skew-normal independent distributions from a Bayesian perspective[END_REF] . Due to technical limitations and as our objective was not to model inter-facility heterogeneity, the case of a skewed distribution was not considered.

Nevertheless, the ranking of the regions according to their adjusted prevalence should be interpreted with caution. With a given dataset, two different models can generate two different rankings and both of them may fit the data equally well [START_REF] Goldstein | League tables and their limitations: Statistical issues in comparisons of institutional performance[END_REF] . An over-interpretation of a set of rankings with wide confidence intervals can thus easily lead to unfairness. This is why we believe that it is more appropriate to keep, as main result, the existence of two clusters of regions that, according to our analysis, were clearly separated. Within each cluster, the median ranks of regions were very similar, which appears to be an inevitable consequence of attempting to rank units with similar performance.

The difference between these two clusters could be explained by different regional policies of healthcare, changes in regional surveillance programmes, or some other regional difference.

Although a relatively high or increased HAI prevalence may suggest a potential problem in the IC programme, and vice versa [START_REF] Suka | Epidemiological approach to nosocomial infection surveillance data: the Japanese Nosocomial Infection Surveillance System[END_REF] , the interpretation of our results is not that simple [START_REF] O'neill | Use of surveillance data for prevention of healthcare-associated infection: risk adjustment and reporting dilemmas[END_REF] .

First, although great efforts were made to ensure that all data were collected in a standardized way by trained and dedicated personnel, we cannot exclude the existence of methodological variations during the implementation of the protocol in some regions, particularly regarding HAI case finding and recording, which may have resulted in a lower PPS in some regions.

Furthermore, there are still unknown characteristics at facility and at patient level that could partially explain the differences in HAI prevalence which were not taken into account in this study (e.g., differences in patient acceptance policies between facilities, types of surgical interventions, patients' comorbidities). Thus, saying that risk-adjusted differences in outcomes can be entirely attributed to the level of healthcare quality would not be true.

It is noteworthy that the French north-south gradient, frequently identified in descriptive studies of various health indicators (cancer, cardiovascular disease) [START_REF] Drees | L'état de santé de la population en France : rapport 2008. La documentation française[END_REF] was found again in this cluster identification (Figure 2,c). This suggests that the level of individual morbidity considered here (i.e., the adjustment on only McCabe score) is probably not sufficient. The 2012 national point-prevalence survey is ongoing in France and in Europe, based partly on the 2006 French protocol. Improvements concerning the comorbidities were made in the new 2012 French protocol by adding information of cancer status [START_REF] Ecdc | Point-prevalence survey of healthcare-associated infections and antimicrobial use in European acute care hospitals[END_REF][START_REF]INVS : Enquête nationale de prévalence 2012 des infections nosocomiales et des traitements anti-infectieux en établissements de santé -Mai-Juin 2012 -Protocole / Guide de[END_REF] .

Furthermore, the obtained results may guide future action plans. For instance, in several regions such as Lorraine and Picardie, the effects of exposure to invasive devices were very similar between facilities suggesting certain homogeneity in care practices, whereas in other regions the heterogeneity of these effects was highlighted (Poitou-Charentes for urinary catheterization). This suggests a need for regional interventions to reduce the effect of invasive devices in some healthcare facilities. Our proposed approach could as well be applied to the data of the 2012 national survey (data available from 2013) to benchmark healthcare facilities and regional healthcare policies according to their quality of care and to evaluate the impact of regional infection control programs.

CONCLUSION

The present study illustrates use of a three-level hierarchical logistic regression model with a Bayesian non-parametric method to model between-region heterogeneity. In the context of the national French HAI prevalence survey, two clusters of regions regarding HAI prevalence were identified. This type of method could be more applied for benchmarking regional healthcare policies according to their quality of care and the estimation of facility-specific effects of invasive devices may orient future regional action plans. l'enquêteur. http://www.cclinparisnord.org/ENP/ENP2012/ENP_protocole2012.pdf Accessed June 7th 2012.

List of tables :

Table I -Description of the patients' characteristics at regional level and of their heterogeneity between regions -French national HAI PPS, 2006.

Characteristic

APPENDICES

The first step of the implementation consisted in the generation of an offset per facility corresponding to the part of HAI prevalence explained by the characteristics of the patients and their exposure to surgical procedures. A predicted conditional probability of HAI was obtained for each patient using a regression logistic model implemented with SAS software. y i was the binary outcome corresponding to the presence or absence of at least one HAI for the patient i.

ρ i represented the risk of HAI for the patient i.

The linear predictions omitting the intercept for the patients of each facility were averaged at the facility level to generate an offset value for each facility. n i j was the number of patients in facility j.

In the second step, a three-level hierarchical logistic regression model was built. To aggregate our dataset at facility-level, we used combinations of three extrinsic risk factors.

Since central venous or peripheral catheterization was the single three-level variable (all others being binary), twelve combinations of risk factors per facility were obtained. A Bayesian analysis was then used to estimate the parameters of the model.

+ Cathet + offset + + u = ) (ρ     2 1 logit Intercept n = offset j i i j   0 i ) ( logit 
y ijk was the number of patients with HAI for the combination k of the factors introduced in the model, in facility j and region i. n ijk was the total number of patients for the combination k, in facility j and region i.

ρ ijk was the probability of HAI. u j was the log of the HAI odds for the facility j. It corresponded to the facility random effect.

i was the random effect of the region i expressed as a log of odds ratio.

Regional random effect distribution:

At the region-level, a non parametric distribution generated by a truncated Dirichlet Process (DP), was used to model the random effect distribution F. The main idea was to consider that the studied population was a mixture of regional clusters with various characteristics which can be represented by different normal distributions to provide larger tolerance to the outliers.

The DP was the prior placing over F:   0 , ~F DP F  with F 0 corresponding to a baseline distribution and  to a positive real parameter which measures the concentration of F around F 0 . The higher the value of α, the closer F is to F 0 . We used the stick-breaking construction of the truncated DP which assumed that our dataset was generated from a finite mixture of mass-points. It consisted in drawing N values of  from the baseline normal distribution F 0 and N-1 values of  from a beta distribution with parameters 1 and :

Facilities' random effect distribution:

At the facility-level, a Student's-t Distribution was used to model the random effect distribution, in order to provide a more flexible distribution than the normal one.

For practical purposes, we considered the t-distribution parameterised as a ratio υ V Z

/

, where Z is normally distributed with mean 0 and variance 1, V has a chi-square distribution with υ degrees of freedom.

For all the parameters, priors were chosen to be weak.

Figure1

  Figure1 Distribution of invasive device exposure effects as estimated by the hierarchical logistic regression model at facility level, according to each region. Panel a: peripheral or arterial catheterization. Panel b: central venous catheterization. Panel c: urinary catheterization. Panel d: intubation (French national HAI PPS, 2006).

to reach the probability 1 . 1 

 11 the prior random effect values of the N mass-points (random effect of region cluster k) and N corresponded to the maximum number of mass-points fixed in our case at 27, the total number of regions. The 1  k s were used to obtain the prior probabilities p k (the probability of a region belonging to the k th cluster) assigned to the mass-points as follow: for k=1 This generation can be considered as a stick-breaking prior, as we can imagine that a stick of length 1 is broken into 27 parts. At time k, a proportion  k of the stick remaining length is broken off, leaving a length equal to  . In this way, the length of the k th part broken off is p k .

Table II -

 II ( a ) Adjusted odds ratios of healthcare associated infection for patients' characteristics and surgery procedure exposure, estimated by a single-level multiple logistic regression model; ( b ) Estimation of the mean adjusted Odds Ratios of health-care associated infection and of their heterogeneity between facilities, for the factors relative to invasive device exposure. Results of the three-level logistic regression model-French national HAI PPS, 2006.

			Mean effect	Inter-facility variation
	Risk factors	OR	95% CI	25 th percentile	75 th percentile
	Age a (per five-year increment)	1.04	1.04 -1.05	-	-
	Male vs. Female a	1.12	1.08 -1.17	-	-
	McCabe Score a			-	-
	1	1.00	-	-	-
	2	1.87	1.82 -1.92	-	-
	3	3.50	3.31 -3.69	-	-
	Immunodeficiency a	1.34	1.27 -1.41	-	-
	Hospitalisation unit a			-	-
	Psychiatry	1.00	-	-	-
	Medicine	1.52	1.39 -1.67	-	-
	Gynecology	0.67	0.57 -0.80	-	-
	Intensive care	8.38	7.50 -9.37	-	-
	Rehabilitation	2.03	1.85 -2.24	-	-
	Long-term care	1.59	1.44 -1.76	-	-
	Cancer	2.29	2.01 -2.61	-	-
	Surgery	1.77	1.60 -1.95	-	-
	Surgery procedure exposure a	1.64	1.56 -1.72	-	-
	Intubation b	1.97	1.78 -2.18	1.95	1.98
	Catheterization b
	None	1.00
	Peripheral or arterial	1.55	1.45 -1.65	1.44	1.58
	Central venous	3.19	2.95 -3.45	3.11	3.22