Maud Delattre

Estelle Kuhn

Computing an empirical Fisher information matrix estimate in latent variable models through stochastic approximation

Keywords: Model-based standard error, moment estimate, Fisher identity, stochastic approximation algorithm

The Fisher information matrix (FIM) is a key quantity in statistics. However its exact computation is often not trivial. In particular in many latent variable models, it is intricated due to the presence of unobserved variables. Several methods have been proposed to approximate the FIM when it can not be evaluated analytically. Different estimates have been considered, in particular moment estimates. However some of them require to compute second derivatives of the complete data log-likelihood which leads to some disadvantages. In this paper, we focus on the empirical Fisher information matrix defined as an empirical estimate of the covariance matrix of the score, which only requires to compute the first derivatives of the log-likelihood. Our contribution consists in presenting a new numerical method to evaluate this empirical Fisher information matrix in latent variable model when the proposed estimate can not be directly analytically evaluated. We propose a stochastic approximation estimation algorithm to compute this estimate as a by-product of the parameter estimate. We evaluate the finite sample size properties of the proposed estimate and the convergence properties of the estimation algorithm through simulation studies.

Introduction

The Fisher information matrix (FIM) is a key quantity in statistics as it is required for examples for evaluating asymptotic precisions of parameter estimates, for building optimality criteria in experimental designs, for computing Wald test statistics or classical asymptotic distributions in statistical testing [START_REF] Van Der Vaart | Asymptotic Statistics[END_REF]. It also appears more recently in post model selection inference [START_REF] Charkhi | Asymptotic Post-Selection Inference for the Akaike Information Criterion[END_REF], in asymptotic distribution of the likelihood ratio test statistics when testing variance component in mixed models (Baey C., Cournède P.-H., and Kuhn E. 2019) or as a particular Riemannian metric on complex manifold (Le [START_REF] Brigant | Fisher-Rao Geometry of Dirichlet Distributions[END_REF]. However its exact computation is often not trivial. This is in particular the case in many latent variables models, also called incomplete data models, due to the presence of the unobserved variables. Though these models are increasingly used in many fields of application, such as in ecophysiology [START_REF] Technow | Integrating Crop Growth Models with Whole Genome Prediction Through Approximate Bayesian Computation[END_REF], in genomic [START_REF] Picard | A Segmentation/Clustering Model for the Analysis of Array CGH Data[END_REF] or in ecology [START_REF] Gloaguen | An Autoregressive Model to Describe Fishing Vessel Movement and Activity[END_REF]. They especially allow a better consideration of the different variability sources and when appropriate, a more precise characterization of the known mechanisms at the origin of the data. When the FIM can not be exactly computed, people either approximate it numerically, for example by using Monte Carlo technics like developed in the R package MIXFIM (Riviere-Jourdan M.-K. and Mentre F. 2018) or focus on an estimate of the FIM. The probably most widely used is the observed FIM [START_REF] Efron | Assessing the Accuracy of the Maximum Likelihood Estimator: Observed Versus Expected Fisher Information[END_REF]. When it can not be directly computed in latent variable models, several methods have been proposed to approximate it. Among the most frequently used approaches are Monte-Carlo methods or iterative algorithms derived from the missing information principle [START_REF] Orchard | A Missing Information Principle: Theory and Applications[END_REF]. Indeed according to this principle, the observed Fisher information matrix can be expressed as the difference between two matrices corresponding to the complete information and the missing information due to the unobserved variables (see e.g. [START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF]) chapter 4). It enables the development of alternative methods to compute the observed FIM: the Louis's method [START_REF] Louis | Finding the Observed Information Matrix When Using the EM Algorithm[END_REF], combined with a Monte Carlo method or a stochastic approximation algorithm by [START_REF] Delyon | Convergence of a Stochastic Approximation Version of the EM Algorithm[END_REF]), the Oakes method [START_REF] Oakes | Direct Calculation of the Information Matrix via the EM Algorithm[END_REF] or the supplemented Expectation Maximization algorithm [START_REF] Meng | Using EM to Obtain Asymptotic Variance-Covariance Matrices: The SEM Algorithm[END_REF]. However as the observed FIM involves the second derivatives of the observed log-likelihood, all these methods require to compute second derivatives of the complete data log-likelihood which leads to some disadvantages from a computational point of view. More recently, [START_REF] Meng | Efficient Computation of the Fisher Information Matrix in the EM Algorithm[END_REF] proposed an accelerated algorithm based on numerical first order derivatives of the conditional expectation of the log-likelihood. Another estimate is the empirical Fisher information matrix. This estimator of the FIM is defined as the moment estimate of the covariance matrix of the score. It is much less used than the observed Fisher information matrix. However it has a nice property since it is positive definite, which is not systematically the case for the latter and it is numerically more interesting because it only requires the calculation of the first derivatives of the log-likelihood.

In this paper, our contribution consists in presenting a new numerical method to evaluate the empirical FIM in latent variables model. Indeed, when the proposed estimate can not be directly analytically evaluated, we propose a stochastic approximation estimation algorithm to compute it, which provides this estimate of the FIM as a by-product of model parameter estimates.

The paper is organized as follows. In Section 2, we recall the three main FIM estimates and discuss their immediate properties. In Section 3, we give practical tools for the computation of the empirical Fisher information matrix in incomplete data models. In particular, we introduce a new stochastic approximation procedure based on the first derivatives of the complete log-likelihood only and state its asymptotic properties. In Section 4, we illustrate the finite sample size properties of both estimators and the convergence properties of the computation algorithm through simulations. The paper ends by a discussion.

Moment estimates of the Fisher information matrix

Let us consider a random vector 𝑌 taking value in 𝒴. Assume 𝑌 admits a density 𝑔(⋅; 𝜃) with respect to a given common measure 𝜇, depending on some parameter 𝜃 taking values in an open subset Θ of ℝ 𝑑 , such that the log-likelihood function log 𝑔 is differentiable on Θ and ‖𝜕 𝜃 log 𝑔(𝑦; 𝜃)(𝜕 𝜃 log 𝑔(𝑦; 𝜃)) 𝑡 ‖ is integrable with respect to 𝑔, where 𝑥 𝑡 stands for the transpose of a vector or a matrix 𝑥. Then, by definition (see [START_REF] Lehmann | Theory of Point Estimation[END_REF])), the Fisher information matrix is given for all 𝜃 ∈ Θ by: 𝐼 (𝜃) = 𝐸 𝜃 [𝜕 𝜃 log 𝑔(𝑌 ; 𝜃)(𝜕 𝜃 log 𝑔(𝑌 ; 𝜃)) 𝑡 ] .

(1)

When this expression can not be analytically evaluated, people are interested in computing an estimate of the Fisher information matrix. Considering this expression, one can derive a first moment estimator of the Fisher information matrix based on a 𝑛-sample 𝑦 = (𝑦 1 , … , 𝑦 𝑛 ) of independent observations:

𝐼 𝑛,𝑠𝑐𝑜 (𝜃, 𝑦) = 1 𝑛 𝑛 ∑ 𝑖=1 𝐼 𝑠𝑐𝑜 (𝜃, 𝑦 𝑖 ) = 1 𝑛 𝑛 ∑ 𝑖=1
𝜕 𝜃 log 𝑔(𝑦 𝑖 ; 𝜃)(𝜕 𝜃 log 𝑔(𝑦 𝑖 ; 𝜃)) 𝑡 .

This estimate is indeed equal to the mean of the Gram matrices of the scores.

Moreover, we can get another expression for the Fisher information (see [START_REF] Lehmann | Theory of Point Estimation[END_REF]). If we assume that the set 𝐴 = {𝑦, 𝑔(𝑦; 𝜃) > 0} is independent of 𝜃, that for 𝜇-almost all 𝑦, 𝑔(𝑦; ⋅) is differentiable on Θ, and that the derivative with respect to 𝜃 on the left side of ∫ 𝑔(𝑦; 𝜃)𝑑𝜇(𝑦) = 1

(2) can be obtained by differentiating under the integral sign, then the Fisher information matrix is given for all 𝜃 ∈ Θ by: 𝐼 (𝜃) = 𝑉 𝜃 [𝜕 𝜃 log 𝑔(𝑌 ; 𝜃)] .

(3)

One can also derive a second estimate from this expression defined as

𝐼 𝑛,𝑐𝑜𝑣 (𝜃, 𝑦) = 1 𝑛 𝑛 ∑ 𝑖=1
𝜕 𝜃 log 𝑔(𝑦 𝑖 ; 𝜃)(𝜕 𝜃 log 𝑔(𝑦 𝑖 ; 𝜃)) 𝑡 -s (𝜃, 𝑦) s (𝜃, 𝑦) 𝑡 , where s (𝜃, 𝑦) = 1 𝑛 ∑ 𝑛 𝑖=1 𝜕 𝜃 log 𝑔(𝑦 𝑖 ; 𝜃) (see e.g. [START_REF] Scott | Maximum Likelihood Estimation Using the Empirical Fisher Information Matrix[END_REF]). We emphasize here that the terminology "empirical Fisher information matrix" is used in the literature for both estimates (see e.g. [START_REF] Kunstner | Limitations of the Empirical Fisher Approximation for Natural Gradient Descent[END_REF]).

Moreover if additionally the second derivative with respect to 𝜃 of log 𝑔(𝑦; 𝜃) exists for all 𝑦 and 𝜃 and the second derivative with respect to 𝜃 of the left side of Equation 2 can be obtained by differentiating twice under the integral sign (see [START_REF] Lehmann | Theory of Point Estimation[END_REF])), we have

𝐼 (𝜃) = -𝐸 𝜃 [𝜕 2 𝜃 log 𝑔(𝑌 ; 𝜃)] . (4)
Considering this third expression, we can derive another moment estimator of the Fisher information matrix based on a 𝑛-sample (𝑦 1 , … , 𝑦 𝑛 ) of observations, called the observed Fisher information matrix defined as:

𝐼 𝑛,𝑜𝑏𝑠 (𝜃, 𝑦) = 1 𝑛 𝑛 ∑ 𝑖=1 𝐼 𝑜𝑏𝑠 (𝜃, 𝑦 𝑖 ) = - 1 𝑛 𝑛 ∑ 𝑖=1 𝜕 2 𝜃 log 𝑔(𝑦 𝑖 ; 𝜃).
Some detailed discussion about the three estimators above can be found in [START_REF] Scott | Maximum Likelihood Estimation Using the Empirical Fisher Information Matrix[END_REF].

Remark. We emphasize that the estimate 𝐼 𝑛,𝑠𝑐𝑜 (𝜃, 𝑦) is always positive semi-definite, since it is a mean of Gram matrices, contrary to the other estimates 𝐼 𝑛,𝑜𝑏𝑠 (𝜃, 𝑦) and 𝐼 𝑛,𝑐𝑜𝑣 (𝜃, 𝑦). Moreover assuming 𝑛 sufficiently large allows to prove positive definiteness of 𝐼 𝑛,𝑠𝑐𝑜 (𝜃, 𝑦). Consider for any nonzero vector 𝑥 the quantity 𝑥 𝑡 𝐼 𝑛,𝑠𝑐𝑜 (𝜃, 𝑦)𝑥. We have that 𝑥 𝑡 𝐼 𝑛,𝑠𝑐𝑜 (𝜃, 𝑦)𝑥 = (∑ 𝑛 𝑖=1 𝑥 𝑡 𝜕 𝜃 log 𝑔(𝑦 𝑖 ; 𝜃)𝜕 𝜃 log 𝑔(𝑦 𝑖 ; 𝜃) 𝑡 𝑥)/𝑛 = ∑ 𝑛 𝑖=1 (𝑥 𝑡 𝜕 𝜃 log 𝑔(𝑦 𝑖 ; 𝜃)) 2 /𝑛. Thus, 𝑥 𝑡 𝐼 𝑛,𝑠𝑐𝑜 (𝜃, 𝑦)𝑥 = 0 implies that 𝑥 𝑡 𝜕 𝜃 log 𝑔(𝑦 𝑖 ; 𝜃) = 0 for all 1 ≤ 𝑖 ≤ 𝑛. If 𝑛 is sufficiently large, there exist 𝑑 indexes 𝑖 1 , ..., 𝑖 𝑑 such that the family of vectors {𝜕 𝜃 log 𝑔(𝑦 𝑖 𝑙 ; 𝜃), 1 ≤ 𝑙 ≤ 𝑑} is linearly independent. Thus this implies that 𝑥 = 0 leading to the results.

Remark. The asymptotical properties of the estimates 𝐼 𝑛,𝑠𝑐𝑜 (𝜃, 𝑦) and 𝐼 𝑛,𝑜𝑏𝑠 (𝜃, 𝑦) are straighforward when considering independent and identically distributed sample (𝑦 1 , … , 𝑦 𝑛 ). In particular, assuming standard regularity conditions on 𝑔, it follows directly from the central limit theorem that 𝐼 𝑛,𝑠𝑐𝑜 (𝜃, 𝑦) and 𝐼 𝑛,𝑜𝑏𝑠 (𝜃, 𝑦) are asymptotically normal. If the variables 𝑌 1 , … , 𝑌 𝑛 are independent not identically distributed, for example if their distributions depend on some individual covariates which is often the case in practice, we can also get asymptotic properties for the estimates assuming more strengthed reguarity conditions by applying for example the Kolmogorov criterion (see e.g. [START_REF] Feller | An Introduction to Probability Theory and Its Applications[END_REF])) for the consistency and the Lindeberg theorem for the normality result (see theorem 27.2 of [START_REF] Billingsley | Convergence of Probability Measures[END_REF])).

