Oriol Avino-Salvado

Hervé Morel

Cyril Buttay
email: cyril.buttay@insa-lyon.fr

Denis Labrousse

Stéphane Lefebvre

Threshold voltage instability in SiC MOSFETs as a consequence of current conduction in their body diode

An antiparallel PiN diode is present in the structure of the SiC Power MOSFET. For simplicity, increased module current ratings, and cost reasons, it sometimes may be preferable to use this diode rather than to add an external Schottky diode to the MOSFET. However, SiC PiN diodes are sensitive to defects in the SiC crystal, and their performance can degrade rapidly.

In this paper, we assess the effect of several stresses on the characteristics of a SiC MOSFET, focusing on the stresses which involve the body diode. We demonstrate that using the body diode does not degrade its performances. However, it may produce a noticeable shift in the threshold voltage of the MOSFET. While this behaviour is observed for MOSFETs from one manufacturer, it is not exhibited for others, which indicates that it may be a consequence of some particular design or manufacturing parameters.

Introduction

SiC power devices are now commercially available from many manufacturers [START_REF] Boomer | Body of knowledge for silicon carbide power electronics[END_REF]. They offer dramatic improvements in performance over the ubiquitous Si power devices [START_REF] Millán | A Survey of Wide Bandgap Power Semiconductor Devices[END_REF], especially regarding high voltage or high temperature [START_REF] Raynaud | Comparison of high voltage and high temperature performances of wide bandgap semiconductors for vertical power devices[END_REF]. SiC technology is progressing rapidly, and although devices still suffer from some reliability issues, they continuously improve.

Many issues regarding manufacturing yield or device reliability can be traced back to defects in the SiC crystal structure, such as screw dislocations [START_REF] Neudeck | Study of bulk and elementary screw dislocation assisted reverse breakdown in low-voltage ( <250 V) 4H-SiC p+-n junction diodes. i. dc properties[END_REF] or micropipes [START_REF] Kimoto | Performance limiting surface defects in SiC epitaxial p-n junction diodes[END_REF]: when they do not alter a device's performance so much that it becomes useless, defects can cause a gradual drift in its performance. An example is the increase in forward voltage of PiN diodes as a result of the propagation of defects such as stacking faults [START_REF] Singh | Reliability and performance limitations in SiC power devices[END_REF]. Similar observations were made on SiC bipolar transistors [START_REF] Muzykov | Physical phenomena affecting performance and reliability of 4H-SiC bipolar junction transistors[END_REF]. These defects also cause increased leakage current, and therefore strongly affect high voltage devices [START_REF] Chen | Effects of different defect types on the performance of devices fabricated on a 4H-SiC homoepitaxial layer[END_REF]. SiC wafers with reduced defect density are now available, with less than 1000 dislocations/cm 2 over a 150 mm wafer [START_REF] Boomer | Body of knowledge for silicon carbide power electronics[END_REF], so crystal quality is no longer considered as a major issue which could impede the commercial success of SiC.

Another set of issues affects SiC MOSFETs. It is related to the electron mobility in the channel [START_REF] Singh | Reliability and performance limitations in SiC power devices[END_REF] and the reliability of the oxide layer which forms the MOS structure [START_REF] Lelis | Basic Mechanisms of Threshold-Voltage Instability and Implications for Reliability Testing of SiC MOSFETs[END_REF]. The consequences are a higher on-state resistance (in the case of poor channel mobility) and a gradual drift in the threshold voltage of the MOSFET. In some extreme cases, such drift could cause the device to become normally-on [START_REF] Yang | High temperature gate-bias and reverse-bias tests on SiC MOSFETs[END_REF]. As with the crystal defects, recent improvements in the technology have largely mitigated these issues, and SiC MOSFETs have de facto become the major SiC controlled switch. This is illustrated by the discontinuation of alternative devices such as Infineon's JFET of Fairchild's BJT. With the improvement in technology, the SiC devices have now reached a point where they can be used in real applications, and the reliability studies have shifted from simple test setups to more complex configurations. An example is the body diode of the MOSFET. In inverter topologies, an antiparallel diode is usually connected to each transistor so current can flow during the "dead-time" periods (when both transistors of a half-bridge are in the off-state). Classically, a SiC Schottky (or Merged PiN-Schottky) diode is used, because it has a low forward voltage (V F ≈ 1 V) and because it does not suffer from material degradation, as presented above. However, this requires two components per switch (a MOSFET and a diode). As the power MOSFET structure already contains a PiN diode (Fig. 1), a possibility is to use it instead of an external diode. This saves space (one die per switch instead of two) and money, at the cost of increased losses: PiN diodes have a higher forward voltage drop (V F ≈ 3 V instead of 1 V for Schottky diodes) and recovery charge (although not nearly as much as comparable Si diodes). However, as the dead-time period is short, usually a few percents of the switching period, the effect remains small. It is worth noting that a SiC MOSFET with an integrated Schottky diode was introduced in [START_REF] Kusumoto | Reliability of Diode-Integrated SiC Power MOSFET (DioMOS)[END_REF], but in most cases, the internal (body) diode has a PiN junction as described in Fig. 1.

A recent study [START_REF] Fayyaz | Body diode reliability investigation of sic power mosfets[END_REF] investigated the degradation of the PiN junction of a SiC MOSFET in realistic use cases (a 3-phase inverter). Only a very small change in forward voltage drop was observed (increase as well as decrease, indicating that stacking faults were not the predominant degradation mechanism at work). Overall, it was found that the body diode could be used as a functional element, in place of external diodes.

Here, we present the results of a similar investigation, where we not only monitored the changes in the forward voltage drop of body diodes (V F ) but also in the threshold voltage of the corresponding MOSFET (V th ). As in [START_REF] Fayyaz | Body diode reliability investigation of sic power mosfets[END_REF], we observe no significant change in V F due to current flowing in the body diode. However, we observe a strong shift in V th in at least one manufacturer's parts, which has not been reported so far.

The next section of this article presents the test protocols based on static (i.e. no switching) or dynamic (more realistic case) operations. The corresponding results are presented and discussed in section 3.

Test Protocol

Except when explicitly mentioned, all tests presented here were performed on C2M0080120D SiC MOSFETs (Wolfspeed), which are rated at 1200 V and 80 mΩ. The devices were purchased in two batches (2014 and 2016): the first batch was used for the HTRB/CMB tests (see below), while the second was used for all the other tests (Q1, Q3, Q4 and HTGB). For the sake of simplicity, these devices are referred to as "Wolfspeed" MOSFETs in the text. All characterizations were performed using a Keysight B1505A curve tracer, at room temperature, 1 h after the end of the corresponding test (to allow enough time for the devices to cool down, but soon enough to minimize relaxation effects).

Static Tests

Five static test circuits are used: the first two are the "High Temperature Reverse Bias" (HTRB) and "High Temperature Gate Bias (HTGB) (Fig. 2, top), described in a JEDEC standard [START_REF]JESD22-A108D, temperature, bias and operating life[END_REF], which are classically used to evaluate the effects of a continuous bias on power devices. For the HTRB, a negative voltage source (here -8 V, corresponding to a value close to the absolute minimum voltage, -10 V, allowed in its datasheet [START_REF] Wolfspeed | D Silicon Carbide Power MOSFET, C2M MOSFET Technology[END_REF]) is connected to the gate and source of the Device Under Test (DUT), and a positive voltage is applied between drain and

DUT -8 V 540 V DUT -8 V/22 V DUT 20 V 16 A DUT 20 V 16 A DUT -8 V 10 A Q1 Q3 Q4
HTRB HTGB

V gs + V gs - source (here 540 V, which corresponds to the bus voltage of the aeronautic HVDC standard [START_REF] Nya | Benefits of higher voltage levels in aircraft electrical power systems[END_REF]). For the HTGB, the drain and source terminals are shorted, and two gate levels were investigated here: -8 V or 22 V (a value very close to the absolute maximum allowed). In these tests, no current flows through the DUT.

I d + I d -
The three remaining test circuits are also depicted in Fig. 2 (configurations Q1, Q3 and Q4). They represent the three configurations which result in current flowing through the drain of the MOSFET: either the transistor is in the on-state (V gs =20 V, the recommended value in the datasheet) so the current flows through its channel (I CH path in Fig. 1) in both directions (Q1 or Q4), or the transistor is in the off-state (V gs =-8 V), so the current flows through its body diode (I BD in Fig. 1). This latter configuration corresponds to Q3 in Fig. 2.

All tests are run at the maximum junction temperature allowed by the manufacturer (150 °C), to accelerate any degrada-