Remark. Since both estimators 𝐼 𝑛,𝑠𝑐𝑜 (𝜃, 𝑦) and 𝐼 𝑛,𝑜𝑏𝑠 (𝜃, 𝑦) are moment estimates of 𝐼 (𝜃), they are unbiased for all 𝜃 ∈ Θ. This is not the case for 𝐼 𝑛,𝑐𝑜𝑣 (𝜃, 𝑦). Regarding the variance, none of both estimators is better than the other one. This can be highlighted through the following examples. First consider a Gaussian sample with unknown expectation and fixed variance. Then, the variance of the estimator 𝐼 𝑛,𝑜𝑏𝑠 (𝜃, 𝑦) is zero whereas the variance of the estimator 𝐼 𝑛,𝑠𝑐𝑜 (𝜃, 𝑦) is positive. Second consider a centered Gaussian sample with unknown variance. Then, the variance of 𝐼 𝑛,𝑠𝑐𝑜 (𝜃, 𝑦) is smaller than the variance of 𝐼 𝑛,𝑜𝑏𝑠 (𝜃, 𝑦). Therefore, none of both estimators is more suitable than the other in general from this point of view.

3 Computing the estimator 𝐼 𝑛,𝑠𝑐𝑜 (𝜃) in latent variable model Let us consider independent random variables 𝑌 1 , … , 𝑌 𝑛 . Assume in the sequel that there exist independent random variables 𝑍 1 , … , 𝑍 𝑛 taking values in 𝒵 and a measure 𝜆 on 𝒵 such that for each 1 ≤ 𝑖 ≤ 𝑛, the random vector (𝑌 𝑖 , 𝑍 𝑖 ) admits a parametric probability density function denoted by 𝑓 parametrized by 𝜃 ∈ Θ with respect to 𝜇 × 𝜆 on 𝒴 × 𝒵. We present in this section dedicated tools to compute the estimator 𝐼 𝑛,𝑠𝑐𝑜 (𝜃) in latent variable model when it can not be evaluated analytically.

Analytical expressions in latent variable models

In latent variable models, the estimator 𝐼 𝑛,𝑠𝑐𝑜 (𝜃, 𝑦) can be expressed using the conditional expectation as stated in the following proposition.

Proposition 3.1. Assume that for all 𝑦 and all 𝜃 ∈ Θ the function 𝑓 (𝑦, .; 𝜃) is integrable with respect to 𝜆, that for all 𝑦 and for 𝜆-almost all 𝑧 the function 𝑓 (𝑦, 𝑧; ⋅) is differentiable on Θ, that there exists a mesurable function 𝑚 such that ∫ 𝑚(𝑧)𝜆(𝑑𝑧) < ∞ and for all 𝜃 ∈ Θ and for 𝜆-almost all 𝑧 |𝜕 𝜃 𝑓 (𝑦, 𝑧; 𝜃)| ≤ 𝑚(𝑧). Then for all 𝜃 ∈ Θ and all 𝑛 ∈ ℕ * :

𝐼 𝑛,𝑠𝑐𝑜 (𝜃) = 1 𝑛 𝑛 ∑ 𝑖=1 E 𝑍 𝑖 |𝑌 𝑖 ;𝜃 (𝜕 𝜃 log 𝑓 (𝑌 𝑖 , 𝑍 𝑖 ; 𝜃))E 𝑍 𝑖 |𝑌 𝑖 ;𝜃 (𝜕 𝜃 log 𝑓 (𝑌 𝑖 , 𝑍 𝑖 ; 𝜃)) 𝑡 , ( 5
)
where E 𝑍 |𝑌 ;𝜃 denotes the expectation under the law of 𝑍 conditionally to 𝑌.

We apply the classical Fisher identity [START_REF] Fisher | Statistical Methods for Research Workers[END_REF] to establish the equality stated in Proposition 3.1. We refer to Proposition 100 of [START_REF] Cappé | Inference in Hidden Markov Models[END_REF] for the statement of the Fisher identity. This statement is indeed in the same spirit as the well-known Louis formulae for the observed Fisher information matrix estimate [START_REF] Louis | Finding the Observed Information Matrix When Using the EM Algorithm[END_REF]. The result follows directly.

Remark. In some specific cases the conditional expectations involved in the previous proposition admit exact analytical expressions for example in mixture models which are developed in Section 4 in some simulation studies.

Computing 𝐼 𝑛,𝑠𝑐𝑜 (𝜃) using stochastic approximation algorithm

When exact computation of the estimator 𝐼 𝑛,𝑠𝑐𝑜 (𝜃) is not possible for all 𝜃 ∈ Θ, we propose to evaluate its value by using a new stochastic algorithm which provides the estimate 𝐼 𝑛,𝑠𝑐𝑜 ( θ𝑀𝐿 ) as a by-product of the maximum likelihood estimate θ𝑀𝐿 . More precisely we provide three algorithms: a first one in the curved exponential family context which requires to simulate the latent variable from a transition kernel of an ergodic Markov chain and assumes less strong assumptions to get theoretical convergence result thanks to a truncation on random boundaries step, a second one in the curved exponential family context which does not include this additional projection step but requires stronger assumptions to ensure theoretical convergence. This second one and the related results are presented in Appendix. Finally we provide a third algorithm dedicated to general latent variables models without any theoretical results as it is usually the case for such kind of methods (see Section 3.2.3).

Description of the algorithm with truncation on random boundaries in curved exponential family model

We develop an extension of the stochastic approximation Expectation Maximization algorithm coupled with a Monte Carlo Markov Chain studied by (Allassonnière S., Kuhn E., and Trouvé A. 2010) which allows to compute simultaneously the maximum likelihood estimate and the FIM estimate proposed in the previous section. We assume in this section that all the individual complete loglikelihoods belong to the curved exponential family (see [START_REF] Bickel | Mathematical Statistics: Basic Ideas and Selected Topics, Volumes i-II Package[END_REF]) for stating the theoretical results. As our estimate involves individual conditional expectations, we have to consider an extended form of sufficient statistics for the model at the individual level. Indeed, it is necessary to compute stochastic approximation of each individual sufficient statistic at individual level since they are required to be able to compute the proposed FIM estimate. This is the main difference with the usual algorithm. Therefore we introduce the following notations and assumptions.

The individual complete data likelihood function is given for all 1 ≤ 𝑖 ≤ 𝑛 by:

𝑓 𝑖 (𝑧 𝑖 ; 𝜃) = exp (-𝜓 𝑖 (𝜃) + ⟨𝑆 𝑖 (𝑧 𝑖 ), 𝜙 𝑖 (𝜃)⟩) , (6)
where ⟨⋅, ⋅⟩ denotes the scalar product, 𝑆 𝑖 is a function on Moreover we assume that there exists a function θ ∶ 𝒮 → Θ, such that ∀𝑠 ∈ 𝒮 , ∀𝜃 ∈ Θ, 𝐿(𝑠; θ (𝑠)) ≥ 𝐿(𝑠; 𝜃).

Initialization step:

Initialize arbitrarily for all 1 ≤ 𝑖 ≤ 𝑛 𝑠 0 𝑖 and 𝜃 0 . Set 𝜅 0 = 𝜁 0 = 𝜈 0 = 0. Repeat until convergence the three steps defined at iteration 𝑘 by:

• Simulation step: for 1 ≤ 𝑖 ≤ 𝑛 simulate a realization Z 𝑖 from a parametric transition kernel Π 𝑖 of a Markov Chain parametrized by the current parameter value 𝜃 𝑘-1 and having the conditional distribution given the observations 𝑌 𝑖 denoted by 𝑝 𝑖 as stationary distribution

• Stochastic approximation step: compute the quantities for all 1 ≤ 𝑖 ≤ 𝑛

s𝑖 = (1 -𝛾 𝑘 )𝑠 𝑘-1 𝑖 + 𝛾 𝑘 𝑆 𝑖 (𝑍 𝑘 𝑖 )
where (𝛾 𝑘 ) is a sequence of positive step sizes satisfying ∑ 𝛾 𝑘 = ∞ and ∑ 𝛾 2 𝑘 < ∞.

• Truncation step: Let us denote Z = ( Z 𝑖 ), s = ( s𝑖 ) and 𝑠 = (𝑠 𝑖 ).

If s ∈ 𝒦 𝜅 𝑘-1 and ‖ s -𝑠 𝑘-1 ‖ ≤ 𝜀 𝜁 𝑘-1 , then set (𝑍 𝑘 , 𝑠 𝑘 ) = ( Z , s ), 𝜅 𝑘 = 𝜅 𝑘-1 , 𝜈 𝑘 = 𝜈 𝑘-1 +1, 𝜁 𝑘 = 𝜁 𝑘-1 +1, else set (𝑍 𝑘 , 𝑠 𝑘 ) = ( Z , s ) ∈ K×𝒦 0 , 𝜅 𝑘 = 𝜅 𝑘-1 + 1, 𝜈 𝑘 = 0, 𝜁 𝑘 = 𝜁 𝑘-1 + Ψ(𝜈 𝑘-1 )
where Ψ ∶ ℕ → ℤ is a function such that Ψ(𝑘) > 𝑘 for any 𝑘 and ( Z , s ) chosen arbitrarily.

• Maximisation step: update of the parameter estimator according to:

𝜃 𝑘 = arg max 𝜃 𝑛 ∑ 𝑖=1 (-𝜓 𝑖 (𝜃) + ⟨𝑠 𝑘 𝑖 , 𝜙 𝑖 (𝜃)⟩) = θ (𝑠 𝑘 )
When convergence is reached, say at iteration 𝐾 of the algorithm, evaluate the FIM estimator according to:

𝐼 𝐾 𝑛,𝑠𝑐𝑜 = 1 𝑛 𝑛 ∑ 𝑖=1 Δ𝑖 (𝑠 𝐾 ) Δ𝑖 (𝑠 𝐾 ) 𝑡
where Δ𝑖 (𝑠) = -𝜕𝜓 𝑖 ( θ (𝑠)) + ⟨𝑠 𝑖 , 𝜕𝜙 𝑖 ( θ (𝑠))⟩ for all 𝑠.

Remark. Note that the projection step which is done through the truncation procedure on random boundaries ensures the stability of the algorithm in particular for the theoretical analysis provided below. More details on this projection step are available in (Andrieu C., Moulines E., and Priouret P. 2005). However, the authors do not provide recommendations on the choice of Z and s . In practice, one could use ( Z , s ) = (𝑍 𝑘-1 , 𝑠 𝑘-1 ) or ( Z , s ) = (𝑍 0 , 𝑠 0 ), as done in the numerical experiments in Section 4.

Theoretical convergence property

The theoretical result provided in this section for the sequence (𝜃 𝑘 ) generated by the algorithm with truncation on random boundaries is based on that of (Allassonnière S., Kuhn E., and Trouvé A. 2010). Indeed it established convergence guarantees for the FIM estimate obtained as a by-product of that for the MLE. To that purpose, in addition to the exponential family assumption for each individual likelihood, we also make the same type of regularity assumptions as those presented in (Allassonnière S., Kuhn E., and Trouvé A. 2010) at each individual level. These assumptions are detailed in the appendix section.

We establish our theoretical result for transition kernels (Π 𝑖 ) corresponding to those of the random walk Metropolis Hastings algorithm [START_REF] Jarner | Geometric Ergodicity of Metropolis Algorithms[END_REF]. We denote by (𝑞 𝑖 ) the family of symmetric densities used to generate the candidate with the proposal distribution. We introduce additional assumptions required to control the stochastic behavior of the algorithm:

• (H1) There exists a constant 𝑀 0 such that

ℒ = {𝑠 ∈ 𝒮 , ⟨∇𝑙( θ (𝑠)), ℎ(𝑠)⟩ = 0} ⊂ {𝑠 ∈ 𝒮 , -𝑙( θ (𝑠)) < 𝑀 0 }.
In addition, there exists 𝑀 1 ∈ (𝑀 0 , ∞] such that {𝑠 ∈ 𝒮 , -𝑙( θ (𝑠)) ≤ 𝑀 1 } is a compact set.

• (H2) For all 𝑠 ∈ 𝒮, lim 𝑧→∞ 𝑛(𝑧).∇ 𝑧 log 𝑝(𝑧; θ (𝑠)) = -∞ and lim 𝑧→∞ sup 𝑛(𝑧).𝑚 𝑠 (𝑧) < 0 where where 𝑛(𝑧) = 𝑧/|𝑧| for 𝑧 ≠ 0, and 𝑚 𝑠 (𝑧) = ∇ 𝑧 𝑝(𝑧; θ (𝑠))/𝑝(𝑧; θ (𝑠)) with 𝑝(𝑧; 𝜃) = ∏ 𝑖 𝑝 𝑖 (𝑧 𝑖 ; 𝜃).

• (H3)

The family {𝑞 𝑖 } 1≤𝑖≤𝑛 of symmetric densities is such that, for 𝑖 = 1, … , 𝑛, there exist constants 𝜂 𝑖 > 0 and 𝛿 𝑖 < ∞ such that 𝑞 𝑖 (𝑧) > 𝜂 𝑖 for all |𝑧| < 𝛿 𝑖 .

• (H4) There exist 𝐶 > 1, 𝜌 ∈ (0, 1) and 𝜃 0 ∈ Θ such that, for all 𝑧 ∈ ℝ 𝑑 , |𝑆(𝑧)| ≤ 𝐶𝑝(𝑧; 𝜃 0 ) -𝜌 .

Remark. Assumption (𝐻 2) is standard and usually called super-exponentiality property in the literature [START_REF] Jarner | Geometric Ergodicity of Metropolis Algorithms[END_REF].