V dd =540 V M1 DUT I sw =10 A V gs = -8 V R V p1 R Figure 3: Circuit diagram of the CMB test.
tion mechanism which may occur. For the HTRB test, where no power is dissipated, this simply requires the DUT to be attached to a hot plate at a set temperature of 150 °C. On the contrary, for the other tests from Fig. 2 (Q1, Q3 and Q4), the DUTs dissipate a large amount of power, because of their large drain current, resulting in a noticeable self-heating. To ensure that they still operate at T J =150 • C, they are mounted on temperature-regulated plates. Their thermal resistance "junction-to-plate" (R T h JP ) is then characterized using the drain current of the DUT in the saturation region (V gs =6 V, V ds =8 V), both to generate selfheating of the device, and as a Temperature-Sensitive Electrical Parameter (TSEP, [START_REF] Avenas | Temperature measurement of power semiconductor devices by thermo-sensitive electrical parameters -a review[END_REF]) to monitor the corresponding temperature rise. R T h JP is found to be 2.19 K/W. During the static tests, the drain current is set to 10 A when it flows through the body diode (quadrant Q3) and 16 A when it flows through the channel of the MOSFET (quadrants Q1 and Q4). Both conditions correspond to a power dissipation P=37 W when the DUT operates at T J =150 °C. As a consequence, the temperature-regulated plates are set to a temperature T P of:

T P = T J -P × R T h JP = 68.3 • C (1)
A summary of the test conditions is given in Tab. 1. Tests in the Q3 quadrant were run for two values of V gs (-5 V and -8 V) as this parameter was found to have a strong effect on the results.

Dynamic Tests

The circuit diagram used for the dynamic tests is presented in Fig. 3. As can be seen in the theoretical waveforms presented in Fig. 4, the DUT experiences both high voltage when M1 is on, and high current (through its body diode, as the DUT remains off) when M1 is off. This operating mode is referred to as "Chopper-Mode Bias" (CMB).

An arbitrary switching frequency of 20 kHz was chosen for the results presented here, with a 50 % duty cycle ratio between M1 and the DUT. As for the static tests, the devices are attached to a temperature-regulated plate, set at the temperature required to achieve T J =150 • C. The temperature of the plate is calculated using eq. ( 1), where P is the power dissipated in the DUT when M1 is off.

Results and Discussion

Static tests

Fig. 5 presents the measurements of V F for quadrants 1, 3 and 4. No significant evolution is observed over 100 h: the changes in V F during the tests are comparable to the experimental scattering, except for Q1, which show a slightly clearer trend, with a 2 % drop in V F after 15 h. For the HTRB and HTGB tests (not presented here because of space limitations), V F exhibits no significant change. As for V F , the evolution of R ds on was also measured, but it is not presented here as this parameter shows no significant change, for any of the test conditions (HTRB, HTGB, Q1, Q3 or Q4).

The evolution of V th is plotted in Figs. 7 and6. It remains essentially constant for the HTRB, HTGB, Q1 and Q4 tests. For Q3, however, V th shows a dramatic drop over time (-40 % in 60 h) when the MOSFETs are driven with V gs = -8 V. The tests were repeated on 3 other devices, with V gs = -5 V, a value which still ensures proper blocking of the transistor (initial V th values were found to be 1.12 -1.45 V). The drop in V th then appears to slow down considerably (orange line in Fig. 7).

Q3 corresponds to the condition in which the biasing current flows through the body diode of the MOSFET (Fig. 2); in Q1 and Q4 the current flows through the channel, while HTRB and HTGB correspond to a voltage-only bias. However, contrary to what was expected, no degradation of the body diode is observed (increase in V F caused by the propagation of stacking faults). Instead, the devices exhibit a large drop in V th , for the Q3 test condition. Such decrease is undesirable, as it makes the device harder to block, and makes it more vulnerable to electromagnetic interference (which might cause spurious turn-ons).

Dynamic Tests

Using the CMB setup (Fig. 3), devices from three manufacturers (Wolfspeed, ST Microelectronics and Rohm Semiconductors) were tested. The results are given in Fig. 8 for V th (V F and R ds on were also monitored, but as with the static tests, they show no significant evolution). For Wolfspeed devices, a large V th drop is observed (Fig. 8). On the contrary, almost no change is observed for ST's or Rohm's MOSFETs. A comparison between the V th drop observed using the CMB test setup and the static Q3 test setup is presented in Fig. 9. Note that the x-axis corresponds to the time during which current actually flows through the body diode (i.e. 100 % of the test time for Q3, and 50 % of the test time for CMB). Tests were performed using two values of V gs : -5 and -8 V. In both cases, a good agreement is found: most of the V th drop can be associated with the use of the body diode of the DUT.

V F [normalized] Vgs = -8 V Vgs = -5 V Q1 Q3 Q4 V gs + V gs - I d + I d -
In an attempt to identify the cause behind the differences in behaviour among components from different manufacturers, additional investigations were conducted. Fig. 10 presents the results of HTGB tests performed on devices identical to those used for Fig. 8. Contrary to the data presented in section 3.1, these tests were run at T J =200 • C, to strongly accelerate ageing. It should be noted that while this temperature is within the allowed range for ST's MOSFETs [17], it exceeds Rohm's (175 • C, [START_REF] Semiconductor | Sct2080ke datasheet[END_REF]) or Wolfspeed's (150 • C, [START_REF] Wolfspeed | D Silicon Carbide Power MOSFET, C2M MOSFET Technology[END_REF]). The high tem- perature was nonetheless required to observe any difference in behaviour: as can be seen in Fig. 10, while all devices exhibit a small shift in V th (≈ 6 %) within 50 h of test, MOSFETs from ST and Rohm tend to stabilize over time at around a ≈ 10 % increase in V th . On the contrary, Wolfspeed MOSFETs exhibit a continuous increase in V th , with 30 % more over 1000 h. Only positive-bias HTGB tests were performed at 200 • C, because of limited sample availability.

Vgs = -8 V Vgs = -5 V Q1 Q3 Q4 V gs + V - I d + I d -
Additionally, cross sections were performed on devices from the three manufacturers. They were prepared mechanically (cut with a diamond saw, and polished using diamond particles down to 1 µm, followed by 50 nm oxide suspension), and observed using a scanning electron microscope (SEM, Vega 3, Tescan). It is found that all three devices exhibit a planar gate structure, with some differences in the size of their features. The resolution of the images is limited, and very fine details such as oxide thickness cannot be observed clearly. There is, however, no obvious geometric parameter which could be correlated with the observed shift in V th .

Discussion and Conclusion

Recent SiC MOSFETs were tested, initially to evaluate the extend of know degradation mechanisms in static tests and realistic (switching) conditions: the drift in V th caused by the voltage bias of the gate, and the increase in V F caused by the propagation of stacking faults in the crystal.

While none of these degradation mechanisms was found to have significant effects in our test conditions (this shows the progress of the SiC MOSFET technology in recent years), we observed a new degradation mechanism: the use of the body diode resulted in a strong drop in V th . In summary, the following conclusions can be drawn from the measurements presented above:

• for the Wolfspeed MOSFETs, static (Q3) and dynamic (CMB) tests have the same consequence, a large V th drop occurring when current flows through the body diode, with comparable time constants in the drop;

• this drop depends on V gs , with a larger change in V th with V gs = -8 V (> 40 % after 80 h) than with V gs = -5 V (20 %);

• in the absence of current biasing of the body diode (i.e HTRB/HTGB conditions, with only voltage bias, or Q1/Q4 conditions, with current biasing of the channel), no noticeable V th drift is observed at T J =150 • C;

• increasing T J up to 200 • C yields a noticeable V th increase in HTGB conditions, but it remains much slower than that observed for Q3 or CMB tests (hundreds of hours vs. tens of hours)

• devices from different manufacturers (here ST Microelectronics and Rohm Semiconductor) appear to be much less sensitive to this V th drift, which indicates that the issue can be mitigated through improvements in the device design or manufacturing processes;

• the V th drift in Q3 or CMB conditions is correlated with a larger drift in HTGB conditions at 200 • C (i.e. it affects Wolfspeed's MOSFETs much more than Rohm's or ST's), but the physical cause remains unknown. Cross sections of the three devices do not show fundamental differences in the geometry of these transistors.

The root cause of this mechanism is not fully understood yet: as depicted in Fig. 1, body diode current and channel current flow through different paths. Charge injection in the channel (which would cause a drift in V th ) is therefore unlikely when using the body diode. For the same reason, it is not expected that the channel area would reach a higher temperature when using the body diode, as self heating occurs in a slightly different location. The power dissipation (and therefore the junction temperature) is identical for Q1, Q3 and Q4 tests, so a possible error in the estimation of the junction temperature would affect Q1, Q3 and Q4 test results identically, and yet the V th drift is very different (see Fig. 7). This degradation is relatively severe, as it may make some devices impossible to turn off, or make them very sensitive to voltage transients on the gate within a few tens of hours of operation.

Qualifications test standards, such as HTRB or HTGB, were developed mainly for Si devices. With the development of SiC MOSFETs, some additional tests should be introduced. Here, we describe a very simple static test circuit (the Q3 test setup) which gives results which are consistent with those obtained in switching conditions. As shown in the paper, this test setup is capable of exposing weaknesses which were not visible after HTGB or HTRB tests, and yet may eventually cause some transistors to fail.

Figure 1 :

 1 Figure 1: Cross section of a power MOSFET, showing the I BD (body diode) and I CH (channel) current paths.

Figure 2 :

 2 Figure 2: The five static test circuits: voltage bias on the drain (HTRB), negative or positive gate bias on the gate, with drain and source shorted (HTGB), and three tests with current bias (Q1 to Q4).

Figure 4 :

 4 Figure 4: Principle of the CMB test bench: When M1 is in the off-stage, the current flows through the body diode of the DUT. When it is on, the full voltage (540 V) is applied between drain and source of the DUT.

Figure 5 :Figure 6 :

 56 Figure 5: Evolution of the forward voltage of the MOSFETs (normalized for t=0, measured at I d = 10 A, V gs = -8 V, T J = 25 • C) during the current-bias tests described in Fig. 2 (3 samples per test for Q1, Q3 and Q4, 10 samples for HTRB). All tests were run at T J =150 • C.

Figure 7 :Figure 8 :

 78 Figure 7: Evolution of the threshold voltage of the MOSFETs (normalized for t=0, measured at I d =100 µA, V dd = 1 V, T J = 25 • C) during the static tests described in Fig. 2 (3 samples per test). All tests were run at T J =150 • C.

Figure 9 : 4 VC2M0080120DFigure 10 :

 9410 Figure 9: Evolution of the V th as a function of time for the Static (3 rd quadrant) and CMB tests (3 samples per test), data normalized for t=0, measured at I d =100 µA, V dd = 1 V, T J = 25 • C. All tests were run at T J =150 • C.

Table 1 :

 1 Summary of the static tests conditions

	Test	samples V gs (N) (V) (A) (V) (°C) I ds V ds T j	Power (W)
	HTRB	10	-8	0	540 150	0
	HTGB	1	-8	0	0	150	0
	HTGB 1 st quadrant 3 rd quadrant 3 rd quadrant 4 th quadrant	3 3 3 3 3	22 20 -8 -5 20 -16 0 16 -10 -10	0 0 0 0 0	150 150 150 150 150	0 37.3 37.3 37.3 37.3

Acknowledgement

This work was funded as part of the GENOME-PREMICES initiative, in the scope of the "Programme d'Investissements d'Avenir" (Investment for the future) of CORAC (French Research Council for Civil Aeronautics).


Mylène Duval
email: mylene.duval@toulouse.inra.fr

Christèle Robert-Granié

SAEM-MCMC: some criteria

Keywords: non-linear mixed models, maximum likelihood, stochastic algorithm, Metropolis-Hastings, SAEM-MCMC. 1

come

Introduction

The main interest of this paper is to obtain good parameter estimates using maximum likelihood estimation in non-linear mixed effects models. Let us recall some information about the linear mixed models. The Expectation-Maximization algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] is a very famous tool used for parameter estimation in this class of models. Let y denote the vector of observed data, φ the vector of unobserved (or missing) data, θ the vector of parameters, p(y; θ) the incomplete likelihood and p(y, φ; θ) the complete likelihood. In the general case, since p(y; θ) is not in a closed form or is hard to compute, it seems to be difficult to maximize p(y; θ). In this sense [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] present the EM algorithm, which tends to maximize E log p(y, φ; θ)|y, θ . At iteration k of the algorithm, there are two steps: the E-step computes the conditional expectation of the complete log-likelihood, noted Q k (θ, θ (k-1) ) equal to E log p(y, φ; θ)|y, θ (k-1) and the M-step determines θ (k) as maximizing Q k (θ, θ (k-1) ). [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] and [START_REF] Wu | On the convergence properties of the EM algorithm[END_REF] proved the convergence of the sequence (θ (k) ) k towards a stationary point of the observed likelihood under general regularity conditions. In the case where we study non-linear mixed models, the E-step leads to an integral that has no closed-form solution. Several approximations of the incomplete log-likelihood have been proposed: the linearization procedure [START_REF] Sheiner | Evaluation of methods for estimating population pharmacokinetic parameters. I. Michaelis-Menten model: routine clinical pharmacokinetic data[END_REF], the LME approximation [START_REF] Lindstrom | Nonlinear mixed effects models for repeated measures data[END_REF], the Laplace approximation [START_REF] Wolfinger | Laplacian's approximation for nonlinear mixed models[END_REF][START_REF] Vonesh | A note on Laplacians approximation in nonlinear mixed effects models[END_REF]. Some of these methods are available on statistical software like the NLME procedure in S-PLUS software (Pinheiro and Bates, 1995a) or the NLMIXED procedure of SAS. Since errors can be large in the approximation of the observed log-likelihood (Davidian and Giltinian, 1995;Pinheiro and Bates, 1995b;[START_REF] Lindstrom | Nonlinear mixed effects models for repeated measures data[END_REF], some exact methods based on Monte Carlo methods have appeared. [START_REF] Wei | A Monte Carlo implementation of the EM algorithm and the Poor's Man's data augmentation algorithms[END_REF] proposed the MCEM algorithm in which the expectation of the E-step is estimated with a mean of some simulated samples from the exact conditionnal distribution of φ|y, θ. The MCEM algorithm requires an increase of the number of simulated data in order to have some accuracy, and so it is highly time consuming. For instance, [START_REF] Booth | Maximizing generalized linear mixed models likelihoods with an automated Monte Carlo EM algorithm[END_REF] report some results from a study on a real data set: they simulated around 60,000 samples for the final iteration. There exists a variant of this method based on the stochastic approximation method of [START_REF] Robbins | A stochastic Approximation Method[END_REF] which promises convergence with fewer simulations: the SAEM algorithm [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF]. When the conditionnal distribution of the missing effects given the observations, that is to say p(φ|y, θ) is unknown, [START_REF] Walker | An EM algorithm for nonlinear random effects models[END_REF] proposes an EM algorithm with importance sampling to estimate the conditional expectations. In his paper, [START_REF] Walker | An EM algorithm for nonlinear random effects models[END_REF] tells that importance sampling is more efficient than other algorithms that use dependent Markov chains to evaluate an expectation under an unknown distribution. In another way, Kuhn andLavielle (2004, 2005) presented a similar method based on the SAEM algorithm: the SAEM-MCMC algorithm, available on the MONOLIX group website (http://www.monolix.org/). In this method, the E-step of EM is replaced by a stochastic approximation and the simulated sample (which is under the unknown conditionnal distribution of φ|y, θ) used for the stochastic approximation is simulated with a Metropolis-Hastings algorithm [START_REF] Robert | Monte Carlo Statistical methods[END_REF]. In this method, several parameters need to be calibrated in order to better estimate the vector of parameters. For example, the Metropolis-Hastings algorithm based on the Monte Carlo and Markov Chains methods, simulates a markov chain using an instrumental distribution. We need to choose the instrumental distribution and the length of the markov chain. The aim of this study was to present some criteria which determine the parameters we need to fix before running the SAEM-MCMC algorithm, and we discuss if this choice of parameters is relevant or not. The paper is organized as follows: In section 2 we introduce the model and the SAEM-MCMC algorithm in its general version. In section 3 we propose some criteria to determine the parameters to run the algorithm. In section 4 we compare our SAEM-MCMC algorithm with the SAEM-MCMC algorithm in its general version applied with different sets of parameters on the well known orange tree data set and on a simulated data set. Section 5 is devoted to the development of the method when the complete data likelihood p(y, φ; θ) does not belong to the curved exponentiel family. Finally a validation of this method by simulation is presented.

2 The nonlinear mixed effects model

The model

We consider the following model:

y ij = f (z ij , φ i , β) + g(z ij , φ i , β, α) ε ij , ∀i ∈ {1, ..., N } ∀j ∈ {1, .., n i },
where y ij is the jth observation of subject i, N is the number of subjects, n i the number of observations of subject i and (z ij ) ij are known covariates. The ε ij 's are supposed to be independent identically distributed centered Gaussian random variables, independent of the φ i , with variance σ 2 . We note φ i (k × 1) the vector of individual random parameters of function f and g. It is modeled by:

φ i = A i µ + η i
where the η i ∼ N (0, Γ) are independent, µ (c × 1) is an unknown vector of fixed effets, the individual matrix A i contains covariates and is assumed to be known, and Γ is the covariance matrix of η i . β (p × 1) corresponds to the vector of the unknown other fixed parameters. The variance function g is dependent on f and a parameter vector α. In general g = f α , reflecting the possible character of intra individual variability.