Remark. We established our results for transition kernels corresponding to random walk Metropolis Hastings algorithms which are of common use in practice. We emphasize that our result can be generalised to more general transition kernels by replacing our assumptions (𝐻 2) and (𝐻 3) by assumption (𝐷𝑅𝐼 ) of (Andrieu C., Moulines E., and Priouret P. 2005) which is more generic. The latter can be verified in practice for more general transition kernels.

Theorem 3.1. Assume that (𝑀1 ′ ) and (𝑀2 ′ ), (𝑀3) to (𝑀5), (𝑆𝐴𝐸𝑀1) and (𝑆𝐴𝐸𝑀2), (𝐻 1) to (𝐻 4) are fulfilled. Let us define ℒ = {𝜃 ∈ Θ, 𝜕 𝜃 𝑙(𝑦; 𝜃) = 0} the set of stationary points of the observed log-likelihood 𝑙 defined as 𝑙(𝑦; 𝜃) = ∑ 𝑛 𝑖=1 log 𝑔(𝑦 𝑖 ; 𝜃). Then, for all 𝜃 0 ∈ Θ, for fixed 𝑛 ∈ ℕ * , we get: lim 𝑘 𝑑(𝜃 𝑘 , ℒ ) = 0 a.s. and lim 𝑘 𝑑(𝐼 𝑘 𝑛,𝑠𝑐𝑜 , ℐ ) = 0 a.s. where ℐ = {𝐼 𝑛,𝑠𝑐𝑜 (𝜃), 𝜃 ∈ ℒ }.

Proof. Let us denote by 𝑆(𝑍 ) = (𝑆 1 (𝑍 1 ), … , 𝑆 𝑛 (𝑍 𝑛 )) the sufficient statistics of the model we consider in our approach. Let us also define 𝐻 (𝑍 , 𝑠) = 𝑆(𝑍 ) -𝑠 and ℎ(𝑠) = E 𝑍 |𝑌 ; θ (𝑠) (𝑆(𝑍 )) -𝑠. The proof is composed of two steps following for example the lines of (Allassonnière S., Kuhn E., and Trouvé A. 2010). First we establish the almost sure convergence of the sequence (𝑠 𝑘 ) generated by the algorithm toward the zero of the function ℎ. Second we deduce the almost sure convergence of the sequences (𝜃 𝑘 ) and (𝐼 𝑘 𝑛,𝑠𝑐𝑜 ) toward the set of critical points of the observed log-likelihood and the set ℐ respectively.

To prove the first step we apply Theorem 5.5 of (Andrieu C., Moulines E., and Priouret P. 2005). Therefore we have to verify that their four conditions denoted (𝐴1) to (𝐴4) are fulfilled. Our proof will follow the same global strategy as for example the one of Theorem 1 in [START_REF] Kuhn | Properties of the Stochastic Approximation EM Algorithm with Mini-Batch Sampling[END_REF]. We get first that condition (𝐴1) is satisfied by applying Lemma 2 of [START_REF] Delyon | Convergence of a Stochastic Approximation Version of the EM Algorithm[END_REF]). Indeed our assumptions (𝑀1 ′ ) and (𝑀2 ′ ) imply that assumptions (𝑀1) and (𝑀2) of [START_REF] Delyon | Convergence of a Stochastic Approximation Version of the EM Algorithm[END_REF]) are satisfied. These assumptions (𝑀1 ′ ) and (𝑀2 ′ ) focus on expressions and regularity properties of the individual likelihood functions and the corresponding sufficient statistics for each index 𝑖 ∈ {1, … , 𝑛}. The implication above follows by linearity of the log-likelihood function. Then we get that assumptions (𝐻 1), ( 𝑀1)-( 𝑀5), (𝑆𝐴𝐸𝑀2) ensured that condition (𝐴1) is satisfied. To prove assumptions (𝐴2) and (𝐴3), we will follow the strategy of (Allassonnière S., Kuhn E., and Trouvé A. 2010) to handle the difficulty of finding a common drift function 𝑉 for the family of posterior distributions indexed by 𝑠 ∈ 𝒮. Therefore we will construct first a family of drift functions (𝑉 𝑠 ) using Proposition 6.1 of (Andrieu C., Moulines E., and Priouret P. 2005), which stated drift conditions, called (𝐷𝑅𝐼 ), easy to check in practice. To prove condition (𝐷𝑅𝐼 1) for each kernels, we use Theorem 4.1 and Theorem 4.3 of [START_REF] Jarner | Geometric Ergodicity of Metropolis Algorithms[END_REF] which stated that assumptions (𝐻 2),(𝐻 3) and (𝐻 4) imply that Equations (6.1) and ( 6.3) of (𝐷𝑅𝐼 1) are satisfied with 𝑚 = 1 and 𝑉 𝑠 (𝑧) = 𝑝(𝑧; θ (𝑠)) -𝜌 with 𝜌 given by (𝐻 4). Then the common drift function is defined by 𝑉 (𝑧) = 𝑝(𝑧; 𝜃 0 ) -𝜌 using assumption (𝐻 4). Thus for any compact 𝒦 of Θ, there exist constants 𝑐 𝒦 > 0 and 𝐶 𝒦 > 0 such that for all 𝜃 ∈ 𝒦 and for all 𝑧, 𝑐 𝒦 𝑉 (𝑧) ≤ 𝑝(𝑧; θ (𝑠)) -𝜌 ≤ 𝐶 𝒦 𝑉 (𝑧). Therefore Equations (6.1) and (6.3) are satisfied for this common drift function 𝑉. Equation (6.2) follows also from Theorem 2.1 of [START_REF] Jarner | Geometric Ergodicity of Metropolis Algorithms[END_REF] which concludes the proof of (𝐷𝑅𝐼 1). The first part of (𝐷𝑅𝐼 2) is ensured by assumption (𝐻 4). The second part is satisfied directly with Lipschitz exponent 𝛽 equal to 1 in our case. Finally assumption (𝐷𝑅𝐼 3) is satisfied also with 𝛽 = 1 in our framework. This proof is obtained by using the usual strategy of splitting the whole space in four parts depending on the acceptance region and on the rejection region (see the proof of Lemma 4.7 in (Fort G. et al. 2015) for example) and the fact that the function θ is twice continuously differentiable. Finally assumption (𝑆𝐴𝐸𝑀1) allows to choose a sequence (𝜀 𝑘 ) such that (𝐴4) is satisfied (see constructive details in (Andrieu C., Moulines E., and Priouret P. 2005) after the statement of assumption (𝐴4)). This concludes the proof of the first step.

The function θ being continuous, we get that lim 𝑘 𝑑(𝜃 𝑘 , ℒ ) = 0 applying Lemma 2 of [START_REF] Delyon | Convergence of a Stochastic Approximation Version of the EM Algorithm[END_REF]). Moreover we get that for 1 ≤ 𝑖 ≤ 𝑛, each sequence (𝑠 𝑘 𝑖 ) converges almost surely toward E 𝑍 𝑖 |𝑌 𝑖 ;𝜃 (𝑆 𝑖 (𝑍 𝑖 )). Since assumption (𝑀2 ′ ) ensures that for all 1 ≤ 𝑖 ≤ 𝑛 the functions 𝜓 𝑖 and 𝜙 𝑖 are twice continuously differentiable and assumption (𝑀5) ensures that the function θ is continuously differentiable, the function

Φ 𝑛 defined by Φ 𝑛 (𝑠 𝑘 ) = 1 𝑛 ∑ 𝑛 𝑖=1 Δ𝑖 (𝑠 𝑘 ) Δ𝑖 (𝑠 𝑘
) is continuous. Therefore we get that lim 𝑘 𝑑(𝐼 𝑘 𝑛,𝑠𝑐𝑜 , ℐ ) = 0 which concludes the whole proof.

Description of the algorithm for general latent variables models

In general settings, the SAEM algorithm can yet be applied to approximate numerically the maximum likelihood estimate of the model parameter. Nevertheless there are no more theoretical guarantees of convergence for the algorithm. However we propose an extended version of our algorithm which allows to get an estimate of the Fisher information matrix as a by-product of the estimation algorithm.

Initialization step: Initialize arbitrarily Δ 0 𝑖 for all 1 ≤ 𝑖 ≤ 𝑛, 𝑄 0 and 𝜃 0 . Repeat until convergence the three steps defined at iteration 𝑘 by:

• Simulation step: for 1 ≤ 𝑖 ≤ 𝑛 simulate a realization 𝑍 𝑘 𝑖 directly from the conditional distribution given the observations 𝑌 𝑖 , denoted by 𝑝 𝑖 , or from a transition kernel of an ergodic Markov Chain having 𝑝 𝑖 as stationary distribution, using the current parameter 𝜃 𝑘-1 .

• Stochastic approximation step: compute the quantities for all 1 ≤ 𝑖 ≤ 𝑛

𝑄 𝑘 (𝜃) = (1 -𝛾 𝑘 )𝑄 𝑘-1 (𝜃) + 𝛾 𝑘 𝑛 ∑ 𝑖=1 log 𝑓 (𝑦 𝑖 , 𝑍 𝑘 𝑖 ; 𝜃) Δ 𝑘 𝑖 = (1 -𝛾 𝑘 )Δ 𝑘-1 𝑖 + 𝛾 𝑘 𝜕 𝜃 log 𝑓 (𝑦 𝑖 , 𝑍 𝑘 𝑖 ; 𝜃 𝑘-1 )
• Maximisation step: update of the parameter estimator according to:

𝜃 𝑘 = arg max 𝜃 𝑄 𝑘 (𝜃).
When convergence is reached, say at iteration 𝐾 of the algorithm, evaluate the FIM estimator according to:

𝐼 𝐾 𝑛,𝑠𝑐𝑜 = 1 𝑛 𝑛 ∑ 𝑖=1 Δ 𝐾 𝑖 (Δ 𝐾 𝑖 ) 𝑡 .
We illustrate through simulations in a nonlinear mixed effects model the performance of this algorithm in Section 4.2.

4 Simulation study 4.1 Asymptotic properties of the estimators 𝐼 𝑛,𝑠𝑐𝑜 (𝜃) and 𝐼 𝑛,𝑜𝑏𝑠 (𝜃)

In this section, we investigate the properties of the estimators 𝐼 𝑛,𝑠𝑐𝑜 (𝜃) and 𝐼 𝑛,𝑜𝑏𝑠 (𝜃) when the sample size 𝑛 grows.

Simulation settings

First we consider the following linear mixed effects model 𝑦 𝑖𝑗 = 𝛽 + 𝑧 𝑖 + 𝜀 𝑖𝑗 , where 𝑦 𝑖𝑗 ∈ ℝ denotes the 𝑗 𝑡ℎ observation of individual 𝑖, 1 ≤ 𝑖 ≤ 𝑛, 1 ≤ 𝑗 ≤ 𝐽, 𝑧 𝑖 ∈ ℝ the unobserved random effect of individual 𝑖 and 𝜀 𝑖𝑗 ∈ ℝ the residual term. The random effects (𝑧 𝑖 ) are assumed independent and identically distributed such that 𝑧 𝑖 ∼ 𝑖.𝑖.𝑑.

𝒩 (0, 𝜂 2 ), the residuals (𝜀 𝑖𝑗 ) are assumed independent and identically distributed such that 𝜀 𝑖𝑗 ∼ 𝑖.𝑖.𝑑.

𝒩 (0, 𝜎 2 ) and the sequences (𝑧 𝑖 ) and (𝜀 𝑖𝑗 ) are assumed mutually independent. Here, the model parameters are 𝜃 = (𝛽, 𝜂 2 , 𝜎 2 ). We set 𝛽 = 3, 𝜂 2 = 2, 𝜎 2 = 5 and 𝐽 = 12.

Second we consider the following Poisson mixture model where the distribution of each observation 𝑦 𝑖 , 1 ≤ 𝑖 ≤ 𝑛, depends on a state variable 𝑧 𝑖 which is latent leading to

𝑦 𝑖 |𝑧 𝑖 = 𝑘 ∼ 𝒫 (𝜆 𝑘 ) with 𝑃(𝑧 𝑖 = 𝑘) = 𝛼 𝑘 and ∑ 𝐾 𝑘=1 𝛼 𝑘 = 1. The model parameters are 𝜃 = (𝜆 1 , … , 𝜆 𝐾 , 𝛼 1 , … , 𝛼 𝐾 -1 ).
For the simulation study, we consider a mixture of 𝐾 = 3 components, and the following values for the parameters 𝜆 1 = 2, 𝜆 2 = 5, 𝜆 3 = 9, 𝛼 1 = 0.3 and 𝛼 2 = 0.5.

For each model, we generate 𝑀 = 500 datasets for different sample sizes 𝑛 ∈ {20, 100, 500}. As a first step, we assume that the true parameter values, denoted by 𝜃 ⋆ , are known in order to investigate the asymptotic properties of both 𝐼 𝑛,𝑠𝑐𝑜 and 𝐼 𝑛,𝑜𝑏𝑠 without adding extra noise induced by the estimation of the parameters. Hence, for each value of 𝑛 and for each 1 ≤ 𝑚 ≤ 𝑀, we derive 𝐼 (𝑚) 𝑛,𝑠𝑐𝑜 (𝜃 ⋆ ) and 𝐼 (𝑚) 𝑛,𝑜𝑏𝑠 (𝜃 ⋆ ). The estimators 𝐼 𝑛,𝑠𝑐𝑜 (𝜃 ⋆ ) and 𝐼 𝑛,𝑜𝑏𝑠 (𝜃 ⋆ ) can be computed explicitly in both models by applying Equation 5and Louis' formula [START_REF] Louis | Finding the Observed Information Matrix When Using the EM Algorithm[END_REF]) (see R functions provided in the Appendix section). We then compute the empirical bias and the root mean squared deviation of each component (ℓ, ℓ ′ ) of the estimated matrix as:

1 𝑀 𝑀 ∑ 𝑚=1 𝐼 (𝑚) 𝑛,𝑠𝑐𝑜,ℓ,ℓ ′ (𝜃 ⋆ ) -𝐼 ℓ,ℓ ′ (𝜃 ⋆ ) and √ 1 𝑀 𝑀 ∑ 𝑚=1 (𝐼 (𝑚) 𝑛,𝑠𝑐𝑜,ℓ,ℓ ′ (𝜃 ⋆ ) -𝐼 ℓ,ℓ ′ (𝜃 ⋆ )) 2 .
In the previous quantities, 𝐼 (𝜃 ⋆ ) is explicit in the linear mixed effects model and approximated by a Monte-Carlo estimation based on a sample of size 10 8 in the Poisson mixture model. The results are presented in Table 1 and Table 2 for the linear mixed effects model and in Table 3 and Table 4 for the mixture model. In a second step, we use the linear mixed effects model to look at what happens when the parameter 𝜃 is unknown and the estimation of the Fisher information matrix requires to compute an estimate θ 𝑛 . We use the datasets simulated with 𝑛 = 500 and we compute the 𝑀 = 500 asymptotic confidence intervals of the three model parameters. We then deduce empirical coverage rates for the following nominal rates 1 -𝛼 ∈ {0.90, 0.95, 0.99} by using the diagonal terms of either the inversed 𝐼

(𝑚) 𝑛,𝑠𝑐𝑜 (𝜃 ⋆ ) (resp. 𝐼 (𝑚) 𝑛,𝑜𝑏𝑠 (𝜃 ⋆ )) or the inversed 𝐼 (𝑚) 𝑛,𝑠𝑐𝑜 ( θ 𝑛 ) (resp. 𝐼 (𝑚) 𝑛,𝑜𝑏𝑠 ( θ 𝑛 )).
The results are depicted in Table 5 andTable 6. ## 1-R script for studying the asymptotic properties of Iobs and Isco in the ## linear mixed effects model

## ----------------------------------------------------------------------------- library(lme4) nsim <-500
# number of replicates seq.n <-c(20,100,500) # number of individuals j <-12 # number of observations per individual ## parameter values beta <-3 sigma2 <-5 eta2 <-2 theta.true <-matrix(c(beta,eta2,sigma2),ncol=1) ## R objects to store estimations resIobs.theta.true <-array(NA,dim=c(3,3,nsim,length(seq.n))) resIsco.theta.true <-array(NA,dim=c(3,3,nsim,length(seq.n))) resIobs.theta.est <-array(NA,dim=c(3,3,nsim,length(seq.n))) resIsco.theta.est <-array(NA,dim=c(3,3,nsim,length(seq.n))) --------------------------------------------------------------------------- ---------------------------------------------------------------------------

EstF11.true <-c() EstF22.true <-c() EstF33.true <-c() EstF12.true <-c() EstF13.true <-c() EstF23.true <-c() EstF11.est <-c() EstF22.est <-c() EstF33.est <-c() EstF12.est <-c() EstF13.est <-c() EstF23.est <-c() beta.est <-c() eta2.est <-c() sigma2.est <-c() ## loop executing the nsim replicates of the experiment for (l in 1:length(seq.n)){ n <-seq.n[l] beta.est.n <-c() eta2.est.n <-c() sigma2.est.n <-c() for (k in 1:nsim){ ## data simulation random <-rnorm(n,0,sqrt(eta2)) residual <-rnorm(n*j,0,sqrt(sigma2)) randompop <-rep(random,j) id <-rep(seq(1,n),

Results

From Table 1, Table 2, Table 3 andTable 4, we observe that whatever the model and whatever the components of 𝐼 𝑛,𝑠𝑐𝑜 (𝜃 ⋆ ) and 𝐼 𝑛,𝑜𝑏𝑠 (𝜃 ⋆ ), the bias is very small even for small values of 𝑛. Note that in the linear mixed effects model the second derivative with respect to parameter 𝛽 is deterministic, which explains why the bias and the dispersion of the estimations 𝐼 𝑛,𝑜𝑏𝑠 (𝜃 ⋆ ) are zero for every value of 𝑛. The bias and the standard deviation decrease as 𝑛 increases overall, which illustrates the consistency of both M-estimators. The distributions of the normalized estimations √ 𝑛 (𝐼

(𝑚) 𝑛,𝑠𝑐𝑜 (𝜃 ⋆ ) -𝐼 (𝜃 ⋆ )) and √ 𝑛 (𝐼 (𝑚)
𝑛,𝑜𝑏𝑠 (𝜃 ⋆ ) -𝐼 (𝜃 ⋆ )) are also represented when 𝑛 = 500 for some components of the matrices in Figure 1 (linear mixed effects model) and Figure 2 (Poisson mixture model). The empirical distributions have the shape of Gaussian distributions and illustrate the asymptotic normality of the two estimators. The numerical results highlight that neither 𝐼 𝑛,𝑠𝑐𝑜 (𝜃 ⋆ ) nor 𝐼 𝑛,𝑜𝑏𝑠 (𝜃 ⋆ ) is systematically better than the other one in terms of bias and asymptotic covariance matrix. In the same model, different behaviors can be observed depending on the components of the parameter vector. Table 5 and Table 6 show that the empirical coverage rates computed from 𝐼 𝑛,𝑠𝑐𝑜 and 𝐼 𝑛,𝑜𝑏𝑠 in the linear mixed effects model are close to the nominal values, which corroborates the relevance of both estimators. Moreover there is little difference between the results obtained when using 𝐼 𝑛,𝑠𝑐𝑜 or 𝐼 𝑛,𝑜𝑏𝑠 to estimate the Fisher information matrix. When the parameter value is unknown, the uncertainty related to the parameter estimation leads to a deterioration of the coverage rates. Still, this deterioration diminishes when 𝑛 increases.

Asymptotic properties of the stochastic approximation algorithm

We now investigate the properties of our algorithm with truncation on random boundaries in the curved exponential family when the number of iterations grows (Section 4.2.1) and the good performance of its extended version in more general latent variable models (Section 4.2.2). We also present a short comparison with existing methods (Section 4.2.3).

In curved exponential family models

We consider the following nonlinear mixed effects model which is widely used in pharmacokinetics for describing the evolution of drug concentration over time:

𝑦 𝑖𝑗 = 𝑔 𝑖 (𝑡 𝑖𝑗 , 𝑧 𝑖 ) + 𝜀 𝑖𝑗 , (7)
where ). In this model, as in a large majority of nonlinear mixed effects models, the likelihood does not have any analytical expression. As a consequence, neither the Fisher Information Matrix, nor the estimators 𝐼 𝑛,𝑠𝑐𝑜 (𝜃), 𝐼 𝑛,𝑜𝑏𝑠 (𝜃) have explicit expressions. However, as the complete data log-likelihood is explicit, stochastic approximations of 𝐼 𝑛,𝑠𝑐𝑜 (𝜃), 𝐼 𝑛,𝑜𝑏𝑠 (𝜃) can be implemented. Note moreover that this model belongs to the curved exponential family as defined in Equation 6with

𝑆 𝑖 (𝑧 𝑖 ) = ( 𝐽 ∑ 𝑗=1 (𝑦 𝑖𝑗 𝑔 𝑖 (𝑡 𝑖𝑗 , 𝑧 𝑖 )), (log 𝑘𝑎 𝑖 ), (log 𝐶𝑙 𝑖 ), (log 𝑉 𝑖 ), (log 𝑘𝑎 𝑖 ) 2 , (log 𝐶𝑙 𝑖 ) 2 , (log 𝑉 𝑖 ) 2 ) ′ 𝜙 𝑖 (𝜃) = ( 1 2𝜎 2 , log 𝑘𝑎 𝜔 2 𝑘𝑎 , log 𝐶𝑙 𝜔 2 𝐶𝑙 , log 𝑉 𝜔 2 𝑉 , - 1 2𝜔 2 𝑘𝑎 , - 1 2𝜔 2 𝐶𝑙 , - 1 2𝜔 2 𝑉 ) , 𝜓 𝑖 (𝜃) = 1 2 ( (log 𝑘𝑎) 2 𝜔 2 𝑘𝑎 + (log 𝐶𝑙) 2 𝜔 2 𝐶𝑙 + (log 𝑉 ) 2 𝜔 2 𝑉
) .

The algorithm described in Section 3.2.1 can therefore be easily implemented to estimate 𝜃 and the Fisher information matrix simultaneously (see R function provided in the Appendix section).

We take the following values for the parameters 𝑉 = 31, 𝑘𝑎 = 1.6, 𝐶𝑙 = 2.8, 𝜔 2 𝑉 = 0.40, 𝜔 2 𝑘𝑎 = 0.40, 𝜔 2 𝐶𝑙 = 0.40 and 𝜎 2 = 0.75. We consider the same dose 𝑑 𝑖 = 320 and the same observation times (in hours): 0.25,0.5, 1, 2, 3.5, 5, 7, 9, 12, 24 for all the individuals. We simulate one dataset with 𝑛 = 100 individuals under model specified by Equation 7. On this simulated dataset, we run 𝑀 = 500 times the stochastic approximation algorithm described in Section 3.2.1 for computing 𝐼 𝑛,𝑠𝑐𝑜 ( θ ) together with θ and the algorithm of (Delyon B., Lavielle M., and Moulines E. 1999) for computing 𝐼 𝑛,𝑜𝑏𝑠 ( θ ). We perform 𝐾 = 3000 iterations in total for each algorithm by setting 𝛾 𝑘 = 0.95 for 1 ≤ 𝑘 ≤ 1000 (burn in iterations) and 𝛾 𝑘 = (𝑘 -1000) -3/5 otherwise, 𝜀 𝑘 = 5. At any iteration, we compute the empirical relative bias and the empirical relative standard deviation of each component (ℓ, ℓ ′ ) of 𝐼 𝑛,𝑠𝑐𝑜 defined respectively as:

1 𝑀 𝑀 ∑ 𝑚=1 Î (𝑘,𝑚) 𝑛,𝑠𝑐𝑜,ℓ,ℓ ′ -𝐼 ⋆ 𝑛,𝑠𝑐𝑜,ℓ,ℓ ′ 𝐼 ⋆ 𝑛,𝑠𝑐𝑜,ℓ,ℓ ′ and √ √ √ √ √ √ √ √ √ 1 𝑀 𝑀 ∑ 𝑚=1 ( Î (𝑘,𝑚) 𝑛,𝑠𝑐𝑜,ℓ,ℓ ′ -𝐼 ⋆ 𝑛,𝑠𝑐𝑜,ℓ,ℓ ′ 𝐼 ⋆ 𝑛,𝑠𝑐𝑜,ℓ,ℓ ′ ) 2
where Î (𝑘,𝑚) 𝑛,𝑠𝑐𝑜 denotes the estimated value of 𝐼 𝑛,𝑠𝑐𝑜 ( θ ) at iteration 𝑘 of the 𝑚 𝑡ℎ algorithm. We compute the same quantities for 𝐼 𝑛,𝑜𝑏𝑠 . As the true values of 𝐼 ⋆ 𝑛,𝑠𝑐𝑜 = 𝐼 𝑛,𝑠𝑐𝑜 (𝜃 ⋆ ) and 𝐼 ⋆ 𝑛,𝑜𝑏𝑠 = 𝐼 𝑛,𝑜𝑏𝑠 (𝜃 ⋆ ) are not known, they are estimated by Monte-Carlo integration based on 10 5 iterations, including 5000 burnin, of a Metropolis-Hastings algorithm.

## R script for studying the properties of the SAEM algorithm in the curved ## exponential family when the number of iterations grow. ## - ----------------------------------------------------------------------------## 1-Data simulation # Sample characteristics set.seed(3005) n <-100 # number of subjects times <-c(0.25,0.5,1,2,3. 5,5,7,9,12,24) [4], omega2v=tm [5],omega2cl=tm [6],sigma2=tm [7])

FisherMC <-FIM_mc(datasim, nbMC, nbMCburnin, thetaMean) iscoMC <-FisherMC$iscoMC iobsMC <-FisherMC$iobsMC # Evaluation of the mean relative bias and of the mean relative standard errors # per iteration.

biasIsco <-array(0,dim=c (7,7,nbsim,nbiterem))

for (j in 1:nbsim){ for (k in 1:nbiterem){ biasIsco[,,j,k] <-(iscoarray[j,,,k] -iscoMC)/iscoMC } } biasIobs <-array(0,dim=c (7,7,nbsim,nbiterem)) for (j in 1:nbsim){ for (k in 1:nbiterem){ biasIobs[,,j,k] <-(iobsarray[j,,,k] -iobsMC)/iobsMC } } rsdIsco <-array(0,dim=c (7,7,nbsim,nbiterem)) for (j in 1:nbsim){ for (k in 1:nbiterem){ rsdIsco[,,j,k] <-(iscoarray[j,,,k] -iscoMC)^2/iscoMC^2 } } rsdIobs <-array(0,dim=c (7,7,nbsim,nbiterem)) The results are displayed in Figure 3 and Figure 4.