The aim of this study was to estimate the complete vector of unknown parameters θ = (β, µ, Γ, σ 2 , α) by maximum likelihood estimation. In the case of a linear model, that is to say when f and g are linear in φ, the estimation of θ can be treated with the analytic EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. However, a non-linear function is often more suitable for modeling the physical problems but requires a specific approach for estimating the parameters because the Expectation step of the EM algorithm can not be in closed form. [START_REF] Kuhn | Coupling a stochastic approximation version of EM with a MCMC procedure[END_REF] propose a method based on the SAEM algorithm coupled with a Monte Carlo Markov Chains method, especially the Metropolis-Hastings algorithm. In this method, some parameters need to be calibrated and we propose some criteria in this sense.

The SAEM algorithm

The Stochastic Approximation version of the EM algorithm was proposed by [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF]. It consists in replacing the E-step of the EM algorithm by two steps: a simulation step of the missing data under the conditional distribution of φ|y, θ and a Stochastic Approximation step.

Given θ (k-1) , estimated parameter value of θ at iteration k -1

• Simulation step: Draw φ (k) under the conditional distribution p(.|y, θ (k-1) ) • Stochastic Approximation step: update Q k (θ) according to Q k (θ) = Q k-1 (θ) + γ k log p(φ (k) |y, θ (k-1) ) -Q k-1 (θ)
where (γ k ) k is a decreasing sequence of positive numbers

• the M-step is the same as the EM algorithm:

θ (k) = arg max θ Q k (θ)
If we assume that the complete data likelihood p(y, φ; θ) belongs to the curved exponentiel family, then we can write it by: p(y, φ; θ) = exp -Φ(θ) + S(y, φ), ψ(θ) where ., . denotes the scalar product and S(y, φ) is the minimal sufficient statistic of the complete model. Then the Stochastic Approximation step of the SAEM algorithm is reduced to compute:

s k = s k-1 + γ k S(y, φ (k-1) ) -s k-1
The maximization step of the SAEM algorithm consists in computing

θ (k) = Arg max θ -Φ(θ) + s k , ψ(θ)
Properties about the convergence of the sequence (θ (k) ) k towards the maximum likelihood under mild conditions are presented in [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF] in the case of the curved exponentiel family.

Remark In order to have the convergence of the sequence (θ (k) ) k , the serie (γ k ) k has to be chosen such that each γ k must belong to [0,1] and the series γ k must diverge, γ 2 k must converge. [START_REF] Kuhn | Coupling a stochastic approximation version of EM with a MCMC procedure[END_REF] propose to take the sequence (γ k ) k such that γ k = 1 for 1 ≤ k ≤ K and γ k = (k -K) -1 else, where K is an integer that can be fixed between 50 and 100. In practice since the sequence (γ k ) k is decreasing quickly, the choice of K seems to be very important and so θ (K) must be close to the maximum likelihood of θ to be sure that the sequence converges towards the maximum likelihood. In section 3 we present a criterion to fix K.

The SAEM-MCMC algorithm

Since (φ i |(y ij ) j , θ) i are independent random variables, we can generate them independently. In the general case, the distribution of (φ i |(y ij ) j , θ) i is not in a closed form, so the Simulation step of the SAEM algorithm cannot be directly performed. In this sense [START_REF] Kuhn | Coupling a stochastic approximation version of EM with a MCMC procedure[END_REF] proposed to combine the SAEM algorithm with a MCMC procedure: the Metropolis-Hastings algorithm [START_REF] Robert | Monte Carlo Statistical methods[END_REF]. The results of convergence of this method of estimation, called SAEM-MCMC, are proposed in [START_REF] Kuhn | Coupling a stochastic approximation version of EM with a MCMC procedure[END_REF]. This algorithm produces an ergodic Markov chain with stationnary distribution p(φ i |(y ij ) j , θ). At this stage some parameters need to be fixed before running this algorithm: the length of the chain (noted itM C) and the instrumental distribution. In this paper, we studied another instrumental distribution than the one proposed by [START_REF] Kuhn | Maximum Likelihood estimation in non linear mixed effects models[END_REF] and the one proposed in the Monolix software. The algorithm proposed is presented in detail in Section 3. As suggested in the user guide of Monolix software, we can improve the convergence of the SAEM-MCMC algorithm by running L independent Markov chains and then by doing the Stochastic Approximation on the mean of the (S(y, φ

(k,l,itM C) i
)) l=1,...,L , where φ (k,l,t) i corresponds to the tth iteration of the lth chain at the kth iteration of the SAEM-MCMC algorithm. In fact if L is large enough,

1 L L l=1 S(y, φ (k,l,itM C) i
) is a good estimator of E S(y, φ

(k,l,itM C) i
)|y, θ . To obtain a better approximation, we prefer doing moreover the mean on several iterations intra chains. However the first iterations obtained in the simulated markov chains with the Metropolis-Hastings algorithm are not under the stationary distribution, so we define a parameter (noted burn), which corresponds to the length of the burn-in period in the Metropolis-Hastings algorithm. Then the Stochastic Approximation step becomes:

s k = s k-1 + γ k 1 L(itM C -burn) L l=1 itM C t=burn+1 S(y, φ (k,l,t) ) -s k-1
where L is the number of independent chains, itM C the length of the chains and burn the length of the burn-in period within chain in the Metropolis-Hastings algorithm.

Estimations of standard errors and log-likelihood

Let θ be the estimator of θ at convergence, and we denote ∂ θ (∂ 2 θ ) the differential (the hessian) with respect to θ. Using the Louis formula [START_REF] Louis | Finding the observed information matrix when using the EM algorithm[END_REF] and following [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF], we approximate the observed information matrix of θ in the following way:

At the Stochastic Approximation step, we calculate

∆ k = ∆ k-1 + γ k [∂ θ log p(y, φ (k) ; θ (k-1) ) -∆ k-1 ] G k = G k-1 +γ k ∂ 2 θ log p(y, φ (k) ; θ (k-1) )+ ∂ θ log p(y, φ (k) ; θ (k-1) ) t ∂ θ log p(y, φ (k) ; θ (k-1) ) -G k-1 H k = G k -∆ k t ∆ k
Then we note H( θ) the matrix H k at the convergence, and we approximate the observed information matrix of θ by the inverse of -H( θ).

Concerning the log-likelihood of the observations, we follow the user guide of the Monolix software, -we draw φ 1 , ..., φ s under the conditional distribution p(φ|y, θ) using a Metropolis-Hastings algorithm -we estimate the log-likelihood of the observations by log(l s ) where

l s = 1 s s j=1 p(y|φ j , θ) p(φ j ; θ)
p(φ j |y, θ) and p(φ; θ) is the distribution of the missing data φ. Since p(φ|y, θ) is not known in the general case, we estimate this distribution with a Gaussian distribution, with mean E(φ s |y, θ) and variance V ar(φ s |y, θ) estimated using the simulated φ j 's.

3 The SAEM-MCMC algorithm with criteria

The Metropolis-Hastings algorithm

At the kth step of the SAEM-MCMC algorithm and at the (t + 1)th step of the lth chain simulated by the Metropolis-Hastings algorithm. We note

m 1 ∈ N, ρ 1 ∈ R + and ρ 2 ∈ R + . Given φ (k,l,t) i = x (t) : • If t < m 1 , then generate W t ∼ N (A i µ (k-1) , Γ (k-1) ) → the value obtained is noted w t .
Let the acceptance rate:

ρ(x (t) , w t ) = min p((y ij ) j |θ (k-1) ,φ i ) |φ i =w t p((y ij ) j |θ (k-1) ,φ i ) |φ i =x (t)
, 1

• If t ≥ m 1 , then generate ρ ∼ U [ρ 1 ,ρ 2 ] and W t ∼ N (x (t) , ρΓ (k-1) ) the value obtained is noted w t .
The acceptance rate is:

ρ(x (t) , w t ) = min p((y ij ) j |θ (k-1) ,φ i ) |φ i =w t ×p(φ i |θ (k-1) ) |φ i =w t p((y ij ) j |θ (k-1) ,φ i ) |φ i =x (t) ×p(φ i |θ (k-1) ) |φ i =x (t) , 1 • φ (k,l,t+1) i = w t with probability ρ(x (t) , w t ) x (t)
with probability 1 -ρ(x (t) , w t )

At this stage some parameters need to be calibrated:

-the value of m 1 (iteration at which we change instrumental distribution), -the value of itM C (the length of the Markov chains), -the value of burn (the length of the burn-in period within chain),

-the parameters of the second instrumental distribution ρ 1 and ρ 2 , -the value of L (the number of independent chains).

Parameters of Metropolis-Hastings

These parameters are fixed as follows: we run a Markov chain with the Metropolis-Hastings algorithm for one individual taking the initial value θ 0 for the vector of parameter θ and we calibrate the parameters with the following methods. Then these parameters are used to run the Metropolis-Hastings algorithm for all individuals.

• The second instrumental distribution:

We decided to cut the markov chain into two parts. In the first part, we simulate W t with the prior distribution of the random effects. In the second part of the chain, W t is simulated with a Gaussian random walk centered on the precedent iteration of the chain.

According to [START_REF] Robert | Monte Carlo Statistical methods[END_REF], the acceptance rate must not be greater than 50%. So we fix ρ 1 and ρ 2 such that the acceptance rate is between 30% and 40 %.

• m1: The simulated markov chain does not have the same behavior in the first and in the second part of the chain. In general, in the first part, since the new simulation W t can be very far from x (t) , W t is often rejected and in this case the chain stays constant.

In the second part, the chain varies more. So we choose to fix m 1 = itM C/10.

• burn: many ways to choose the burn parameter are presented in detail in [START_REF] Robert | Monte Carlo Statistical methods[END_REF]. We fixed it at itM C/2.

• itM C: is it necessary that the chain reaches the stationary distribution to better estimate the conditional expectations? Indeed in the MCEM algorithm [START_REF] Wei | A Monte Carlo implementation of the EM algorithm and the Poor's Man's data augmentation algorithms[END_REF], the longer the chains are, the more precise are the estimations of the conditional expectations. Nevertheless in the SAEM-MCMC algorithm the Stochastic Approximation step may be a way to exempt the convergence of the chain towards the stationnary distribution. In this study, we prefer to run enough iterations and test the convergence of the chains towards the stationnary distribution thanks to the [START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF] criterion. This criterion provides a diagnosis of the Markov chain convergence by comparing within-chain and between-chain variances.

• L: We know that the more chains are run, the better is the estimation of E S(y, φ

(k,l,itM C) i
)|y, θ .

In practice we can run several chains and compare the behavior of each of them. If their properties (mean, variance) are very different, we may choose L = 5 or L = 10 chains, else, the L = 1 chain can be enough.

3.3 Convergence and "smoothing" criteria 3.3.1 A convergence criteria for the SAEM algorithm [START_REF] Booth | Maximizing generalized linear mixed models likelihoods with an automated Monte Carlo EM algorithm[END_REF] propose to use in the MCEM algorithm the same standard stopping rule as in the deterministic EM algorithm:

e(k) = max j |θ (k) j -θ (k-1) j | |θ (k) j | + δ 1 < δ 2 (1)
where δ 1 and δ 2 have small values and can be δ 1 = 0.001 and δ 2 = 0.0001 (Searle, 1992 p.296). However, [START_REF] Booth | Maximizing generalized linear mixed models likelihoods with an automated Monte Carlo EM algorithm[END_REF] realized that the algorithm can satisfy the criterion thanks to random chance, and so they recommend to stop the algorithm when the criterion has satisfied three consecutive iterations.

On the contrary, [START_REF] Jank | Implementing and diagnosing the stochastic approximation EM algorithm[END_REF] has shown that criterion (1) is not adapted if the sequence (γ k ) k is such that γ k ∝ (1/k) α with α ≈ 1. Indeed in this case the slow convergence of the algorithm implies that the criterion may lead to a bad estimation of the parameters. In this sense he proposes a stopping rule based on an adaptation of the well known property of increase of the sequence (Q k (θ, θ (k-1) )) k in the EM algorithm.

In our method, we begin the Approximation Step at the Kth iteration and the choice of K depends on the variation of the parameters (see Section 3.3.2). Moreover we run enough iterations in the Metropolis-Hastings algorithm to obtain a good estimation of the conditional expectations. In these conditions we think that the algorithm has almost converged at iteration K and finally we just smooth the estimations thanks to the Approximation

Step. Finally we propose to use criterion (1) with δ 1 = 0, δ 2 = 0.0001, and we note e l (i) the lth component of e(i).

A "smoothing" criterion

In order to determine the parameter K, the iteration at which the sequence (γ k ) k is decreasing in the SAEM-MCMC algorithm, we studied a "smoothing" criterion that is based on the variation of the e(i)'s.

Heuristic Figure 1 represents an illustration of the evolution of e 1 (for β a fixed parameter of the model) during the iterations of the SAEM-MCMC algorithm. K was fixed at 43 iterations. The slope of the curve is high at the beginning and then after iteration around 20, e 1 varies slightly around a small positive value. At this moment, the curve can be smoothed.

At each iteration t of the SAEM algorithm, we fit a linear regression on the ten last points (e(s)) s=t-9,...,t and when the slope of the linear regression is not increasing, we add 15 iterations to ensure that we are really in the neighborhood of the maximum likelihood estimator. We define:

slope(l, k) = 10 m=1 m -1 10 10 s=1 s e l (k -m + 1) -1 10 9 s=0 e l (k -s) 10 m=1 (m -1 10 10 s=1 s) 2 K = 15 + min k, ∀l slope(l, k) -slope(l, k -1) < 0

Differences between our SAEM-MCMC and the classical SAEM-MCMC

The two algorithms are based on the same theory presented in [START_REF] Kuhn | Coupling a stochastic approximation version of EM with a MCMC procedure[END_REF]. Nevertheless many parameters need to be calibrated by the user in the classical SAEM-MCMC algorithm implemented for example in Monolix software (Lavielle 2005, Monolix user guide manuel at http://www.monolix.org/). [START_REF] Kuhn | Maximum Likelihood estimation in non linear mixed effects models[END_REF] proposed to do the Approximation Stochastic step on L = 10 chains maximum with length itM C = 10 iterations. In general it is well known that MCMC methods need Markov chains with millions of iterations until they converge towards the stationary distribution [START_REF] Robert | Monte Carlo Statistical methods[END_REF]. About the parameters proposed in [START_REF] Kuhn | Maximum Likelihood estimation in non linear mixed effects models[END_REF], probably 10 iterations within chain are not enough to reach the stationary distribution. So we want to analyse if the difference of iterations can deteriorate the estimation of θ or if the Stochastic Approximation step is a way for the algorithm to converge towards the maximum likelihood estimator of θ despite short Markov chains.

Application and simulation

The SAEM-MCMC algorithm proposed with some criteria was applied to the orange tree data set and to a simulated data set in order to study its properties.

4.1 The Orange tree data set

The model

We applied our SAEM-MCMC algorithm and the classical SAEM-MCMC algorithm with several sets of parameters to the well known orange tree data, studied by Pinheiro andBates (1995a, 1995b). This data set is available for example on S-Plus and R software and consists of seven measurements of the trunk circumference of each of five orange trees. Following Pinheiro andBates (1995a, 1995b) and [START_REF] Kuhn | Maximum Likelihood estimation in non linear mixed effects models[END_REF], we suggest using the following non-linear mixed model:

y ij = φ i 1 + exp - t ij -β 2 β 3 + ε ij , ∀i ∈ {1, ..., 5} ∀j ∈ {1, .., 7}, (2)
where y ij is the jth measurement at age t ij (in days) of the ith tree, the (ε ij ) ij 's are assumed independent identically distributed with distribution N (0, σ 2 ), and the φ i 's are assumed independent identically distributed with distribution N (β 1 , τ 2 ). The φ i 's are independent from the (ε ij ) ij 's. Let θ = (β 1 , β 2 , β 3 , τ 2 , σ 2 ) be the vector of parameters.

The parameters φ i , β 2 , β 3 have a physical interpretation: φ i corresponds to the asymptotic trunk circumference, β 2 represents the age at which the tree attains half of its asymptotic trunk circumference, and β 3 is the growth scale. The parameters β 2 , β 3 are treated as fixed effects.

The aim of this study was to obtain a good estimator of θ by maximum likelihood estimation. Since the model is linear on φ, we can compare our results with the exact EM algorithm results. Details about the minimal sufficient statistic functions used to estimate θ can be obtained in [START_REF] Kuhn | Maximum Likelihood estimation in non linear mixed effects models[END_REF]. We compare our SAEM-MCMC algorithm with the classical SAEM-MCMC algorithm using six different sets of parameters (Table 1). Time of running of the SAEM-MCMC algorithm is noted time in seconds and let itSAEM the number of iterations in the SAEM-MCMC algorithm, fixed with the stopping rule (1).

For example, set 1. is the set of parameters given in [START_REF] Kuhn | Maximum Likelihood estimation in non linear mixed effects models[END_REF]: 10 iterations with the first instrumental distribution, 0 with the second one, i.e. itM C = 10 and m 1 = 11, L = 10 independent chains, K fixed at 100 iterations, the algorithm stopped at iteration 303 (3 s). Parameters (itM C, m 1 , burn, L, K, itSAEM ) in sets 2 to 6 have been chosen to study the behavior of the estimator of θ. In sets 1 and 2, only the first instrumental distribution is used, with longer Markov chains in set 2 than in set 1. In sets 3 and 4, we used the two instrumental distributions, with longer chains in set 4 than in set 3. Set 5 is composed of the same parameters of set 4 but with a burn-in period of 50 iterations for each Markov chain. Set 6 is composed of the same parameters of set 4 with the value of K (iteration at which the sequence (γ k ) k is decreasing) fixed at 50 iterations.

Results and discussion

Our criteria determined the following parameters for the SAEM-MCMC algorithm: itM C = 500, m1 = 50, burn = 250, L = 10, K = 46 and the convergence criterion stopped the SAEM-MCMC algorithm at the 61st iteration (19 s).

For each set (set 1. to set 6. and our set of parameters), the following initial values of θ were used:

β (0) 1 =100, β (0) 2 =650, β (0)
3 =250, τ 2(0) =500, σ 2(0) =10. Since the model is linear in φ, we can apply the EM algorithm and obtain the exact maximum likelihood estimator of θ. Table 2 presents the estimation of θ and the value of the known log-likelihood of the observations at these points for the classical SAEM-MCMC applied with the six sets of parameters presented in Table 1, our SAEM-MCMC algorithm, and the EM algorithm. The standard error (noted σ( θ)) of the estimator of θ were calculated for the EM algorithm and our method. Except for sets 1 and 3, we obtain good estimators of θ.

Comparison between the results of sets 1. and 2. If we use only the first instrumental distribution (sets 1 and 2) to simulate new iterations in short chains in the Metropolis-Hastings algorithm then the sequence (θ (k) ) k does not converge towards the maximum likelihood estimator (set 1.). We obtain this result because the markov chains do not vary enough, and so the estimation of τ 2 is too small. If the chains are longer (set 2), they are more variable and so the variances are estimated better.

Comparison between the results of sets 1. and 3., sets 2. and 4. Adding a second instrumental distribution in the Metropolis-Hastings algorithm (sets 3 and 4), the estimations for all parameters are better but the estimators for the variances are not exactly the same as those obtained with the EM algorithm.

Comparison between the results of sets 1. and 2., sets 3. and 4. The longer are the Markov chains in the Metropolis-Hastings algorithm, the more precise are the estimators.

Comparison between the results of sets 4. and 5. When we do not run the Stochastic Approximation step with all the iterations of the chains simulated in the Metropolis-Hastings algorithm, that is to say when burn > 0, the estimators are closer to the EM estimates.

Comparison between the results of sets 4. and 6. The choice of the parameter K does not seem to be important in this study, we obtained almost the same results with the two sets of parameters.

The results obtained with our method indicate that it is not necessary to run many iterations with the first instrumental distribution, but Gelman and Rubin's criterion implies to run more iterations with the second one. Our "smoothing" criterion fixed K to 46 iterations. Finally we obtained results close to the EM results.

The standard error of the estimator of θ obtained with our method was close to the EM standard error. In this real data set we only have a few individuals, that's why the standard errors of the estimators are very large, in particular: σ(τ 2 )=650.3 for our method. Moreover the log-likelihood of the observations remains stable around the estimator of the maximum likelihood so we cannot have precise estimators. Nevertheless we have seen that set 5. and our method give the best values for the estimators. In order to validate our method we simulated a data set on the same model (2) with more observations.

A simulated data set 4.2.1 The model

The model used for the simulated data set is similar to (2).

y ij = φ i 1 + exp - t ij -β 2 β 3 + ε ij , ∀i ∈ {1, ..., 100} ∀j ∈ {1, .., 15},
The parameters for generating the 100 simulated data sets are the following: θ = (β 1 = 20, β 2 = 70, β 3 = 30, τ 2 = 10, σ 2 = 0.5), the 15 points of observations t ij are the follow: 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 120, 140, 160, 180, 200. In this section we carry out 100 simulated data sets to compare performances of our SAEM-MCMC algorithm, the classical SAEM-MCMC algorithms previously defined by 6 sets of parameters (set 1. to set 6. presented in Table 1) and the EM algorithm. The choice of the parameters of our method is described in the following section.

Choice of the parameters of SAEM-MCMC

All the parameters were fixed with the criteria presented in Section 3 on the first simulated data set. Then the same parameters were used to apply the SAEM-MCMC algorithm to the other data sets.

• The choice of m 1 : Let us remind the two different instrumental distributions: the first one is the Gaussian distribution with mean β 1 and variance τ 2 . The second instrumental distribution is a Gaussian random walk with mean equal to the value of the precedent iteration and with variance equal to ρτ 2 where ρ ∼ U[ρ 1 , ρ 2 ]. When the initial values of the parameters are quite far from the maximum likelihood estimator, the two instrumental distributions do not provide the same behavior to the markov chain simulated in the Metropolis-Hastings algorithm. To analyse that, we simulated a long chain (100,000 iterations) for the individual 1, with parameters β 1 = 10, β 2 =60, β 3 =20, σ 2 =1, τ 2 =3, and m 1 = 50, 000. The choice of ρ 1 et ρ 2 did not have an incidence on the difference of behavior between the two parts of the chain.

Figure 2 shows that in the first 50,000 iterations, with the first instrumental distribution, the markov chain stays during many iterations at the same place on the contrary to the second part of the markov chain (from iteration 50,000 to 100,000) where it moves around a certain mean value. This fact can be explained by the following: the first instrumental distribution simulates larger variables than the second instrumental distribution and so these variables are more rejected in the first case than in the second case. Since the model is linear in φ i , the stationnary distribution of the markov chain is known:

p(φ 1 |Y 1 , θ) = N (u 1 , V )
where u 1 and V are given in [START_REF] Kuhn | Maximum Likelihood estimation in non linear mixed effects models[END_REF]. Figure 2 clearly shows that the distribution of the markov chain from iterations 40,000 to 50,000 is not Gaussian and so the chain has not reached the stationnary distribution before iteration m1.

On the contrary, the histogram of the iterations from 50,000 to 60,000 of the markov chain and the graph of the density of the distribution N (u 1 , V ) is presented in Figure 3.

It clearly shows that the chain has yet reached the stationnary distribution at the 60,000 th iteration.

The lack of variability of the chain with the first instrumental distribution can lead us to use only the second instrumental distribution. Nevertheless the second instrumental distribution imposes to the chain to stay in a neighboorhood of the element obtained at iteration m 1 . So we prefer to fix a small value for m 1 , in particular: m 1 = itM C/10.

• The choice of ρ 1 and ρ 2 : Several values of ρ 1 and ρ 2 were tested and to obtain an average acceptance rate close to 35 %, the values retained are ρ 1 = 1/4, ρ 2 = 1/2.

• The choice of itM C and burn:

As suggested in section 4.1., the burn-in period was fixed at burn = itM C/2 and we fixed itMC at 300 iterations. With these values, we obtained a Gelman and Rubin coefficient equal to 1.01.

• The choice of L: Several independent chains were compared and their behavior was similar, so we decided to fix L at 1.

The other parameters of the SAEM-MCMC algorithm were fixed during the run of the algorithm.

Results and discussion

The same six sets of parameters presented in Section 4.1 were considered to compare our SAEM-MCMC algorithm with the classical SAEM-MCMC algorithm.

Table 3 presents the mean value of the estimators, the mean computing time to run the SAEM-MCMC algorithm in seconds (noted time), the mean of the number of iterations for the estimation of θ (noted itm) and the mean value of K (noted Km), based on the 100 simulated data sets. Table 4 presents the bias and the Mean Squared Error (MSE) of the estimator of θ obtained with the seven SAEM-MCMC algorithms and the EM algorithm.

The trajectories of the SAEM-MCMC algorithm strongly depend on the value of the initial parameter θ 0 , and particularly when itM C is very small. Indeed in this case, if θ 0 is far from θ, the algorithm may not converge towards the maximum likelihood estimator.

Since the aim of the paper was to study the bias of the estimator with different sets of parameters for the SAEM-MCMC algorithm and not to study the divergence or the convergence of this algorithm, θ 0 was drawn with a Gaussian distribution of mean θ equal to (β1 = 20, β 2 = 70, β 3 = 30, τ 2 = 10, σ 2 = 0.5) and standard deviation 0.1θ for each of the 100 simulated data sets.

Concerning the criteria chosen to compare the eight methods, the estimation of the variances of estimators obtained by the classical SAEM-MCMC algorithm may not be positive because of a bad estimation of the conditional expectations (in particular for set 1. and set 2.). In this case, the observed information matrix is not well estimated and the estimation of the matrix is not always definite positive. Consequently we could not compare the different methods with the MQE (Mean Quadratic Error) criterion. Nevertheless the variances of estimates may be a good indicator to analyse and to see if the curve of the complete likelihood is flat close to the maximum likelihood. With our method, the variances of estimators are the following: σ( θ) = (0.32, 0.27, 0.23, 1.43, 0.02), indicating that the complete likelihood is not flat around the maximum likelihood and we can hope to have precise estimators of the maximum likelihood.

Concerning the mean value of the estimators of the fixed effects, we obtained good results for all the methods. Since the markov chains do not vary enough during the first part of the Metropolis-Hastings algorithm (with the first instrumental distribution), the variance τ 2 estimated as the variability within chains, is underestimated in sets 1. and 2. In the same way results from sets 3. and 4. show that it is necessary to run long chains to have good estimators. The comparison between the results of sets 1. and 3. indicates that adding a second instrumental distribution is relevant for the estimation of the variances τ 2 and σ 2 . The mean value of τ 2 (equal to 8.8) and σ 2 (equal to 0.8) for set 2. is better than the ones for set 3. (τ 2 = 7.4 and σ2 = 1.3), which indicates that we must use several instrumental distributions and also run long chains in the Metropolis-Hastings algorithm.

Concerning the choice of the parameter burn, the results from sets 4. and 6. show that we obtain a light improvement with the burn-in period in the estimation of τ 2 and σ 2 . As in the orange tree data set, the results from sets 4. and 6. indicate that the value of K does not have a large incidence since we obtain the same estimators. This may be because we simulate long chains and so the conditional expectations are good estimates as the first iterations of the SAEM-MCMC algorithm. If we had varied the parameter itM C along the SAEM-MCMC algorithm, that is to say for example short chains at the beginning and long chains at the end, we think that the parameter K would have had more impact on the results.

Our method provides the same results as the EM algorithm, which gives the exact value of the maximum likelihood estimator. Set 5. also provides good estimators close to the EM algorithm estimators but when compared to our method, no criteria are used to choose the parameters of simulation. In terms of computing time, our method is the best one.

Concerning the biases and the MSE, set 5. and our method provide the best results close to the ones of the EM algorithm. In particular, our SAEM-MCMC algorithm provides the best results on biases, especially for the variance τ 2 , and the results are quite similar for MSE. Running long chains with two instrumental distributions reduce the MSE and the bias of the estimators for all the parameters.

To conclude, we need to run long chains and use several instrumental distributions in the Metropolis-Hastings algorithm to have a good estimator for θ, and in particular for the variances. The choice of the parameter burn is relevant and the burn-in period is necessary. Concerning the parameter K, we think that it depends on the length of the chain. The longer the chains are, the smaller is the impact of parameter K on the estimation of θ. One advantage of our method is to propose criteria to fix all these parameters, even if some preliminary statistical study must be done before. However some of these criteria may be stringent. For example, the Gelman and Rubin criterion that diagnose the convergence of Markov chains, imposes generally to run long chains.

In the next section we present the case where the model does not belong to the exponential family class, for which [START_REF] Kuhn | Coupling a stochastic approximation version of EM with a MCMC procedure[END_REF] did not prove the convergence of the SAEM-MCMC algorithm. [START_REF] Wang | EM algorithms for nonlinear mixed effects models[END_REF] presents a method of estimation based on the EM algorithm and applied the SAEM-MCMC algorithm on several non-exponential models.

5 The SAEM-MCMC algorithm in a wider class of models [START_REF] Wang | EM algorithms for nonlinear mixed effects models[END_REF] used the SAEM-MCMC algorithm to estimate parameters of Model (2) with the two parameters φ and β 2 considered as random effects. In this case the model does not belong to the exponential family. Although we cannot use the minimal sufficient statistic to estimate the maximum likelihood estimator of θ, the same strategy as the EM algorithm was adopted to estimate the conditional expectations with the Stochastic Approximation step of the SAEM-MCMC algorithm. In this section, our procedure and the exact SAS NLMIXED procedure based on the quadrature method, were compared on the orange trees data set and on a simulated data set.

The Model

The model is defined by:

y ij = φ 1i 1 + exp - t ij -φ 2i β 3 + ε ij , ∀i ∈ {1, ..., N } ∀j ∈ {1, .., N i }, (3)
where y ij is the jth measurement of the ith tree, the (ε ij ) ij 's are assumed independent identically distributed as N (0, σ 2 ), and the φ i = (φ 1i , φ i2 )'s are random variables assuming independent identically distributed with distribution N (β, D) where β = (β 1 , β 2 ) and

D = τ 2 a τ ab τ ab τ 2 b
the covariance matrix of the φ i 's. Let θ = (β 1 , β 2 , β 3 , D, σ 2 ) the vector of parameters.

The estimation of the parameters

The heuristic of this method of estimation is the following: for the two random parameters φ 1 and φ 2 , the mean and the variance are estimated from the simulations of the missing data, and more precisely from the chains of the metropolis-Hastings algorithm. Concerning the other parameters, in the general case, we need to evaluate the gradient and the hessien of the complete log-likelihood and we apply Newton's method to obtain the estimation of the parameters. The method used to estimate the vector of parameters θ is presented in [START_REF] Wang | EM algorithms for nonlinear mixed effects models[END_REF].

Since log p(y, φ; θ) = log p(y|φ; θ) + log p(φ; θ) where p(φ; θ) is the log-likelihood of the missing data, the EM algorithm at iteration t leads to these estimations:

We note ω = (β 3 , σ 2 ).

β(t) 1 = N i=1 E φ 1i |y,θ (t-1) N τ 2(t) a = N i=1 E (φ 1i - β(t) 1 ) 2 |y,θ (t-1) N β(t) 2 = N i=1 E φ 2i |y,θ (t-1) N τ 2(t) b = N i=1 E (φ 2i - β(t) 2 ) 2 |y,θ (t-1) N τ (t) ab = N i=1 E (φ 1i - β(t) 1 )(φ 2i - β(t)
2 )|y,θ (t-1) N ∂ log p(y,φ;θ (t-1) ) ∂ω = ∂ log p(y|φ,θ (t-1) ) ∂ω can not be solved in closed form. To obtain an estimation of ω (t) , we consider a single iteration of Newton's method:

ω (t) = ω (t-1) -H(φ, ω (t-1) ) -1 P (φ, ω (t-1) )
where H(φ, ω (t-1) ) = ∂ 2 log p(y|φ,θ (t-1) ) ∂ω 2 and P (φ, ω (t-1) ) = ∂ log p(y|φ,θ (t-1) ) ∂ω . When H is not definite positive, that is to say when the estimator is not in the neighborhood of the maximum likelihood, we can use a Marquardt algorithm [START_REF] Marquardt | An Algorithm for Least-Squares Estimation of Nonlinear Parameters[END_REF]: we replace H(φ, ω (t-1) ) by (H(φ, ω (t-1) ) + λ t diag(H(φ, ω (t-1) ))), where λ t is a positive real. In this sense we take λ t large and then (H(φ, ω (t-1) ) + λ t diag(H(φ, ω (t-1) ))) -1 P (φ, ω (t-1) ) ∼ 1 λt P (φ, ω (t-1) ): the algorithm is a descent of gradient.

5.3

The NLMIXED procedure (Pinheiro and Bates, 1995b) The NLMIXED procedure of SAS software estimates parameters in the nonlinear mixed models by maximizing an approximation of the likelihood of the observations. This one can be written under an integral form over the random effects. Different integral approximations are available, the two principal ones are adaptive Gaussian quadrature (exact approximation) and a first-order Taylor series approximation (approximation by linearization). In this paper we used the adaptive Gaussian quadrature approximation, because it seems to be the best method in terms of estimation of the value of the parameters (Pinheiro and Bates, 1995b). It is used to approximate integrals thanks to a given kernel by a weighted average of the integrand evaluated at pre-determined abscissas. In this sense, this method can be viewed as a deterministic version of importance sampling method in which the sample and the weights are fixed beforehand. One disadvantage of this method is that it gives accurate results for a large number of abscissas, and so the computing time can be very large.

Results on the Orange tree data set

In this model seven measurements (Ni=7) of the trunk circumference of each of five orange trees (N=5) were made. The following initial values were taken:

β (0) 1 =150, β (0) 2 =600, β (0) 3 =200, τ 2(0) a =500, τ 2(0) b
=200, τ ab (0)=0 and σ 2(0) =10. Our criteria fixed the parameters of the SAEM-MCMC algorithm for this real data set and we obtained: itM C = 300, m 1 = 30, burn = 150, with a Gelman and Rubin criterion at 1.06, ρ 1 = 0, ρ 2 = 0.4 with an average acceptation rate at 25%. It was difficult to increase this rate because either we would have reduced the interval [ρ 1 ; ρ 2 ], which may affect the behavior of the markov chains, either we may have reduced the length of the chain, which would have generated a problem to satisfy the Gelman and Rubin criterion. Then we fixed L = 5. The other criteria of the algorithm that are setting up during the algorithm runnig fixed K = 37, and the stopping rules stopped the algorithm at 172 iterations, the computing time was equal to 28 seconds. The results are presented in Table 5.

Concerning the estimators and the standard errors for the fixed effects, and the variances τ 2 a and σ 2 , we obtained similar results with the two methods. However, the estimators of τ 2 b and τ ab are very different. We noted that in this study, the NLMIXED procedure was unstable and the estimators strongly depended on the initial values of the parameters. In terms of value of the log-likelihood of the observations, the estimator obtained with our method was the best.

Since the estimated value for the variances are very different between the two methods, we study in the next paragraph a simulated data set based on Model(3).

Results on a simulated data set

We simulated 100 data sets based on Model (3) with the following parameters: N = 100, N i = 15, β 1 = 20, β 2 = 70, β 3 = 30, t ij defined in Section 4.2.1, D = 10 -1 -1 40 and σ 2 = 0.5.

After applying our criteria to fix the parameters of simulations of the SAEM-MCMC algorithm for the first simulated data set, we obtained: itM C = 800, m1 = 80, burn = 400, with Gelman and Rubin criterion at 1.04, L = 1, ρ 1 = 0, ρ 2 = 0.4 with an average rate of acceptation at 30%. The mean computing time was 360 seconds, the mean number of iterations was 143 and the mean K was 38.

Results

The values of the estimators obtained with the two methods were similar for all the parameters. These values were close to the exact parameters, except for τ ab , which indicates that the two methods give coherent results. Concerning biases, MSE, σ( θ) and MQE the results were also similar.

On the contrary to the study of the Orange tree data in the same model, here the NLMIXED method was stable for all initial values. Since we obtained the same results as the exact NLMIXED method, our method seems to be adapted and efficient in maximum likelihood estimation in the general class of the non-linear mixed models.

Conclusion

In summary, we propose some criteria to fix the different parameters of the SAEM-MCMC algorithm presented by Kuhn andLavielle (2004, 2005) in maximum likelihood estimation. We show on the orange tree data and on a simulated data set that we need to run long chains in the Metropolis-Hastings algorithm to obtain precise estimates. These chains must be simulated using several instrumental distributions and taking a burn-in period improves the estimator's value.

In this study our method provides estimates similar to the EM algorithm estimates and the computing time was comparable to that of other software used in maximum likelihood estimation.

Moreover we test our algorithm on a wider class of non-linear models, according to the method of estimation presented by [START_REF] Wang | EM algorithms for nonlinear mixed effects models[END_REF]. We obtain the same estimates as the NLMIXED method and close to the exact parameters on a simulated data set. Our algorithm seems to be adapted to the study of several kinds of non-linear models.

In quantitative genetics and animal breeding, heteroscedasticity has generated much interest. In fact the assumption of homogeneous variances may not always be appropriate.

In linear mixed models, there is now a large amount of experimental evidence of hetero-geneous variances for most important livestock production traits [START_REF] Garrick | Variance heterogeneity in direct and maternal weight by sex and percent purebred for Simmental-sired calves[END_REF][START_REF] Visscher | Estimation of genetic and environmental variances for fat yield in individual herds and an investigation into heterogeneity of variance between herds[END_REF][START_REF] Visscher | Heterogeneity of variance and dairy cattle breeding[END_REF][START_REF] Robert-Granié | Accounting for variance heterogeneity in french dairy cattle genetic evaluation[END_REF]. Major theoretical and applied work has been carried out for estimating and testing the source of heterogeneous variances arising in univariate linear mixed models [START_REF] Foulley | A method for assessing extent and sources of heterogeneity of residual vriances in mixed linear models[END_REF][START_REF] Gianola | Estimation of heterogeneous variances using empirical Bayes methods: theoretical considerations[END_REF][START_REF] San Cristobal | Inference about multiplicative heteroscedastic components of variance in a mixed linear Gaussian model with an application to beef cattle breeding[END_REF][START_REF] Foulley | Heterogeneous variances in Gaussian linear mixed models[END_REF][START_REF] Foulley | A link function approach to heterogeneous variance components[END_REF]Robert et al., 1995a;Robert et al., 1995b). [START_REF] Foulley | A method for assessing extent and sources of heterogeneity of residual vriances in mixed linear models[END_REF] formalized the heterogeneous variances using a mixed model approach on log variances in linear mixed models. According to this modelisation, it would be interesting to adapt our algorithm to the study of heterogeneous variances in non-linear mixed models.