for (j in 1:nbsim){ for (k in 1:nbiterem){ rsdIobs[,,j,k] <-(iobsarray[j,,,k] -iobsMC)^2/iobsMC^2 } } MbiasIsco <-apply(biasIsco[,,, (nbiterburnin+1)
We observe that the bias and the standard deviations of the estimates of the components of both ( ω Cl 2 , ω Cl 2 )

Figure 4: Non linear mixed effects model. Representation over iterations of the mean relative standard error of the diagonal components of the estimated Fisher information matrix computed from the 𝑀 = 500 runs of the stochastic algorithm. Red line corresponds to 𝐼 𝑛,𝑠𝑐𝑜 (𝜃) and blue line corresponds to 𝐼 𝑛,𝑜𝑏𝑠 (𝜃). The burn-in iterations of the algorithme are not depicted matrices decrease over iterations, and that for both estimates the bias is nearly zero when the convergence of the algorithm is reached. According to these simulation results, there is no evidence that one method is better than the other in terms of bias or standard deviation.

In general latent variable models

We use model specified by Equation 7again, but we now consider that individual parameter 𝑉 𝑖 is fixed, i.e. 𝑉 𝑖 ≡ 𝑉 ∀𝑖 = 1, … , 𝑛. The model is no longer exponential in the sense of Equation 6. We must therefore use the general version of the stochastic approximation algorithm from Section 3.2.3 to compute 𝐼 𝑛,𝑠𝑐𝑜 ( θ ) (see R function provided in the Appendix section). We simulate 500 datasets according to this model and we estimate 𝐼 𝑛,𝑠𝑐𝑜 ( θ ) and θ for each one. We perform 𝐾 = 3000 iterations of the algorithm by setting 𝛾 𝑘 = 𝑘 -0.501 . We compute the 500 asymptotic confidence intervals of the model parameters, by using either the inversed 𝐼 𝑛,𝑠𝑐𝑜 ( θ 𝑘 )'s or the inversed 𝐼 𝑛,𝑜𝑏𝑠 ( θ 𝑘 )'s and then deduce from them empirical coverage rates.

## R script for studying the relevance of the SAEM algorithm out of the curved ## exponential family ## - ----------------------------------------------------------------------------## 1-Data simulation ## Sample characteristics n <-100 # number of subjects times <-c(0.25,0.5,1,2,3.5,5,7,9,12,24) (6,6,nbsim)) iobs.est <-array(0,dim=c (6,6,nbsim)) for (kk in 1:nbsim){ set.seed(kk*2500+10) ## Simulation of individual parameters vi <-rep(vpop,n) kai <-exp(rnorm(n,log(kapop),sd=sqrt(omega2ka))) cli <-exp(rnorm(n,log(clpop),sd=sqrt(omega2cl))) ## Simulation of the observations ypred <-c() for (k in 1:n){ ypred <-c(ypred,model1cptsim(cbind(kai,vi,cli),k,times,dose)) } y <-ypred + rnorm(n*j,0,sd=sqrt(sigma2)) We obtain for the six parameters (𝑘𝑎, 𝑉 , 𝐶𝑙, 𝜔 2 𝑘𝑎 , 𝜔 2 𝐶𝑙 , 𝜎 2 ) empirical covering rates of 0.96, 0.948, 0.948, 0.932, 0.948, 0.948 respectively for a nominal covering rate of 0.95. This highlights that our estimate accurately quantifies the precisions of parameter estimates. Note that empirical coverage rates computed from 𝐼 𝑛,𝑜𝑏𝑠 are similar (here 0.952, 0.93, 0.942, 0.924, 0.952, 0.946) but that the real advantage of our method is that it requires stochastic approximation only on the first-order derivatives of the complete log-likelihood, contrary to 𝐼 𝑛,𝑜𝑏𝑠 which requires deriving the complete log-likelihood at the second order and thus implies more complicated formulas since the model does not belong to the exponential family.

Convergence graphs obtained from a simulated data set are shown in Figure 5. Although theoretical guarantee is missing in non exponential models, the stochastic approximation algorithm proposed in Section 3.2.3 converges in practice on this example for both the estimation of the model parameters and the estimation of the Fisher information matrix.

Comparison with other methods

To the best of our knowledge, although there exists contributions focusing on the estimation of the Fisher information matrix in latent variable models, there is currently no method based on the first iscoest <-fisher_estimation_gaussian_mixture(y,est$m1,est$m2,est$prob) -list(isco=round(apply(isco,c(2,3),mean),3),recouvprob=recouvprob/nrep, recouvm1=recouvm1/nrep,recouvm2=recouvm2/nrep) save(res,file="Rfiles/ResGaussianMixture.Rdata")

ICinf <- c(est$prob,est$m1,est$m2) -qnorm(1-(1-rate)/2)*sqrt(diag(solve(iscoest))) ICsup <- c(est$prob,est$m1,est$m2) + qnorm(1-(1-rate)/2)*sqrt(diag(solve(iscoest))) if ((probtrue>=ICinf[1])&(probtrue<=ICsup[1])){recouvprob <-recouvprob + 1} if ((m1true>=ICinf[2])&(m1true<=ICsup[2])){recouvm1 <-recouvm1 + 1} if ((m2true>=ICinf[3])&(m2true<=ICsup[3])){recouvm2 <-recouvm2 + 1} isco[j,,] <-iscoest } res <
We simulate 10000 datasets according to this Gaussian mixture model, using the same setting as [START_REF] Meng | Efficient Computation of the Fisher Information Matrix in the EM Algorithm[END_REF], i.e. 𝑛 = 750, 𝜋 = 2/3, 𝜇 1 = 3 and 𝜇 2 = 0. For each dataset 𝑘 = 1, … , 10000, we compute the parameter maximum likelihood estimate θ 𝑘 = ( π 𝑘 , μ1𝑘 , μ2𝑘 ) with an EM algorithm and then we derive 𝐼 𝑛,𝑠𝑐𝑜 ( θ 𝑘 ) directly according to Equation 5(see R function provided in the Appendix section) contrary to [START_REF] Meng | Efficient Computation of the Fisher Information Matrix in the EM Algorithm[END_REF] who used an iterative method. We compute the empirical mean of the 10000 estimated matrices leading to: 2687.873 -210.795 -251.634 -210.795 170.9 -61.546 -251.634 -61.546 393.115 ) .

1 10000 ∑ 𝑘 𝐼 𝑛,𝑠𝑐𝑜 ( θ 𝑘 ) = (
Comparison with the results of (Meng L. and Spall J. C. 2017) is delicate since their numerical illustration of their method is based on a single simulated dataset thus potentially sensitive to sampling variations. However, they provide an estimation of the Fisher information matrix from this unique dataset 2591.3 -237.9 -231.8 -237.9 155.8 -86.7 -231.8 -86.7 394.5 ) .

𝐼 𝑀𝑒𝑛𝑔 = (
Our results are coherent with their ones. To check the reliability of our results, we then compute as above the 10000 asymptotic confidence intervals of the three model parameters. We obtain for the three parameters (𝜋, 𝜇 1 , 𝜇 2 ) empirical covering rates of 0.9477, 0.9499, 0.9523 respectively for a nominal covering rate of 0.95. Thus 𝐼 𝑛,𝑠𝑐𝑜 accurately quantifies the precisions of parameter estimates.

Conclusion and discussion

In this work, we address the estimation of the Fisher information matrix in general latent variable models. We focus on the empirical Fisher information matrix which is a moment estimate of the covariance matrix of the score. We propose stochastic approximation algorithms to compute this estimate when it can not be calculated analytically and establish theoretical convergence properties in the curved exponential family setting. We carry out a simulation study in mixed effects model and in a Poisson mixture model to compare the performances of several estimates, namely the considered empirical Fisher information matrix and the observed Fisher information matrix. We emphasize that the empirical FIM requires less regularity assumptions than the observed FIM. From a computational point of view, the implementation of the algorithm for evaluating the empirical FIM only involves the first derivatives of the log-likelihood, in contrary to the one for evaluating the observed FIM which involves the second derivatives of the log-likelihood.

The main perspective of this work is to adapt the procedure for statistical models whose derivatives of the log-likelihood have no tractable expressions, coupling the algorithm with numerical derivative procedures.

Appendix

Description of the algorithm without truncation on random boundaries in curved exponential family model

We provide here a simpler algorithm based on an extension of the stochastic approximation Expectation Maximization algorithm proposed by [START_REF] Delyon | Convergence of a Stochastic Approximation Version of the EM Algorithm[END_REF]) using a simulation step performed from the conditional distribution and without truncation on random boundaries. Theoretical results are established assuming a stability condition which is usually quite difficult to check. However, this algorithm can be easily applied in practice.

Initialization step:

Initialize arbitrarily for all 1 ≤ 𝑖 ≤ 𝑛 𝑠 0 𝑖 and 𝜃 0 . Repeat until convergence the three steps defined at iteration 𝑘 by:

• Simulation step: for 1 ≤ 𝑖 ≤ 𝑛 simulate a realization 𝑍 𝑘 𝑖 from the conditional distribution given the observations 𝑌 𝑖 denoted by 𝑝 𝑖 using the current parameter value 𝜃 𝑘-1 .

• Stochastic approximation step: compute the quantities for all 1 ≤ 𝑖 ≤ 𝑛

𝑠 𝑘 𝑖 = (1 -𝛾 𝑘 )𝑠 𝑘-1 𝑖 + 𝛾 𝑘 𝑆 𝑖 (𝑍 𝑘 𝑖 )
where (𝛾 𝑘 ) is a sequence of positive step sizes satisfying ∑ 𝛾 𝑘 = ∞ and ∑ 𝛾 2 𝑘 < ∞. • Maximisation step: update of the parameter estimator according to:

𝜃 𝑘 = arg max 𝜃 𝑛 ∑ 𝑖=1 (-𝜓 𝑖 (𝜃) + ⟨𝑠 𝑘 𝑖 , 𝜙 𝑖 (𝜃)⟩) = θ (𝑠 𝑘 )
When convergence is reached, say at iteration 𝐾 of the algorithm, evaluate the FIM estimator according to:

𝐼 𝐾 𝑛,𝑠𝑐𝑜 = 1 𝑛 𝑛 ∑ 𝑖=1 Δ𝑖 (𝑠 𝐾 ) Δ𝑖 (𝑠 𝐾 ) 𝑡
where Δ𝑖 (𝑠) = -𝜕𝜓 𝑖 ( θ (𝑠)) + ⟨𝑠 𝑖 , 𝜕𝜙 𝑖 ( θ (𝑠))⟩ for all 𝑠.

Remark. In the cases where the latent variables can not be simulated from the conditional distribution, one can apply the extension coupling the stochastic algorithm with a Monte Carlo Markov Chain procedure as presented in [START_REF] Kuhn | Coupling a Stochastic Approximation Version of EM with an MCMC Procedure[END_REF]. All the following results can be extended to this case.

Theoretical convergence properties

The theoretical following results provide convergence guarantees for the FIM estimate obtained as a by-product of the MLE. Therefore they extend those of [START_REF] Delyon | Convergence of a Stochastic Approximation Version of the EM Algorithm[END_REF].

To that purpose, in addition to the exponential family assumption for each individual likelihood, we also make the same type of regularity assumptions as those presented in [START_REF] Delyon | Convergence of a Stochastic Approximation Version of the EM Algorithm[END_REF]) at each individual level. These regularity assumptions on the model are detailed at the end of the appendix section. 𝑠 𝑙 /𝑘 using the results stated by [START_REF] Delyon | Convergence of a Stochastic Approximation Version of the EM Algorithm[END_REF]). Let us denote by 𝑉 𝑒𝑐𝑡(𝐴) the vector composed of the elements of the triangular superior part of matrix 𝐴 ordered by columns. Theorem 6.2. Assume that (𝑀1 ′ ) and (𝑀2 ′ ), (𝑀3) to (𝑀5), (𝑆𝐴𝐸𝑀1), (𝑆𝐴𝐸𝑀2), (𝑆𝐴𝐸𝑀3), (𝑆𝐴𝐸𝑀4), (𝑆𝐴𝐸𝑀4 ′ ) and (𝐿𝑂𝐶1) to (𝐿𝑂𝐶3) are fulfilled. Then, there exists a regular stable stationary point 𝜃 * ∈ Θ such that lim 𝑘 𝜃 𝑘 = 𝜃 * a.s. Moreover the sequence (√𝑘(𝑉 𝑒𝑐𝑡( Ī 𝑘 𝑛,𝑠𝑐𝑜 ) -𝑉 𝑒𝑐𝑡( Ī𝑛,𝑠𝑐𝑜 (𝜃 * ))))𝟙 lim ‖𝜃 𝑘 -𝜃 * ‖=0 converges in distribution toward a centered Gaussian random vector when 𝑘 goes to infinity. The asymptotic covariance matrix is characterized.

Proof. The proof follows the lines of this of Theorem 7 of [START_REF] Delyon | Convergence of a Stochastic Approximation Version of the EM Algorithm[END_REF]. Assumptions (𝐿𝑂𝐶1) to (𝐿𝑂𝐶3) are those of [START_REF] Delyon | Convergence of a Stochastic Approximation Version of the EM Algorithm[END_REF] and ensure the existence of a regular stable stationary point 𝑠 * for ℎ and therefore of 𝜃 * = θ (𝑠 * ) for the observed log-likelihood 𝑙. Then applying Theorem 4 of [START_REF] Delyon | Convergence of a Stochastic Approximation Version of the EM Algorithm[END_REF], we get that:

√𝑘( s 𝑘 -𝑠 * )𝟙 lim ‖𝑠 𝑘 -𝑠 * ‖=0 ℒ → 𝒩 (0, 𝐽 (𝑠 * ) -1 Γ(𝑠 * )𝐽 (𝑠 * ) -1 )𝟙 lim ‖𝑠 𝑘 -𝑠 * ‖=0
where the function Γ defined in assumption (𝑆𝐴𝐸𝑀4 ′ ) and 𝐽 is the Jacobian matrix of the function ℎ.