Figure 1 :

 1 Figure 1: Illustration of the evolution of e 1 during the iterations of the SAEM-MCMC algorithm.

Figure 2 :Figure 3 :

 23 Figure 2: Simulation of a markov chain of the Metropolis-Hastings algorithm.

Table 1 :

 1 ). Presentation of the six sets of parameters used in the SAEM-MCMC algorithm, and the set of parameters fixed thanks to the criteria.

	Set number	itM C m 1 burn L K itSAEM time
	set 1.	10	11	0	10 100	303	3
	set 2.	100	101	0	10 100	195	18
	set 3.	20	11	0	10 100	222	5
	set 4.	200	101	0	10 100	120	22
	set 5.	200	101	50	10 100	137	20
	set 6.	200	101	0	10 50	102	18
	SAEM-MCMC with criteria	500	51	250 10 46	61	19

LogL = log p(y| θ)

Table 2 :

 2 Values of θ for the EM algorithm, the classical SAEM-MCMC algorithm with six sets of parameters and our SAEM-MCMC algorithm on the Orange tree data.The biases and MSE of the estimates, based on 100 simulations.

	θ	β 1	β 2	β 3	τ 2	σ 2
	Bias
	EM	-0.61 1.46 -1.31 0.03 0.09
	set 1.	-0.75 0.74 -1.84 -5.67 1.64
	set 2.	-0.66 1.22 -1.48 -1.23 0.32
	set 3.	-0.66 1.18 -1.51 -2.60 0.85
	set 4.	-0.63 1.35 -1.39 -0.54 0.20
	set 5.	-0.62 1.40 -1.35 -0.14 0.10
	set 6.	-0.63 1.35 -1.39 -0.54 0.20
	SAEM-MCMC with criteria -0.62 1.44 -1.33 0.03 0.09
	MSE
	EM	0.46 2.21 1.77 2.35 0.01
	set 1.	0.75 1.95 4.00 32.84 2.78
	set 2.	0.53 1.75 2.32 3.33 0.11
	set 3.	0.55 1.81 2.48 8.04 0.74
	set 4.	0.49 1.96 1.99 2.36 0.04
	set 5.	0.48 2.08 1.88 2.28 0.01
	set 6.	0.49 1.96 1.99 2.36 0.04
	SAEM-MCMC with criteria 0.47 2.17 1.80 2.35 0.01

Table 4 :

 4 Estimation of the MSE and MQE of the estimators on the simulation data set using Model (2).LogL( θ) is an estimation of the log-likelihood of the observations to the point θ

	θ	β 1	β 2	β 3	τ 2 a	τ 2 b	τ ab	σ 2 log L( θ)
	initial value θ 0	150 600 200	500	200	0	10	-
	SAEM-MCMC with criteria 191 714 344 1169 984	877	57	-130.89
	σ( θ)	16.2 31.3 23.3 761.7 1895 951	16	-
	NLMIXED	192 725 348 1176 193	313	59	-131.2
	σ( θ)	16.7 37.4 26.7 905.3 2180 838.3 18.9	-

Table 5 :

 5 Value of θ and the standard error of θ for our SAEM-MCMC algorithm and the NLMIXED procedure on the Orange tree data.LogL( θ) is an estimation of the log-likelihooh of the observations to the point θ

	θ	β 1	β 2	β 3	τ 2 a	τ 2 b	τ ab	σ 2	log L( θ)
	exact	20	70	30	10	40	-1	0.5	-
	initial values θ 0	10	35	15	5	30	0	0.1	-
	SAEM-MCMC with criteria 19.41 71.51 28.67 9.77 40.17 -2.09 0.59 -2071.25
	biais	-0.59 1.51 -1.33 -0.23 0.17 -1.09 0.09	-
	MSE	0.43 2.66 1.81 2.48 36.94 5.06 0.01	-
	σ( θ)	0.32 0.69 0.23 1.41 6.52 2.16 0.02	-
	MQE	0.53 3.15 1.86 4.50 80.29 9.80 0.01	-
	NLMIXED	19.41 71.49 28.68 9.81 40.15 -2.12 0.59 -2071.00
	biais	-0.59 1.49 -1.32 -0.19 0.15 -1.12 0.09	-
	MSE	0.43 2.59 1.80 2.36 35.77 5.20 0.01	-
	σ( θ)	0.32 0.69 0.23 1.41 6.52 2.16 0.02	-
	MQE	0.53 3.08 1.86 4.40 78.97 9.94 0.01	-

Table 6 :

 6 Value of estimates, biases, MSE and MQE of θ for our SAEM-MCMC procedure and the NLMIXED procedure on a simulated data.

The mean estimates are based on 100 simulations.

Table 3: Estimation of the mean estimates for the estimators on the simulation data set using Model (2).