Applying the Delta method, we get that:

√𝑘(𝑉 𝑒𝑐𝑡(Φ 𝑛 ( s 𝑘 )) -𝑉 𝑒𝑐𝑡(Φ 𝑛 (𝑠 * )))𝟙 lim ‖𝑠 𝑘 -𝑠 * ‖=0 ℒ → 𝑊 𝟙 lim ‖𝑠 𝑘 -𝑠 * ‖=0
where 𝑊 ∼ 𝒩 (0, 𝜕𝑉 𝑒𝑐𝑡(Φ 𝑛 (𝑠 * ))𝐽 (𝑠 * ) -1 Γ(𝑠 * )𝐽 (𝑠 * ) -1 𝜕𝑉 𝑒𝑐𝑡(Φ 𝑛 (𝑠 * )) 𝑡 ) which leads to the result.

Note that as usually in stochastic approximation results, the rate √𝑘 is achieved when considering an average estimator (see Theorem 7 of (Delyon B., Lavielle M., and Moulines E. 1999) e.g).

It is assumed that the random variables 𝑠 0 , 𝑧 1 , 𝑧 2 , ⋯ are defined on the same probability space (Ω, 𝒜 , 𝑃).

We denote ℱ = {ℱ 𝑘 } 𝑘≥0 the increasing family of 𝜎-algebras generated by the random variables 𝑠 0 , 𝑧 1 , 𝑧 2 , ⋯ , 𝑧 𝑘 . We assume the following conditions:

• (M1') • (M5) There exists a continuously differentiable function θ ∶ 𝒮 → Θ, such that:

∀𝑠 ∈ 𝒮 , ∀𝜃 ∈ Θ, 𝐿(𝑠; θ (𝑠)) ≥ 𝐿(𝑠; 𝜃).

In addition, we define: We also define assumptions required for the normality result:

• (SAEM1) For all 𝑘 in ℕ, 𝛾 𝑘 ∈ [0, 1], ∑ ∞ 𝑘=1 𝛾 𝑘 = ∞ and ∑ ∞ 𝑘=1 𝛾 2 𝑘 < ∞. • (
• (SAEM1') There exist 𝛾 * > 0 and 1/2 < 𝛼 < 1 such that lim 𝑘 𝛼 /𝛾 𝑘 = 𝛾 * , and 𝛾 𝑘 /𝛾 𝑘+1 = 1 + 𝑂(𝑘 -1 ).

• (SAEM4') For some 𝜀 > 0, sup 𝜃 E 𝜃 (‖𝑆(𝑍 )‖ 2+𝜀 ) < ∞ and 𝜃 → Γ(𝜃) is continuous w.r.t. 𝜃.

• (LOC1) The stationary points of 𝑙 are isolated: any compact subset of Θ contains only a finite number of such points.

• : value of the random effects variance # j : number of observations per individual crochet <-2*j*eta2/(sigma2+eta2*j)^2/sigma2^3 + 4/(sigma2+eta2*j)^2/sigma2^2 + 2/(sigma2+eta2*j)^3/sigma2 alpha <-j*(eta2+sigma2)/sigma2^3 -eta2/2*j*(j*eta2+sigma2)*crochet -(j-1)/2/sigma2^2 -1/2/(sigma2+eta2*j)^2 fisher <-cbind(c(j/(sigma2+eta2*j),0,0), c(0,j^2/2/(sigma2+eta2*j)^2,j/2/(sigma2+eta2*j)^2), c(0,j/2/(sigma2+eta2*j)^2,alpha)) return(fisher) }

(K-1)){ deriv1ind[k,] <- exp(-lambda[k])*lambda[k]^y*alpha[k]/denom*(y/lambda[k]-1) deriv1ind[K+k,] <-exp(-lambda[k])*lambda[k]^y/denom - exp(-lambda[K])*lambda[K]^y/denom } deriv1ind[K,] <- exp(-lambda[K])*lambda[K]^y*(1-sum(alpha))/denom*(y/lambda[K]-1)
## computation of conditional expectation of the second derivatives of the ## complete data log-likelihood for (k in 1: (kapop,vpop,clpop,omega2ka,omega2v,omega2cl,sigma2) # variances of the proposal kernels of the MCMC procedure eta2v <-kRW*omega2v eta2ka <-kRW*omega2ka eta2cl <-kRW*omega2cl

(K-1)){ deriv2[k,k] <- sum(exp(-lambda[k])*lambda[k]^y*alpha[k]/denom*(-y/lambda[k]^2)) deriv2[K+k,K+k] <- sum(-exp(-lambda[k])*
<-theta0$vpop kapop <-theta0$kapop clpop <-theta0$clpop omega2v <-theta0$omega2v omega2ka <-theta0$omega2ka omega2cl <-theta0$omega2cl sigma2 <-theta0$sigma2 p <-length(theta0) thetaest <-matrix(0,p,nbiterem) thetaest[,1] <-c
# sequence of step sizes gamma <-c(rep(0.95,nbiterburnin), 1/(2:nbiterem)^0. <-array(0,c(p,p,nbiterem)) tempderiveeas <-matrix(0,p,n) tempderiveeas2 <-matrix(0,p,p) dimstatexh <-7 statexh <-array(0,dim=c(dimstatexh,n,nbiterem)) psi <-array(0,dim=c(n,nb.psi,nbiterem))

# initial values for the individual parameters currentv <-log(vpop) + rnorm(n,0,sqrt(eta2v)) currentka <-log(kapop) + rnorm(n,0,sqrt(eta2ka)) currentcl <-log(clpop) + rnorm(n,0,sqrt(eta2cl)) currentpsi <-cbind(exp(currentka),exp(currentv),exp(currentcl)) psi [,,1] <-psiinit <-currentpsi while(compact(statexh[,,1],kappa)==0){ currentv <-log(vpop) + rnorm(n,0,sqrt(eta2v)) currentka <-log(kapop) + rnorm(n,0,sqrt(eta2ka)) currentcl <-log(clpop) + rnorm(n,0,sqrt(eta2cl)) currentpsi <-cbind(exp(currentka),exp(currentv),exp(currentcl)) psi [,,1] <-psiinit <-currentpsi mco <-matrix((y-model1cpt(psi,id,xidep))^2,n,j,byrow=TRUE) statexh[1:3,,1] <-t(log(psi)) statexh [4:6, <-statexh [4,l+1]-statexh[2,l+1]^2 sigma2 <-statexh [5,l+1]/n/j thetaest[,l+1] <-c (kapop, vpop, clpop, omega2ka, omega2cl, sigma2) eta2cl <-kRW*omega2cl eta2ka <-kRW*omega2ka ## Stochastic approximation of the derivatives of the complete log-likelihood

  -nMC <-100000000 alpha <-c(0.3,0.5) # mixture weights of the first K-1 components lambda <-c(2,5,9) # parameter values of the K Poisson distributions y <-sim_poisson_mixture(nMC,lambda,alpha) trueFIM <-fisher_estimation_poisson_mixture(y, nMC, lambda, alpha) trueFIM <-(trueFIM$Isco+trueFIM$Iobs)/2 save(trueFIM,file='Rfiles/PoissonMixtureTrueFIM.Rdata') ## 2-R script for studying the asymptotic properties of Iobs and Isco in the ## Poisson mixture model ## -

Figure 1 :

 1 Figure 1: Linear mixed effects model. Kernel density estimates of the normalized values of some components of the estimated Fisher information matrix based on the score (𝐼 𝑛,𝑠𝑐𝑜 ) and of the observed Fisher information matrix (𝐼 𝑛,𝑜𝑏𝑠 ) computed in the true parameter values from the 500 simulated datasets with n=500.

Figure 2 :

 2 Figure 2: Poisson mixture model. Kernel density estimates of the normalized values of some components of the estimated Fisher information matrix based on the score (𝐼 𝑛,𝑠𝑐𝑜 ) and of the observed Fisher information matrix (𝐼 𝑛,𝑜𝑏𝑠 ) computed in the true parameter values from the 500 simulated datasets with n=500.

𝑘

  and 𝒦 𝜅 = [-20-𝜅, 20+𝜅] 6 ×[0, 5.10 4 +𝜅].

  :nbiterem],c(1,2,4),mean) MbiasIobs <-apply(biasIobs[,,,(nbiterburnin+1):nbiterem],c(1,2,4),mean) MsdIsco <-apply(rsdIsco[,,,(nbiterburnin+1):nbiterem],c(1,2,4),mean) MsdIobs <-apply(rsdIobs[,,,(nbiterburnin+1):nbiterem],c(1,2,4),mean) save(MbiasIsco,file='Rfiles/ResNLMEexponentialBiasIsco.Rdata') save(MbiasIobs,file='Rfiles/ResNLMEexponentialBiasIobs.Rdata') save(MsdIsco,file='Rfiles/ResNLMEexponentialSdIsco.Rdata') save(MsdIobs,file='Rfiles/ResNLMEexponentialSdIobs.Rdata') load('Rfiles/ResNLMEexponentialBiasIsco.Rdata') load('Rfiles/ResNLMEexponentialBiasIobs.Rdata') load('Rfiles/ResNLMEexponentialSdIsco.Rdata') load('Rfiles/ResNLMEexponentialSdIobs.Rdata') MbiasIobs <-MbiasIobs[,,seq(1,2000,10)] MbiasIsco <-MbiasIsco[,,seq(1,2000,10)] MsdIobs <-MsdIobs[,,seq(1,2000,10)] MsdIsco <-MsdIsco[,,seq(1,2000,10)] save(MbiasIsco,file='Rfiles/ResNLMEexponentialBiasIsco.Rdata') save(MbiasIobs,file='Rfiles/ResNLMEexponentialBiasIobs.Rdata') save(MsdIsco,file='Rfiles/ResNLMEexponentialSdIsco.Rdata') save(MsdIobs,file='Rfiles/ResNLMEexponentialSdIobs.Rdata')

Figure 5 :

 5 Figure 5: Non linear mixed effects model. Convergence plot for some parameter estimates and for some diagonal components of 𝐼 𝑛,𝑠𝑐𝑜 ( θ ) over iterations of the stochastic approximation algorithm.

  SAEM2) 𝑙 ∶ Θ → ℝ and θ ∶ 𝒮 → Θ are 𝑚 times differentiable, where 𝑚 is the integer such that 𝒮 is an open subset of ℝ 𝑚 . • (SAEM3) For all positive Borel functions Φ, we have 𝐸[Φ(𝑧 𝑘+1 )|ℱ 𝑘 ] = ∫ Φ(𝑧)𝑝(𝑧; 𝜃 𝑘 )𝜇(𝑑𝑧). • (SAEM4) For all 𝜃 ∈ Θ, E 𝜃 (‖𝑆(𝑍 )‖ 2 ) < ∞, and the function Γ(𝜃) ≜ Cov 𝜃 [𝑆(𝑧)] ≜ ∫ 𝑆(𝑧) 𝑡 𝑆(𝑧)𝑝(𝑧; 𝜃)𝜇(𝑑𝑧) -[ ∫ 𝑆(𝑧)𝑝(𝑧; 𝜃)𝜇(𝑑𝑧)] 𝑡 [ ∫ 𝑆(𝑧)𝑝(𝑧; 𝜃)𝜇(𝑑𝑧)] is continuous w.r.t. 𝜃, where E 𝜃 stands for the expectation with respect to the posterior distribution 𝑝(⋅; 𝜃).

  LOC2) For every stationary point 𝜃 * , the matrices E * 𝜃 (𝜕 𝜃 𝐿(𝑆(𝑍 ), 𝜃 * )(𝜕 𝜃 𝐿(𝑆(𝑍 ), 𝜃 * )) 𝑡 ) and 𝜕 2 𝜃 𝐿(E * 𝜃 (𝑆(𝑍 )), 𝜃 * ) are positive definite. • (LOC3) The minimum eigenvalue of the covariance matrix 𝑅(𝜃) = E 𝜃 ((𝑆(𝑍 )-s (𝜃))(𝑆(𝑍 )-s (𝜃)) 𝑡 ) is bounded away from zero for 𝜃 in any compact subset of Θ. computation of the Fisher information matrix in the linear mixed effects model ## Computation of the exact Fisher Information matrix in the linear mixed-effects ## model Fisher_LMM <-function(beta,sigma2,eta2,j){ # beta : value of the fixed-effect # sigma2 : value of the residual variance # eta2

  information matrix extimation in the linear mixed effects model ## Computation of the estimator of the Fisher information matrix based on ## the score Isco_LMM <-function(datamat,beta,sigma2,eta2){ # datamat : observations organized in a matrix. Each row is an individual # vector of observations. information matrix estimation in the Poisson mixture model ## Function for computing Isco and Iobs for Fisher Information matrix estimation ## in Poisson mixture models fisher_estimation_poisson_mixture <-function(y, lambda, alpha) { # y : vector of observations # lambda : vector of K Poisson parameters for each component of the mixture #(in ascending order) # alpha : vector of (K-1) mixture proportions K <-length(lambda) # number of components of the mixture n <-length(y) # sample sizederiv1ind <-matrix(0,2*K-1,n) deriv2 <-matrix(0,2*K-1,2*K-1) covderiv <-matrix(0,2*K-1,2*K-1) ## computationof conditional expectation of the first derivatives of the ## complete data log-likelihood denom <-0 for (k in 1:(K-1)){ denom <-denom + exp(-lambda[k])*lambda[k]^y*alpha[k] } denom <-denom + exp(-lambda[K])*lambda[K]^y*(1-sum(alpha)) for (k in 1:

  ,1] <-t(log(psi)^2) candidatka[k] <-candidatka[k] + rnorm(1,0,sqrt(eta2ka)) psicandidat <-cbind(exp(candidatka),exp(currentv),exp(currentcl)) logs <--1/2/sigma2*sum((y-model1cpt(psicandidat,id,xidep))^2)+ 1/2/sigma2*sum((y-model1cpt(currentpsi,id,xidep))^2) logs <-logs-1/2/omega2ka*((candidatka[k]-log(kapop))^2-(currentka[k]-log(kapop)candidatv[k] <-candidatv[k] + rnorm(1,0,sqrt(eta2v)) psicandidat <-cbind(exp(currentka),exp(candidatv),exp(currentcl)) logs <--1/2/sigma2*sum((y-model1cpt(psicandidat,id,xidep))^2)+ 1/2/sigma2*sum((y-model1cpt(currentpsi,id,xidep))^2) logs <-logs -1/2/omega2v*((candidatv[k]-log(vpop))^2-(currentv[k]-log(vpop)candidatcl[k] <-candidatcl[k] + rnorm(1,0,sqrt(eta2cl)) psicandidat<-cbind(exp(currentka),exp(currentv),exp(candidatcl)) of iterations of the SAEM algorithm # nbiterburnin : number of burn-in iterations of the algorithm # theta0: initial parameter values # kRW : coefficient used to adjust the variance of the proposal kernel # of the MCMC procedure# Q quantity floglik <-function(v,y,psi,xidep,id,alpha){ l <-length(alpha) psi <-array(psi,dim=c(n,2,l)) value <-0 for (ll in 1:l){ moyij <-model1cptV(psi[,,ll],id,xidep,v) value <-value + alpha[ll]*sum((y-moyij)<-c(kapop, vpop, clpop, omega2ka, omega2cl, sigma2)# variances of the proposal kernels of the MCMC procedure eta2ka <-kRW*omega2ka eta2cl <-kRW*omega2cl# sequence of step sizes gamma <-1/(1:(nbiterem))^0.501 cumgamma <-matrix(0,nbiterem,nbiterem)diag(cumgamma) <-gamma for (l in 2:nbiterem){ for (m in 1:(l-1)){ cumgamma[l,m] <-(1-gamma[l])*cumgamma[l-1(NA,dim=c(nbiterem,n,j)) mco3 <-array(NA,dim=c(nbiterem,n,j)) mcos <-rep(NA,n) mcos2 <-matrix(NA,nbiterem,n) STATEXH <-rep(NA,dimstatexh) psisauv <-array(0,c(n,2,nbiterem)) H <-array(0,c(p,p,n,nbiterem)) G2<-array(0,c(p,p,nbiterem)) tempderiveeas2 <-matrix(0,p,p) # initial values for the individual parameters currentka <-log(kapop) + rnorm(n,0,sqrt(eta2ka)) currentcl <-log(clpop) + rnorm(n,0,sqrt(eta2cl)) currentpsi <-cbind(exp(currentka),exp(currentcl))## Start of the EM loop for (l in 1:(nbiterem-1)){ ## Simulation step for (k in 1:(n)){ # Parameter ka candidatka <-currentka candidatka[k] <-candidatka[k] + rnorm(1,0,sqrt(eta2ka)) psicandidat <-cbind(exp(candidatka),exp(currentcl)) logs <--1/2/sigma2*sum((y-model1cptV(psicandidat,id,xidep,vpop))^2)+ 1/2/sigma2*sum((y-model1cptV(currentpsi,id,xidep,vpop))^2) logs <-logs-1/2/omega2ka*((candidatka[k]-log(kapop))^2-(currentka[k]-log(kapop))^2) u <-runif(1) logu <-log(u) ind <-(logu<logs) currentpsi <-psicandidat*ind + currentpsi*(1-ind) currentka <-candidatka*ind + currentka*(1-ind) # Parameter cl candidatcl <-currentcl candidatcl[k] <-candidatcl[k] + rnorm(1,0,sqrt(eta2cl)) psicandidat<-cbind(exp(currentka),exp(candidatcl)) logs <--1/2/sigma2*sum((y-model1cptV(psicandidat,id,xidep,vpop))^2)+ 1/2/sigma2*sum((y-model1cptV(currentpsi,id,xidep,vpop))^2) logs <-logs -1/2/omega2cl*((candidatcl[k]-log(clpop))^2-(currentcl[k]-log(clpop)*ind + currentpsi*(1-ind) currentcl <-candidatcl*ind + currentcl*(1-ind) # saving simulated data psisauv[,,l] <-currentpsi } psi <-psisauv[,,l] ## Parameter estimation update # estimation of the fixed effect by numerical optimization resvpop <-optimize(interval=c(0.001,50),f=floglik,y=y,psi=psisauv[,,1:l], xidep=xidep,id=id,alpha=cumgamma[l,1:l]) vpop <-resvpop$minimum # stochastic approximation of exhaustive statistics and estimation of the # other parameters mco <-matrix((y-model1cptV(psi,id,xidep,vpop))^2,n,j,byrow=TRUE) mcos <-apply(mco,1,sum) STATEXH <-c(apply(log(psi),2,mean), apply(log(psi)^2,2,mean), sum(mcos)) statexh[,l+1] <-statexh[,l]*(1-gamma[l])+gamma[l]

  ℝ 𝑑 𝑖 taking its values in a subset 𝒮 𝑖 of ℝ 𝑚 𝑖 .Let us denote for all 1 ≤ 𝑖 ≤ 𝑛 by 𝐿 𝑖 the function defined on 𝒮 𝑖 × Θ by 𝐿 𝑖 (𝑠 𝑖 ; 𝜃) ≜ -𝜓 𝑖 (𝜃) + ⟨𝑠 𝑖 , 𝜙 𝑖 (𝜃)⟩ and by 𝐿 ∶ 𝒮 × Θ → ℝ the function defined as 𝐿(𝑠, 𝜃) = ∑ 𝑖 𝐿 𝑖 (𝑠 𝑖 ; 𝜃) with 𝒮 = ∏ 𝑖 𝒮 𝑖 and 𝑠 = (𝑠 1 , … , 𝑠 𝑛 ). For sake of simplicity, we omitted all dependency on the observations (𝑦 𝑖 ) 1≤𝑖≤𝑛 since the considered stochasticity relies here on the latent variables.Finally let us denote by (𝛾 𝑘 ) 𝑘≥1 and (𝜀 𝑘 ) 𝑘≥1 sequences of positive step sizes, by K a compact set of ℝ 𝑑 with 𝑑 = ∑ 𝑑 𝑖 and by (𝒦 𝑘 ) a sequence of increasing compact sets of 𝒮 such that ∪𝒦 𝑘 = 𝒮 and for all 𝑘 𝒦 𝑘 ⊂ 𝑖𝑛𝑡(𝒦 𝑘+1 ).

  -R script for studying the asymptotic properties of Iobs and Isco in the ## Poisson mixture model -Monte-Carlo estimation of the true Fisher information ## matrix based on a very large sample ## -

	## 2
		j)
	obs	<-beta+randompop+residual
	datamat	<-matrix(obs,n,j)
	## evaluation of the FIM estimators in the true parameter values
	resIobs.theta.true[,,k,l] <-Iobs_LMM(datamat,beta,sigma2,eta2)
	resIsco.theta.true[,,k,l] <-Isco_LMM(datamat,beta,sigma2,eta2)
	## evaluation of the FIM estimators in the estimated parameter values
	est.mle	<-lmer(obs~(1|id),REML = F)
	variances <-as.data.frame(VarCorr(est.mle))
	beta.est.n	<-c(beta.est.n,est.mle@beta)
	eta2.est.n	<-c(eta2.est.n,variances[1,'vcov'])
	sigma2.est.n <-c(sigma2.est.n,variances[2,'vcov'])
	resIobs.theta.est[,,k,l] <-Iobs_LMM(datamat,est.mle@beta,
		variances[2,'vcov'],
		variances[1,'vcov'])
	resIsco.theta.est[,,k,l] <-Isco_LMM(datamat,est.mle@beta,
		variances[2,'vcov'],
		variances[1,'vcov'])

Table 1 :

 1 Linear mixed effects model. Empirical bias to the Fisher Information matrix of 𝐼 𝑛,𝑠𝑐𝑜 and 𝐼 𝑛,𝑜𝑏𝑠 computed in the true parameter values for different values of n.

		n=20		n=100	n=500
		Isco	Iobs	Isco	Iobs	Isco	Iobs
	(β,β)	0.00552	0.00000	0.00075	0.00000	-0.00064	0.00000
	(η2,η2)	0.00325	0.00228	-0.00138	0.00031	-0.00122	-0.00027
	(σ2,σ2)	-0.00314	0.00009	0.00002	0.00077	-0.00027	-0.00013
	(β,η2)	-0.00241	0.00255	-0.00040	-0.00062	-0.00058	-0.00133
	(β,σ2)	0.00360	0.00021	-0.00006	-0.00005	-0.00043	-0.00011
	(η2,σ2)	-0.00050	0.00019	0.00013	0.00003	0.00002	-0.00002

Table 2 :

 2 Linear mixed effects model. Empirical squared deviation to the Fisher Information matrix of 𝐼 𝑛,𝑠𝑐𝑜 and 𝐼 𝑛,𝑜𝑏𝑠 computed in the true parameter values for different values of n.

		n=20		n=100	n=500
		Isco	Iobs	Isco	Iobs	Isco	Iobs
	(β,β)	0.13433	0.00000	0.05843	0.00000	0.02634	0.00000
	(η2,η2)	0.07916	0.05559	0.02931	0.02418	0.01372	0.01090
	(σ2,σ2)	0.08201	0.04177	0.03982	0.01946	0.01729	0.00819
	(β,η2)	0.09730	0.05970	0.04380	0.02786	0.01898	0.01145
	(β,σ2)	0.07115	0.00497	0.03100	0.00232	0.01378	0.00095
	(η2,σ2)	0.02907	0.00463	0.01308	0.00201	0.00639	0.00091

Table 3 :

 3 Poisson mixture model. Empirical bias to the Fisher Information matrix of 𝐼 𝑛,𝑠𝑐𝑜 and 𝐼 𝑛,𝑜𝑏𝑠 computed in the true parameter values for different values of n.

		n=20	n=100	n=500
		Isco	Iobs	Isco	Iobs	Isco	Iobs
	(λ2,λ2)	0.00009	-0.00025	-0.00003	0.00025	0.00008	-0.00029
	(λ3,λ3)	0.00005	0.00047	-0.00023	-0.00035	0.00008	0.00015
	(α1,α1)	0.05981	0.05981	-0.04072	-0.04072	0.01849	0.01849
	(α2,α2)	0.04756	0.04756	-0.04006	-0.04006	0.01205	0.01205
	(λ2,λ3)	0.00009	0.00009	-0.00007	-0.00007	0.00002	0.00002
	(λ3,α2)	-0.00220	0.00082	0.00284	-0.00061	-0.00077	0.00041

Table 4 :

 4 Poisson mixture model. Empirical squared deviation to the Fisher Information matrix of 𝐼 𝑛,𝑠𝑐𝑜 and 𝐼 𝑛,𝑜𝑏𝑠 computed in the true parameter values for different values of n.

		n=20		n=100	n=500
		Isco	Iobs	Isco	Iobs	Isco	Iobs
	(λ2,λ2)	0.00717	0.02238	0.00310	0.00996	0.00141	0.00463
	(λ3,λ3)	0.01523	0.00872	0.00664	0.00403	0.00299	0.00167
	(α1,α1)	1.20192	1.20192	0.52483	0.52483	0.23129	0.23129
	(α2,α2)	1.05566	1.05566	0.46762	0.46762	0.20510	0.20510
	(λ2,λ3)	0.00295	0.00295	0.00132	0.00132	0.00059	0.00059
	(λ3,α2)	0.11013	0.03428	0.04614	0.01561	0.02137	0.00712

Table 5 :

 5 Linear mixed effects model. Comparison of the coverage rates computed from both estimates of the Fisher information matrix in either the true or the estimated parameter values when n=100.

	1-α	Fisher est. θ	β	η2	σ2
		Isco	Known	0.992	0.98	0.984
	0.99		Estimated 0.988	0.974	0.992
		Iobs	Known	0.99	0.982	0.982
			Estimated 0.986	0.972	0.988
		Isco	Known	0.964	0.946	0.948
	0.95		Estimated 0.944	0.936	0.944
		Iobs	Known	0.96	0.944	0.948
			Estimated 0.938	0.942	0.944
		Isco	Known	0.904	0.914	0.876
	0.9		Estimated 0.89	0.906	0.884
		Iobs	Known	0.892	0.916	0.89
			Estimated 0.882	0.89	0.872

Table 6 :

 6 Linear mixed effects model. Comparison of the coverage rates computed from both estimates of the Fisher information matrix in either the true or the estimated parameter values when n=500.

	1-α	Fisher est. θ	β	η2	σ2
		Isco	Known	0.996	0.986	0.984
	0.99		Estimated 0.996	0.986	0.992
		Iobs	Known	0.996	0.986	0.984
			Estimated 0.996	0.988	0.986
		Isco	Known	0.952	0.954	0.952
	0.95		Estimated 0.95	0.956	0.944
		Iobs	Known	0.952	0.956	0.952
			Estimated 0.95	0.95	0.946
		Isco	Known	0.922	0.914	0.906
	0.9		Estimated 0.916	0.898	0.91
		Iobs	Known	0.918	0.912	0.908
			Estimated 0.914	0.888	0.912

  𝑧 𝑖 = (log 𝑘𝑎 𝑖 , log 𝐶𝑙 𝑖 , log 𝑉 𝑖 ) ′ are individual random parameters such thatlog 𝑘𝑎 𝑖 = log(𝑘𝑎) + 𝜂 𝑖,1 , log 𝐶𝑙 𝑖 = log(𝐶𝑙) + 𝜂 𝑖,2 , log 𝑉 𝑖 = log(𝑉 ) + 𝜂 𝑖,3 ,and𝑔 𝑖 (𝑡 𝑖𝑗 , 𝑧 𝑖 ) = 𝑑 𝑖 𝑘𝑎 𝑖 𝑉 𝑖 𝑘𝑎 𝑖 -𝐶𝑙 𝑖 [𝑒 𝑡 𝑖𝑗 -𝑒 -𝑘𝑎 𝑖 𝑡 𝑖𝑗 ] .For all 1 ≤ 𝑖 ≤ 𝑛 and all 1 ≤ 𝑗 ≤ 𝐽, 𝑦 𝑖𝑗 denotes the measure of drug concentration on individual 𝑖 at time 𝑡 𝑖𝑗 , 𝑑 𝑖 the dose of drug administered to individual 𝑖, and 𝑉 𝑖 , 𝑘𝑎 𝑖 and 𝐶𝑙 the residuals (𝜀 𝑖𝑗 ) are assumed independent and identically distributed such that 𝜀 𝑖𝑗 ∼ and the sequences (𝜂 𝑖 ) and (𝜀 𝑖𝑗 ) are assumed mutually independent. Here, the model parameter is 𝜃 = (𝑘𝑎, 𝑉 , 𝐶𝑙, 𝜔 2 𝑘𝑎 , 𝜔 2 𝑉 , 𝜔 2 𝐶𝑙 , 𝜎 2

		𝑖.𝑖.𝑑.	𝒩 (0, 𝜎 2 )
	-	𝐶𝑙 𝑖 𝑉 𝑖

𝑖 respectively denote the volume of the central compartment, the drug's absorption rate constant and the drug's clearance of individual 𝑖. The terms 𝜂 𝑖 = (𝜂 𝑖,1 , 𝜂 𝑖,2 , 𝜂 𝑖,3 ) ′ ∈ ℝ 3 are unobserved random effects which are assumed independent and identically distributed such that 𝜂 𝑖 ∼ 𝑖.𝑖.𝑑. 𝒩 (0, Ω), where Ω = diag(𝜔 2 𝑘𝑎 , 𝜔 2 𝐶𝑙 , 𝜔 2 𝑉 ),

  Figure 3: Non linear mixed effects model. Representation over iterations of the mean relative bias of the diagonal components of the estimated Fisher information matrix computed from the 𝑀 = 500 runs of the stochastic algorithm. Red line corresponds to 𝐼 𝑛,𝑠𝑐𝑜 (𝜃) and blue line corresponds to 𝐼 𝑛,𝑜𝑏𝑠 (𝜃).The burn-in iterations of the algorithm are not depicted.

		( ka , ka )			( V , V )					( Cl , Cl )
	0.20						0.06						0.15
	0.15						0.04						0.10
	0.10						0.02						0.05
	0.05	1000	1500	2000	2500	3000	0.00	1000	1500	2000	2500	3000	0.00	1000	1500	2000	2500	3000
	( ω ka 2 , ω ka 2 )			( ω V 2 , ω V 2 )				( ω Cl 2 , ω Cl 2 )
	0.8						0.20						0.4
	0.6						0.15						0.3
	0.4						0.10						0.2
	0.2						0.05						0.1
		1000	1500	2000	2500	3000	0.00	1000	1500	2000	2500	3000	0.0	1000	1500	2000	2500	3000
		( ka , ka )				( V , V )					( Cl , Cl )
	0.08						0.008						0.04
	0.06						0.006						0.03
	0.04						0.004						0.02
	0.02						0.002						0.01
	0.00												0.00
		1000	1500	2000	2500	3000		1000	1500	2000	2500	3000		1000	1500	2000	2500	3000
	( ω ka 2 , ω ka 2 )			( ω V 2 , ω V 2 )
	0.8						0.06
	0.2 0.4 0.6						0.02 0.04						0.1 0.2
	0.0												0.0
		1000	1500	2000	2500	3000		1000	1500	2000	2500	3000		1000	1500	2000	2500	3000

  Theorem 6.1. Assume that (𝑀1 ′ ) and (𝑀2 ′ ), (𝑀3) to (𝑀5) and (𝑆𝐴𝐸𝑀1) to (𝑆𝐴𝐸𝑀4) are fulfilled. Assume also that with probability 1 clos({𝑠 𝑘 } 𝑘≥1 ) is a compact subset of 𝒮. Let us define ℒ = {𝜃 ∈ Θ, 𝜕 𝜃 𝑙(𝑦; 𝜃) = 0} the set of stationary points of the observed log-likelihood 𝑙 defined as 𝑙(𝑦; 𝜃) = ∑ 𝑛 𝑖=1 log 𝑔(𝑦 𝑖 ; 𝜃). Then, for all 𝜃 0 ∈ Θ, for fixed 𝑛 ∈ ℕ * , we get: lim 𝑘 𝑑(𝜃 𝑘 , ℒ ) = 0 and lim 𝑘 𝑑(𝐼 𝑘 𝑛,𝑠𝑐𝑜 , ℐ ) = 0 a.s. where ℐ = {𝐼 𝑛,𝑠𝑐𝑜 (𝜃), 𝜃 ∈ ℒ }.Proof. Let us denote by 𝑆(𝑍 ) = (𝑆 1 (𝑍 1 ), … , 𝑆 𝑛 (𝑍 𝑛 )) the sufficient statistics of the model we consider in our approach. Note as recalled in(Delyon B., Lavielle M., and Moulines E. 1999), these are not unique.Let us also define 𝐻 (𝑍 , 𝑠) = 𝑆(𝑍 ) -𝑠 and ℎ(𝑠) = E 𝑍 |𝑌 ; θ (𝑠) (𝑆(𝑍 )) -𝑠. Assumptions (𝑀1 ′ ) and (𝑀2 ′ ) imply that assumptions (𝑀1) and (𝑀2) of Theorem 5 of (Delyon B., Lavielle M., and Moulines E. 1999) are fulfilled. Indeed these assumptions focus on expressions and regularity properties of the individual likelihood functions and the corresponding sufficient statistics for each index 𝑖 ∈ {1, … , 𝑛}. Then by linearity of the log-likelihood function and of the stochastic approximation and applying Theorem 5 of (Delyon B., Lavielle M., and Moulines E. 1999), we get that lim 𝑘 𝑑(𝜃 𝑘 , ℒ ) = 0. Moreover we get that for 1 ≤ 𝑖 ≤ 𝑛, each sequence (𝑠 𝑘 𝑖 ) converges almost surely toward E 𝑍 𝑖 |𝑌 𝑖 ;𝜃 (𝑆 𝑖 (𝑍 𝑖 )). Since assumption (𝑀2 ′ ) ensures that for all 1 ≤ 𝑖 ≤ 𝑛 the functions 𝜓 𝑖 and 𝜙 𝑖 are twice continuously differentiable and assumption (𝑀5) ensures that the function θ is continuously differentiable, the function Φ 𝑛 defined by Φ 𝑛 (𝑠 𝑘 ) = 1 𝑛 ∑ 𝑛 𝑖=1 Δ𝑖 (𝑠 𝑘 ) Δ𝑖 (𝑠 𝑘 ) is continuous. Therefore we get that lim 𝑘 𝑑(𝐼 𝑘 𝑛,𝑠𝑐𝑜 , ℐ ) = 0.

	We now establish the asymptotic normality of the estimate Ī 𝑘 𝑛,𝑠𝑐𝑜 defined as Ī 𝑘 𝑛,𝑠𝑐𝑜 = Φ 𝑛 ( s 𝑘 ) with
	s 𝑘 = ∑	𝑘 𝑙=1

  The parameter space Θ is an open subset of ℝ 𝑝 . The individual complete data likelihood function is given for all 𝑖 = 1, … , 𝑛 by:𝑓 𝑖 (𝑧 𝑖 ; 𝜃) = exp (-𝜓 𝑖 (𝜃) + ⟨𝑆 𝑖 (𝑧 𝑖 ), 𝜙 𝑖 (𝜃)⟩) ,where ⟨⋅, ⋅⟩ denotes the scalar product, 𝑆 𝑖 is a Borel function on ℝ 𝑑 𝑖 taking its values in an open subset 𝒮 𝑖 of ℝ 𝑑 𝑖 , 𝜙 𝑖 and 𝜓 𝑖 are measurable function of Θ taking values in open subsets of ℝ 𝑑 𝑖 and ℝ respectively. Moreover, the convex hull of 𝑆(ℝ ∑ 𝑑 For all 1 ≤ 𝑖 ≤ 𝑛 the function 𝑙 𝑖 ∶ Θ → ℝ defined as 𝑙 𝑖 (𝜃) = log ∫ 𝑓 𝑖 (𝑧 𝑖 ; 𝜃)𝜇 𝑖 (𝑑𝑧 𝑖 ) is continuously differentiable on Θ and 𝜕 𝜃 ∫ 𝑓 𝑖 (𝑧 𝑖 ; 𝜃)𝜇 𝑖 (𝑑𝑧 𝑖 ) = ∫ 𝜕 𝜃 𝑓 𝑖 (𝑧 𝑖 ; 𝜃)𝜇 𝑖 (𝑑𝑧 𝑖 ).

𝑖

) is included in 𝒮 and for all 𝜃 ∈ Θ ∫ 𝑆(𝑧) ∏ 𝑝 𝑖 (𝑧 𝑖 ; 𝜃)𝜇(𝑑𝑧) < ∞

• (M2') Define for each 𝑖 𝐿 𝑖 ∶ 𝒮 𝑖 × Θ → ℝ as 𝐿 𝑖 (𝑠 𝑖 ; 𝜃) ≜ -𝜓 𝑖 (𝜃) + ⟨𝑠 𝑖 , 𝜙 𝑖 (𝜃)⟩.The functions 𝜓 𝑖 and 𝜙 𝑖 are twice continuously differentiable on Θ.

• (M3)

The function s ∶ Θ → 𝒮 defined as s (𝜃) ≜ ∫ 𝑆(𝑧)𝑝(𝑧; 𝜃)𝜇(𝑑𝑧) is continuously differentiable on Θ.

• (M4)

4 SAEM algorithm in the PK model belonging to the curved exponential family #

  # R function implementing the saem algorithm to compute the parameter estimates ## and the FIM estimates simultaneously in the PK nonlinear mixed-effects model ## with three random effects, thus belonging to the curved exponential family

	## computation of Isco and Iobs
	Isco <-deriv1ind%*%t(deriv1ind)/n
	Iobs <-deriv1ind%*%t(deriv1ind)/n -deriv2/n -covderiv/n
	res <-list(Isco = Isco, Iobs = Iobs)
	return(res)
	}
	6.3.saem <-function(data, nbiterem, nbiterburnin, theta0, kRW=0.5) {
	# data	: dataset
	# nbiterem	: total number of iterations of the saem algorithm
	# nbiterburnin : number of burn-in iterations of the algorithm
	# theta0	: initial parameter values
	# kRW		: coefficient used to adjust the variance of the proposal
	#		kernel of the MCMC procedure
	# data processing
	xidep <-cbind(data$dose,data$time)
	y	<-data$y
	id	<-as.matrix(data$subject)
	n	<-length(unique(id))
	j	<-length(unique(data$time))
	nb.psi	<-3
	dimstatexh	<-7
	# initial parameter values
	vpop
			lambda[k]^y/denom/alpha[k] -

# beta

: value of the fixed-effect # sigma2 : value of the residual variance # eta2

: value of the random effects variance n <-dim(datamat)[1] j <-dim(datamat) [2] derivative <-matrix(0,3,n) derivative[1,] <-apply(datamat-beta,1,sum)/(sigma2+eta2*j) derivative[2,] <-apply(datamat-beta,1,sum)^2/2/(sigma2+eta2*j)^2j/2/(sigma2+eta2*j) derivative[3,] <-apply((datamat-beta)^2,1,sum)/2/sigma2^2apply(datamat-beta,1,sum)^2*eta2*(j*eta2+2*sigma2)/2/sigma2^2/(sigma2+eta2*j)^2 -1/2/(sigma2+eta2*j)-(j-1)/2/sigma2

Isco <-derivative%*%t(derivative)/n return(Isco) } ## Computation of the observed information matrix in the linear mixed-effects ## model Iobs_LMM <-function(datamat,beta,sigma2,eta2){ # datamat : observations organized in a matrix. Each row is an individual # vector of observations. # beta : value of the fixed-effect # sigma2 : value of the residual variance # eta2

: value of the random effects variance

## computation of the conditional covariance matrix of the first derivatives ## of the complete data log-likelihood

## Computation of the FIM estimations isco <-array(0,c(p,p,nbiterem)) iobs <-array(0,c(p,p,nbiterem)) SH <-vector("list",nbiterem)

]<=(50000))) res <-as.numeric(res) return(res) }

SAEM algorithm in the PK model not belonging to the curved exponential family

## R function implementing the saem algorithm to compute the parameter estimates ## and the FIM estimates simultaneously in the PK nonlinear mixed-effects model ## with two random effects, thus not belonging to the curved exponential family saem_non_exp <-function (data,nbiterem,nbiterburnin,theta0,kRW=0.5 Kuhn E., and Trouvé A. 2010