Thu Thuy Nguyen
email: thu-thuy.nguyen@inserm.fr

France Mentré

Evaluation of the Fisher information matrix in nonlinear mixed effect models using adaptive Gaussian quadrature

Keywords: Design, Dose-response studies, Fisher information matrix, adaptive Gaussian quadrature, Linearisation, Nonlinear mixed effect model

Nonlinear mixed effect models (NLMEM) are used in model-based drug development to analyse longitudinal data. To design these studies, the use of the expected Fisher information matrix (M F ) is a good alternative to clinical trial simulation. Presently, M F in NLMEM is mostly evaluated with first-order linearisation. The adequacy of this approximation is, however, influenced by model nonlinearity. Alternatives for the evaluation of M F without linearisation are proposed, based on Gaussian quadratures. The M F , expressed as the expectation of the derivatives of the log-likelihood, can be obtained by stochastic integration. The likelihood for each simulated vector of observations is approximated by Gaussian quadrature centred at 0 (standard quadrature) or at the simulated random effects (adaptive quadrature). These approaches have been implemented in R. Their relevance was compared with clinical trial simulation and linearisation, using dose-response models, with various nonlinearity levels and different number of doses per patient. When the nonlinearity was mild, three approaches based on M F gave correct predictions of standard errors, when compared with the simulation. When the nonlinearity increased, linearisation correctly predicted standard errors of fixed effects, but over-predicted, with sparse designs, standard errors of some variability terms. Meanwhile, quadrature approaches gave correct predictions of standard errors overall, but standard Gaussian quadrature was very time-consuming when there were more than two random effects. To conclude, adaptive Gaussian quadrature is a relevant alternative for the evaluation of M F for models with stronger nonlinearity, while being more computationally efficient than standard quadrature.

Introduction

Nonlinear mixed effect models (NLMEM) are frequently used in model-based drug development to analyse longitudinal data obtained during clinical trials [START_REF] Lalonde | Model-based drug development[END_REF][START_REF] Smith | Biostatistics and pharmacometrics: quantitative sciences to propel drug development forward[END_REF]. They were introduced about 40 years ago [START_REF] Sheiner | Modelling of individual pharmacokinetics for computer-aided drug dosage[END_REF][START_REF] Sheiner | Estimation of population characteristics of pharmacokinetic parameters from routine clinical data[END_REF] and were initially used in pharmacokinetic analyses as an alternative to the non-compartmental approach (NCA) [START_REF] Gabrielson | Pharmacokinetic and Pharmacodynamic Data Analysis: Concepts and Applications[END_REF]. This modelling approach is more complex than NCA, but allows for the analysis of few samples per subject. It accounts for within and between subject variability, and is appropriate for exploiting the richness of repeated measurements. Consequently, this approach is increasingly used in the biomedical field, not only for pharmacokinetic analyses [START_REF] Sheiner | Modelling of individual pharmacokinetics for computer-aided drug dosage[END_REF][START_REF] Sheiner | Estimation of population characteristics of pharmacokinetic parameters from routine clinical data[END_REF], but also for analyses of viral loads [START_REF] Perelson | Estimating drug efficacy and viral dynamic parameters: HIV and HCV[END_REF], of bacterial resistance to antibiotics [START_REF] Nielsen | Semimechanistic pharmacokinetic/pharmacodynamic model for assessment of activity of antibacterial agents from time-kill curve experiments[END_REF], and of the dose-response relationship.

This approach has become the main statistical tool in pharmacometrics, the science of quantitative pharmacology [START_REF] Van Der Graaf | CPT: Pharmacometrics and Systems Pharmacology[END_REF]. Parameters of these models are commonly estimated by likelihood maximisation [START_REF] Dartois | Overview of model-building strategies in population PK/PD analyses: 2002-2004 literature survey[END_REF]. However, the nonlinearity of the structural model prevents a closed form solution for the integration over the random effects in the expression of the likelihood function. Many approaches have been proposed over the years to overcome this difficulty, and implemented in several estimation software packages. These are first-order marginal quasi-likelihood or first-order linearisation [START_REF] Lindstrom | Nonlinear mixed effects models for repeated measures data[END_REF] in NONMEM, R, Splus, Laplace approximation [START_REF] Wolfinger | Laplace's approximation for nonlinear mixed models[END_REF] in NONMEM and SAS, adaptive Gaussian quadrature [START_REF] Pinheiro | Approximations to the log-likelihood function in the nonlinear mixed-effects model[END_REF] in SAS and in the R package lme4, Stochastic Approximation Expectation Maximisation (SAEM) [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF] in MONOLIX and NONMEM. [START_REF] Pillai | Non-linear mixed effects modeling -from methodology and software development to driving implementation in drug development science[END_REF] described these estimation methods in a review paper and recently, [START_REF] Plan | Performance comparison of various maximum likelihood nonlinear mixed-effects estimation methods for dose-response models[END_REF] compared their performance, showing that adaptive Gaussian quadrature, although the slowest, was generally the best method.

Before the modelling step to estimate parameters, it is important to define an appropriate design, which consists in determining a balance between the number of subjects and the number of samples per subject, as well as the allocation of times and doses, according to experimental conditions. The choice of design is crucial for an efficient estimation of model parameters [START_REF] Al-Banna | Experimental design and efficient parameter estimation in population pharmacokinetics[END_REF][START_REF] Hashimoto | Designs for population pharmacodynamics: value of pharmacokinetic data and population analysis[END_REF][START_REF] Jonsson | Comparison of some practical sampling strategies for population pharmacokinetic studies[END_REF], especially when the studies are conducted in children or in patients where only a few samples can be taken per subject. The main approach for design evaluation has long been based on clinical trial simulation (CTS), but it is a cumbersome method and so the number of designs that can be evaluated is limited. An alternative approach has been described in the general theory of optimum experimental design used for classical nonlinear models [START_REF] Atkinson | Optimum Experimental Designs, with SAS[END_REF][START_REF] Walter | Identification of Parametric Models from Experimental Data[END_REF][START_REF] Atkinson | Elemental information matrices and optimal experimental design for generalized regression models[END_REF], relying on the Rao-Cramer inequality which states that the inverse of the Fisher information matrix (M F ) is the lower bound of the variance-covariance matrix of any unbiased estimate of the parameters and its diagonal elements are the expected variances of the parameters. Several criteria based on M F have been developed to evaluate designs. One of the criteria widely used is the criterion of D-optimality, which consists in maximising the determinant of M F . The computation of this criterion requires a priori knowledge of the model and its parameters, which can usually be obtained from previous experiments. This leads to the concept of "locally optimal designs", which has been studied in several publications [START_REF] Chernoff | Locally optimal designs for estimating parameters[END_REF][START_REF] Box | Design of experiments in nonlinear situations[END_REF][START_REF] D'argenio | Optimal sampling times for pharmacokinetic experiments[END_REF]. Since there is no closed form of the likelihood in NLMEM, there is no analytical expression of M F . An approximation of the expected M F has been proposed for NLMEM, using first order linearisation of the model around the random effect expectation [START_REF] Mentré | Optimal design in random effect regression models[END_REF][START_REF] Retout | Fisher information matrix for non-linear mixed-effects models: evaluation and application for optimal design of enoxaparin population pharmacokinetics[END_REF][START_REF] Bazzoli | Fisher information matrix for nonlinear mixed effects multiple response models: evaluation of the appropriateness of the first order linearization using a pharmacokinetic/pharmacodynamic model[END_REF]. This approach has been implemented in several software programs [START_REF] Bazzoli | Design evaluation and optimisation in multiple response nonlinear mixed effect models: PFIM 3.0[END_REF][START_REF] Leonov | Optimal design for population PK/PD models[END_REF][START_REF] Gueorguieva | A program for individual and population optimal design for univariate and multivariate response pharmacokineticpharmacodynamic models[END_REF][START_REF] Nyberg | PopED: an extended, parallelized, nonlinear mixed effects models optimal design tool[END_REF] such as PFIM (INSERM, University Paris Diderot), POPED (University of Uppsala), POPDES (University of Manchester), and POPT (University of Otago), frequently used to design new studies in academia as well as in pharmaceutical companies [START_REF] Mentré | Current use and developments needed for optimal design in pharmacometrics: A study performed among DDMoRe's European Federation of Pharmaceutical Industries and Associations members[END_REF].

However, it has been shown that the use of the linearisation (LIN) approach is only appropriate if the variances of the random effects are small, or the nonlinearity is mild (Jones and Wang, 1999;Jones et al., 1999). Consequently, as pointed out by [START_REF] Han | Design of population studies of HIV dynamics[END_REF], when an optimal design is found using this approximation but the estimation is carried out using a true NLMEM, the performance of the design needs further investigation. The nonlinearity of a model with respect to its parameters is defined from the behaviour of the first order derivatives of the model function. Its consequences on the structural identifiability of a model have been studied by [START_REF] Walter | Identifiabilities and nonlinearities[END_REF]. The notion of level of nonlinearity ("mild" or "strong" as mentioned in this paper) is derived from the term "close to linear" introduced when evaluating a nonlinear model's behaviour [START_REF] Fletcher | A rapidly convergent descent method for minimization[END_REF]. Measures of nonlinearity have been studied in several publications [START_REF] Bates | Relative curvature measures of nonlinearity[END_REF][START_REF] Cook | Curvatures for parameter subsets in nonlinear regression[END_REF][START_REF] Smyth | Nonlinear regression[END_REF] in order to evaluate whether the "close to linear" condition is satisfied and to indicate if the linear approximation is reasonable or questionable.

When first-order linearisation is to be avoided, alternative approaches are necessary. Various new approaches have been proposed, but these are not always better than the usual first-order linearisation. For instance, linearisation of the model around the individual values of the random effects [START_REF] Retout | Further developments of the Fisher information matrix in nonlinear mixed effects models with evaluation in population pharmacokinetics[END_REF] or around the expected mode of the marginal likelihood [START_REF] Nyberg | PopED: an extended, parallelized, nonlinear mixed effects models optimal design tool[END_REF] has been proposed but is quite time consuming because Monte Carlo simulations are needed.

Other approaches based on the Laplace approximation [START_REF] Vong | Handling below limit of quantification data in optimal trial design[END_REF] or Monte Carlo integration [START_REF] Mielke | Approximations of the Fisher information for the construction of efficient experimental designs in nonlinear mixed effects models[END_REF] give correct predictions for the precision of parameter estimation, but they are also very time consuming. Another possible alternative for computing the Fisher information matrix in designing studies is the use of Gaussian quadrature rules. This consists in approximating integrals of functions with respect to a given probability density by a weighted sum of function values at abscissas chosen within the integration domain. It has been shown that adaptive Gaussian quadrature (AGQ) performs better than standard Gaussian quadrature (GQ) in estimation [START_REF] Pinheiro | Approximations to the log-likelihood function in the nonlinear mixed-effects model[END_REF]; the difference between the two approaches is that the grid of nodes is centred at the expectation of the random effects in GQ while it is centred at the conditional modes of the random effects in AGQ. Neither approach has ever been proposed for designing studies with different types of models, except in an example of an HIV dynamic model written with ordinary differential equations [START_REF] Guedj | Practical identifiability of HIV dynamics models[END_REF].

In this context, we aim to propose alternatives to linearisation for evaluating the predicted Fisher information matrix, based on GQ and AGQ. In order to challenge and investigate the performance of both new approaches as well as linearisation compared with CTS, we use examples of doseresponse trials inspired by the article of [START_REF] Plan | Performance comparison of various maximum likelihood nonlinear mixed-effects estimation methods for dose-response models[END_REF] comparing different estimation methods. Dose-response studies are of critical importance in drug development and need to be planned carefully [START_REF] Bretz | Practical considerations for optimal designs in clinical dose finding studies[END_REF][START_REF] Pronzato | Penalized optimal designs for dose-finding[END_REF][START_REF] Mcgree | Adaptive Bayesian compound designs for dose finding studies[END_REF] but little has been done to study their design in the context of NLMEM. We consider the sigmoid E max model, with various degrees of nonlinearity (i.e. different sigmoidicity coefficients), and in addition a linear model where M F can be calculated exactly.

We introduce the necessary notation for the design and model and present the new GQ and AGQ approaches developed for computing M F in Section 2. The performance of both approaches is evaluated and compared with CTS and LIN for different scenarios in Section 3 and is discussed in Section 4.

Computing Fisher information matrix in NLMEM with Gaussian quadratures

Design

The elementary design ξ i of individual i (i = 1, . . . , N) is defined by the number n i of observations and the design variables (x i1 , . . . , x in i ). Consequently, the population design for N individuals can be defined as Ξ = {ξ 1 , . . . , ξ N }. Usually, population designs are composed of a limited number Q of groups of individuals with identical designs in each group. Each of these groups is composed of a global elementary design ξ q and is performed in a number N q of subjects. The population design can thus be written as

Ξ = {[ξ 1 , N 1 ]; . . . ; [ξ Q , N Q ]}.

Nonlinear mixed effect model

We denote by y i the n i -vector of observations for the individual i obtained with the design ξ i and by f the known function describing the nonlinear structural model. The NLMEM linking the response y i to the samples ξ i = (x i1 , . . . , x in i ) can be written as

y i = f (φ i , ξ i ) + ε i , (1)
where ε i is the vector of random errors which follows a standard normal distribution N (0, σ 2 I n i )

and I n i is an n i × n i identity matrix. The vector of individual parameters φ i can be expressed as function g of µ, the P-vector of fixed effects and of b i , the P-vector of random effects for individual i. g can be additive (for normal parameters), so that the p th component of φ i is written as

φ ip = g(µ p , b ip ) = µ p + b ip , (2)
or exponential (for lognormal parameters), so

φ ip = g(µ p , b ip ) = µ p exp(b ip ). (3)
It is assumed that b i ∼ N (0, Ω), with Ω defined here as a diagonal variance-covariance matrix of size P × P. Each element ω 2 p of Ω represents the variance of the p th component of b i . As usual, the following assumptions are made: ε i |b i are independent between subjects, and ε i and b i are independent for each subject. The model can also be written as

y i = f (g(µ, b i ), ξ i ) + ε i . (4)
Let λ = (ω 1 , ..., ω P , σ ) be the vector of standard deviations of random effects and standard deviation of random error. Ψ, the vector of all population parameters to be estimated, is Ψ = (µ , λ ) .

Fisher information matrix

The individual Fisher information matrix is defined from the likelihood of the observations y i of individual i, which is expressed by the following integral

L(y i ; Ψ) = R P p(y i |b i ; Ψ) p(b i ) db i , (5)
where p(y i |b i ; Ψ) is the conditional density of the observations y i given the random effects b i , and

p(b i ) is the density of b i ∼ N (0, Ω).
Then, the Fisher information matrix for the elementary design ξ i is the following expectation taken with respect to the distribution of the observations p(y; Ψ)

M F (Ψ; ξ i ) = E ∂ log L(y i ; Ψ) ∂ Ψ ∂ log L(y i ; Ψ) ∂ Ψ . ( 6
)
As the individuals are independent, the Fisher information matrix M F (Ψ, Ξ) for a population design Ξ is defined by the sum of the N elementary matrices M F (Ψ, ξ i ), so that

M F (Ψ, Ξ) = N ∑ i=1 M F (Ψ, ξ i ). (7)
In the case of a limited number Q of elementary designs ξ q , we have

M F (Ψ, Ξ) = Q ∑ q=1 N q M F (Ψ, ξ q ). (8)
However, there is no analytical expression for the likelihood L(y i ; Ψ) in NLMEM. Approximations such as linearisation can be used to approximate the Fisher information matrix [START_REF] Mentré | Optimal design in random effect regression models[END_REF][START_REF] Retout | Fisher information matrix for non-linear mixed-effects models: evaluation and application for optimal design of enoxaparin population pharmacokinetics[END_REF][START_REF] Bazzoli | Fisher information matrix for nonlinear mixed effects multiple response models: evaluation of the appropriateness of the first order linearization using a pharmacokinetic/pharmacodynamic model[END_REF]. The predicted M F after linearisation is a block matrix with a block corresponding to derivatives of the log-likelihood with respect to the fixed effects, a block for derivatives with respect to the standard derivation terms and a block containing mixed derivatives with respect to all parameters. In our work, the block of mixed derivatives was set to 0 for linearisation, based on publications showing the better performance of the block diagonal expression compared with the full one [START_REF] Mielke | Some considerations on the Fisher information in nonlinear mixed effects models[END_REF][START_REF] Fedorov | Optimal Design for Nonlinear Response Models[END_REF][START_REF] Nyberg | Methods and software tools for design evaluation for population pharmacokineticspharmacodynamics studies[END_REF].

Gaussian quadrature and adaptive Gaussian quadrature

When linearisation is to be avoided, alternatives such as Gaussian quadratures can be used to express the likelihood analytically. For simplicity, we omit the index i for the individual in this section. The elementary Fisher information matrix M F (Ψ, ξ ) in ( 6) for an individual with elementary design ξ of n observations can also be written as

M F (Ψ, ξ ) = R n h(y; Ψ) p(y; Ψ) dy, (9)
where p(y; Ψ) is the marginal density of y and

h(y; Ψ) = ∂ log L(y; Ψ) ∂ Ψ ∂ log L(y; Ψ) ∂ Ψ = ∂ L(y;Ψ) ∂ Ψ ∂ L(y;Ψ) ∂ Ψ L(y; Ψ) 2 . ( 10
)
First, we propose to evaluate the integral of h(y; Ψ) with respect to p(y; Ψ) in ( 9) by stochastic integration. We generate M vectors of observations y m (m = 1, . . . , M), each vector with the same design ξ of n samples, under the same probability distribution p(y; Ψ). Thus, an approximation of ( 9) can be obtained by

M F (Ψ, ξ ) = 1 M M ∑ m=1 h(y m ; Ψ). (11)
Then, in order to obtain an analytical expression of h(y m ; Ψ) as expressed in (10) for each simulated y m , we need to compute the likelihood L(y m ; Ψ). We define η = Ω -1/2 b, then η ∼ N (0, I P ), where I P is the identity matrix of dimension P × P. Consequently, the likelihood in (5) can also be written as

L(y m ; Ψ) = R P p(y m |η; Ψ)p(η)dη, (12)
where p(η) is the density of N (0, I P ). This integral can be evaluated numerically by the Gaussian quadrature approach, as a weighted sum of p(y m |η; Ψ) evaluated at pre-determined nodes chosen within the distribution of η [START_REF] Pinheiro | Approximations to the log-likelihood function in the nonlinear mixed-effects model[END_REF]. Quadrature weights and nodes for approximating one-dimensional integrals are those of the standard Gauss-Hermite rule [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF][START_REF] Golub | Calculation of Gauss quadrature rules[END_REF][START_REF] Golub | Some modified matrix eigenvalue problems[END_REF]. The nodes are the roots of the Hermite polynomial of degree κ, obtained from successive derivatives of exp(η 2 ) [START_REF] Press | Numerical Recipes in C, The Art of Scientific Computing[END_REF]. As in [START_REF] Pinheiro | Approximations to the log-likelihood function in the nonlinear mixed-effects model[END_REF], the P-dimension integral here was transformed into P successive one-dimensional integrals. The distribution of nodes and weights used in standard GQ does not take into account the nature of the integrand, which is why standard GQ requires a high number of nodes in order to provide a correct estimation of parameters. AGQ, by centring and scaling the standard quadrature nodes, places the nodes according to the areas of high density and improves the approximation. Consequently it requires fewer nodes compared with standard GQ, while giving better performance in estimation [START_REF] Pinheiro | Approximations to the log-likelihood function in the nonlinear mixed-effects model[END_REF].

In order to approximate each one of these integrals in L(y m ; Ψ) using standard GQ, for all simulated y m , κ quadrature nodes and weights are selected from the standard normal distribution N (0, 1). Then the integral of dimension P can be approximated by the weighted sum over K nodes η k , where K = κ P and η k are vectors in R P , associated with weights w k (k = 1, . . . , K). So, we obtain the following analytical expression

L(y m ; Ψ) = K ∑ k=1 w k p(y m |η k ; Ψ), (13)
where

y m |η k ∼ N ( f k , σ 2 I n ) with f k = f (g(µ, Ω 1/2 η k ), ξ
). Thus, we can write

p(y m |η k ; Ψ) = (2πσ 2 ) -n/2 exp -1 2σ 2 (y m -f k ) (y m -f k ) , and (14)
log p(y m |η k ; Ψ) = - n 2 log 2π - n 2 log σ 2 - 1 2σ 2 (y m -f k ) (y m -f k ). ( 15
)
The likelihood L(y m ; Ψ) can also be approximated by AGQ. Since a vector of standardised random effects η m is simulated when generating each y m , we select K = κ P nodes η mk and weights w k based on p(η) which is the N (η m , I P ) density. We can write

L(y m ; Ψ) = R P p(y m |η; Ψ) p(η) p(η) p(η)dη (16) = K ∑ k=1 w k p(η mk ) p(η mk ) p(y m |η mk ; Ψ), (17)
where y m |η mk ∼ N ( f mk , σ 2 I n ) with f mk = f (g(µ, Ω 1/2 η mk ), ξ ), and the weights for AGQ can be written as

w k p(η mk ) p(η mk ) = w k (2π) -P/2 exp(-1 2 η mk η mk ) (2π) -P/2 exp(-1 2 (η mk -η m ) (η mk -η m )) (18) = w k exp - 1 2 η mk η mk -(η mk -η m ) (η mk -η m ) . (19)
Finally, using ( 13) or ( 17) in ( 10) and in ( 11), the individual Fisher information matrix can be approximated as

M F (Ψ, ξ ) = 1 M M ∑ m=1 ∑ K k=1 w mk ∂ p mk (y m ;Ψ) ∂ Ψ ∑ K k=1 w mk ∂ p mk (y m ;Ψ) ∂ Ψ ∑ K k=1 w mk p mk (y m ; Ψ) 2 (20) = 1 M M ∑ m=1 ∑ K k=1 w mk p mk (y m ; Ψ) ∂ log p mk (y m ;Ψ) ∂ Ψ ∑ K k=1 w mk p mk (y m ; Ψ) ∂ log p mk (y m ;Ψ) ∂ Ψ ∑ K k=1 w mk p mk (y m ; Ψ) 2 . ( 21
)
When using GQ, w mk = w k and p mk (y m ; Ψ) = p(y m |η k ; Ψ) as defined in (14). When using AGQ,

w mk = w k p(η mk )
p(η mk ) as defined in ( 19) and p mk (y m ; Ψ) = p(y m |η mk ; Ψ) as defined in ( 14) with f mk instead of f k .

As Ψ = (µ , λ ) , the Fisher information matrix M F (Ψ, ξ ) can be written as a block matrix

M F (Ψ; ξ ) =   M F (µ; ξ ) M F (µ, λ ; ξ ) M F (µ, λ ; ξ ) M F (λ ; ξ )   . ( 22
)
To compute M F (µ; ξ ), we derive the log-likelihood with respect to each element µ p of the vector of fixed effects µ = (µ 1 , ..., µ P ) , using

∂ log p mk (y m ; Ψ) ∂ µ p = 1 σ 2 ∂ f k ∂ µ p (y m -f k ), (23)
To compute M F (λ ; ξ ), we derive the log-likelihood with respect to each element ω p or σ of the vector of standard deviations λ = (ω 1 , ..., ω P , σ ) , using

∂ log p mk (y m ; Ψ) ∂ ω p = 1 σ 2 ∂ f k ∂ ω p (y m -f k ), and (24)
∂ log p mk (y m ; Ψ) ∂ σ = - n σ + 1 σ 3 (y m -f k ) (y m -f k ), (25)
with f mk instead of f k when using AGQ instead of GQ. Finally, using (23, 24, 25) in ( 21), we obtain each element of M F (µ; ξ ) in ( 22).

Evaluation by simulation

Evaluation models and scenarios

Our evaluation examples were inspired by a previous dose-response simulation study [START_REF] Plan | Performance comparison of various maximum likelihood nonlinear mixed-effects estimation methods for dose-response models[END_REF] which mimicked clinical trials including 100 subjects with identical individual design ξ and investigating several dose levels among 0, 100, 300 and 1000 dose units. First, we considered a linear model of the dose-response relationship describing, for subject i,

an effect E i j for the j th dose d i j with two structural parameters: baseline E 0 and slope S.

E i j = E 0 i + S i × d i j + ε i j . ( 26
)
The vector of the fixed effects µ was composed of µ E 0 = 5 and µ S = 0.06. The random effect model was additive with standard deviations ω E 0 = 1.5 and ω S = 0.018. The standard deviation σ of random error ε i j was equal to 2. We studied two designs, a rich one ξ = (0, 100, 300, 1000) and a sparse one ξ = (100, 300). We chose to evaluate first a simple linear mixed effect model (denoted by M1) because in this case the analytical form of the Fisher information matrix is known and so we could easily check whether our calculation and implementation are correct.

Next, we considered a sigmoid E max model of the dose-response relationship, describing for subject i, an effect E i j for the j th dose d i j with the following structural parameters: baseline E 0 , maximal effect E max , dose ED 50 at which 50% of the maximal effect is achieved and sigmoidicity coefficient γ, i.e. the degree of nonlinearity of the function shape.

E i j = E 0 i + E max i × d γ i j ED γ 50 i + d γ i j + ε i j . ( 27
)
We considered a model with γ = 1 and a model with γ = 3. In each model, γ was known with the same value in all subjects (i.e. µ γ not estimated, ω γ = 0). We first studied models with two structural parameters: E max and ED 50 where E 0 was fixed at the same value 5 in all subjects (i.e. µ E 0 not estimated, ω E 0 = 0). The fixed effects µ were µ E max = 30, µ ED 50 = 500 and the random effects were exponential with ω E max = ω ED 50 = 0.3. We denote these models M2 and M3 for γ = 1 and γ = 3, respectively. Then, we studied models with three structural parameters: E max , ED 50 and E 0 , where E 0 can vary from one subject to another and was estimated. The fixed effects µ were µ E max = 30, µ ED 50 = 500, µ E 0 = 5 and the random effects were exponential with ω

E max = ω ED 50 = ω E 0 = 0.3.
We denote these models M4 and M5 for γ = 1 and γ = 3 respectively. The standard deviation σ of random error ε i j was equal to 2 in all models. Designs ξ with 2, 3, or 4 doses among (0, 100, 300, 1000) were studied. The list of the 5 models above and the designs studied are detailed in Table 1.

The rich design is referred to as R and the sparse design as S. Examples of one dataset simulated We plotted the first and second order derivatives to examine graphically the model nonlinearity (Figure 2). The first derivatives with respect to E max and ED 50 vary over design samples and are presented in Figure 2A. We notice that the derivative with respect to E max depends on ED 50 and the one with respect to ED 50 depends on both E max and ED 50 . Thus the optimal design depends on values of E max and ED 50 but not on E 0 . One class of nonlinearity measures is based on the second order derivatives [START_REF] Smyth | Nonlinear regression[END_REF]: if the second derivatives of model f 1 are smaller in absolute values than f 2 , then model f 1 has "closer to linear" behaviour. The second derivatives that are different from 0 are plotted against doses in Figure 2B. The magnitudes of the curves are greater when γ = 3 than when γ = 1, indicating that model nonlinearity increases with γ. This observation was confirmed by [START_REF] Plan | Performance comparison of various maximum likelihood nonlinear mixed-effects estimation methods for dose-response models[END_REF] who showed decreasing performance of the linearisation approach along the γ-increase when estimating parameters for E max models.

Comparison of standard errors between predictions and clinical trial simulation

Our aim was to evaluate and compare the standard errors (SE) predicted using M F by GQ and AGQ with those obtained by CTS and LIN.

We performed clinical trial simulations in R 2.14.1 with the model, parameters and designs above. For each scenario, 1000 datasets of 100 subjects were simulated. Each simulated dataset was then analysed in MONOLIX 3.2 (www.lixoft.eu). Population parameters were estimated by the SAEM algorithm [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF]. Next, the SE by CTS, defined as the unbiased sample estimate of the standard deviation from the 1000 parameter estimates, were calculated. The relative standard errors (RSE) by CTS were defined as the ratio of the SE to the true value of parameters.

As in the publication by [START_REF] Retout | Further developments of the Fisher information matrix in nonlinear mixed effects models with evaluation in population pharmacokinetics[END_REF], the 95% confidence intervals for each RSE by CTS are given, computed as q 1 999 RSE 2 ; q 2 999 RSE 2 where q 1 and q 2 are, respectively, the 2.5% and 97.5% quantiles of the χ 2 distribution with 999 degrees of freedom. We also computed the determinant of the variance-covariance matrix of all parameter estimates, normalised according to the standard definition of D-criterion, i.e. det(M F ) -1/(2P+1) .

In parallel, we first used PFIM 3.2.2 (www.pfim.biostat.fr) with the model, parameters and designs above to predict M F by LIN. Second, we implemented the calculations for GQ and AGQ (Section 2.4) in R 2.14.1 in a working version of PFIM, using the function gauss.quad.prob of the R package statmod and the function as.weight of the R package plink for multidimensional quadratures, which is an extension of gauss.quad.prob. We used M= 1000 simulations to compute the expectation of h (see equation 11). We found that κ = 100 nodes are needed for GQ and only 30 nodes for AGQ to guarantee stable results with all studied models. The predicted SE of parameters were calculated as the square root of the diagonal terms of M -1 F and the predicted RSE were defined as the ratio of the SE to the true value of parameters. We also calculated the normalised determinant of the predicted variance-covariance matrix of parameter estimates (as defined above).

Finally, we compared the RSE and the normalised determinant of the variance-covariance matrix predicted using LIN, GQ and AGQ with those obtained by CTS.

Results

For the linear dose-response model (Figure 3), the RSE predicted by LIN, GQ and AGQ were close and in the same range as the ones calculated from CTS. In this case, the linearisation approach provided the exact calculation of M F and of the true SE. As expected, the RSE were higher in the sparse design with 2 doses (Figure 3, M1S) than in the rich design with 4 doses (Figure 3, M1R), not so much for the fixed effects but more for the standard deviations of the random effects and of the random error.

For the E max model with 2 structural parameters and 2 random effects, when γ = 1 (Figure 4, M2R-M2S), the three prediction approaches gave RSE close to those obtained by CTS. However, when γ = 3 (Figure 4, M3R-M3S), with the sparse design, the linearisation approach over-predicted the RSE for ω ED 50 and especially for σ (42% predicted by LIN versus 15% by GQ, 14% by AGQ

M3S

Figure 4: Relative standard error RSE (%) for the E max dose-response model with fixed baseline, obtained by clinical trial simulation (white bar) versus those predicted by linearisation (light grey bar), Gaussian quadrature (dark grey bar) and adaptive Gaussian quadrature (black bar), for the fixed effects (µ), the standard deviation of the random effects (ω), and of the random error (σ ), with rich or sparse designs: M2R (γ = 1, 4 doses per subject), M2S (γ = 1, 2 doses per subject), M3R (γ = 3, 4 doses per subject), M3S (γ = 4, 2 doses per subject). The 95% confidence intervals of the RSE obtained by clinical trial simulation are displayed on top of the white bars.

and 17% obtained by CTS). With the rich design, the level of nonlinearity seemed to have less impact on the performance of linearisation: like GQ and AGQ, LIN adequately predicted RSE of all model parameters, in the same range as the ones obtained by CTS.

RSE (%)

0 10 20 30 40

µ E max µ ED 50 µ E 0 ω E max ω ED 50 ω E 0 σ M4R CTS LIN GQ AGQ RSE (%) 0 10 20 30 40 µ E max µ ED 50 µ E 0 ω E max ω ED 50 ω E 0 σ M4S RSE (%) 0 10 20 30 40 µ E max µ ED 50 µ E 0 ω E max ω ED 50 ω E 0 σ M5R RSE (%) 0 10 20 30 40 µ E max µ ED 50 µ E 0 ω E max ω ED 50 ω E 0 σ

M5S

Figure 5: Relative standard error RSE (%) for the E max dose-response model with estimated baseline, obtained by clinical trial simulation (white bar) versus those predicted by linearisation (light grey bar), Gaussian quadrature (dark grey bar) and adaptive Gaussian quadrature (black bar), for the fixed effects (µ), the standard deviation of the random effects (ω), and of the random error (σ ), with rich or sparse designs: M4R (γ = 1, 4 doses per subject), M4S (γ = 1, 3 doses per subject), M5R (γ = 3, 4 doses per subject), M5S (γ = 4, 3 doses per subject). The 95% confidence intervals of the RSE obtained by clinical trial simulation are displayed on top of the white bars.

For the E max model with 3 structural parameters and 3 random effects, when γ = 1 (Figure 5, M4R-M4S), the three prediction approaches gave RSE in the same range as those obtained by CTS. However, as in the previous model, when γ = 3 (Figure 5, M5R-M5S), the performance of the linearisation approach deteriorated with the sparse design and it over-predicted the RSE of ω E 0 (27% predicted by LIN versus 19% by GQ, 21% by AGQ and 20% obtained by CTS). With the rich design, the four approaches gave RSE that were close.

It is of note, as reported in Figures 34-5, that the RSE obtained by CTS were representative of the true RSE because rather small 95% confidence intervals were found. Figure 6 represents the normalised determinants of the variance-covariance matrix obtained by the four approaches for

M1R M1S M2R M2S M3R M3S M4R M4S M5R M5S Det(M F ) -1 (2P+1)
0.000 0.001 0.002 0.003 0.004 0.005

CTS LIN GQ AGQ
Figure 6: Determinant of the variance-covariance matrix of parameter estimates, normalised by inverse of the total number of parameters (2P + 1), obtained by clinical trial simulation (white bar) versus those predicted by linearisation (light grey bar), Gaussian quadrature (dark grey bar) and adaptive Gaussian quadrature (black bar), for different models as defined in Table 1.

different scenarios. Overall, the predicted variances were lower than the ones observed from CTS.

However, for the scenario M3S, LIN predicted a determinant of the variance-covariance matrix that is larger than the one obtained by CTS.

Regarding computing time, the linearisation approach was much faster than the two other prediction approaches. For the models with 2 random effects, the runtimes were less than 1 s by LIN in PFIM 3.2.2, about 6 h by GQ, 20 minutes by AGQ, and 10 h by CTS. For the models with 3 random effects, the runtimes were less than 1 s by LIN, about 110 h by GQ, 11 h by AGQ and 12 h by CTS.

Discussion

The present work proposes new approaches to the evaluation of the Fisher information matrix in NLMEM, using Gaussian quadrature and adaptive Gaussian quadrature when designing longitudinal studies. This is an alternative to first-order linearisation, the most commonly used approach in the field of optimal design in NLMEM. We investigated the performance of these approaches in predicting parameter standard errors and determinants of the Fisher matrix for dose-response models, as compared with linearisation and clinical trial simulation. These predictions are important for design evaluation which is always the first step when planning an upcoming study and correct values of the D-criterion are needed for design optimisation.

When the nonlinearity was mild (γ = 1), the linearisation approach gave correct predictions of the RSE, close to the empirical ones obtained by simulation. When nonlinearity increased (γ = 3), linearisation correctly predicted the RSE of fixed effects but over-predicted, in the case of the sparse design, the RSE of ω ED 50 and σ . This poor performance could be partly explained by the behaviour of the first derivatives of the model function, which are used in the calculation of the M F block corresponding to the standard deviation terms [START_REF] Bazzoli | Fisher information matrix for nonlinear mixed effects multiple response models: evaluation of the appropriateness of the first order linearization using a pharmacokinetic/pharmacodynamic model[END_REF]. As shown in Figure 2A, the derivatives of f when γ = 3 (especially with respect to ED 50 ) are much more sensitive to design samples as compared with γ = 1. This might be why when γ = 3, with the sparse design of 2 doses (300, 1000) only, RSE of some standard deviation terms were not well predicted by linearisation. Further evaluations of other examples are needed to understand better the difference in performance of this approach between the sparse and rich designs.

GQ and AGQ gave correct predictions of RSE overall, close to the RSE obtained by CTS, in spite of a slight discrepancy versus CTS for some variability terms. When experimenters design population studies, they are more interested in the magnitude of the RSE than the exact value, so these approaches are relevant to evaluate M F when linearisation is to be avoided. Here, we used κ = 100 nodes for GQ and κ = 30 nodes for AGQ to approximate the integrals correctly, when compared with CTS, in all studied models (so the total number of model evaluations was K = κ P , P is the number of random effects). Even fewer nodes were sufficient to obtain stable results by AGQ for models with milder nonlinearity. [START_REF] Clarkson | Using spherical-radial quadrature to fit generalized linear mixed effects models. journal of computational and graphical statistics[END_REF] needed κ = 15 nodes for spherical-radial quadrature when estimating parameters of a logistic mixed model with two random effects. [START_REF] Gotwalt | Fast computation of designs robust to parameter uncertainty for nonlinear settings[END_REF] also used spherical-radical quadrature to compute the Bayesian design criteria (in a nonlinear model without random effects) and showed that the number of function evaluations to be performed increases as the square of the number of parameters P. In order to evaluate each elementary Fisher information matrix, we also needed to perform an integration over the marginal distribution p(y, Ψ) of the observations y, which depends on the design Ψ (equation 9). Theoretically, this could have been evaluated by Gaussian quadrature rules as well, based on the distribution of y. However, this integral has the same dimension as the number of samples in the individual design ξ , and in practice Gaussian quadrature rules can only be applied when the number of design samples is small. We proposed an alternative method based on stochastic integration, simulating M = 1000 vectors of observations under the same probability distribution as y, which was less time consuming than using a second quadrature, especially for rich designs. [START_REF] Guedj | Practical identifiability of HIV dynamics models[END_REF] also found in their study, with a similar approach using a simulation step, that M = 1000 was sufficient to evaluate M F in a model of HIV dynamics. Further evaluations are needed to select the best balance between M and κ.

As a consequence, AGQ is much more computationally efficient than GQ for designing trials.

For instance, in our study, for the models with 2 random effects, the runtimes were about 6 h by GQ, 20 min by AGQ and 10 h by CTS; and for the models with 3 random effects, about 110 h by GQ, 11 h by AGQ and 12 h by CTS. Here, we considered variability in all parameters; the number of random effects is the number of structural parameters, which is not always the case. Indeed, there are usually only one or two random effects in discrete data mixed models [START_REF] Savic | Implementation and evaluation of the SAEM algorithm for longitudinal ordered categorical data with an illustration in pharmacokineticspharmacodynamics[END_REF][START_REF] Abebe | Bayesian D-optimal designs for the two parameter logistic mixed effects model[END_REF] or survival frailty models [START_REF] Vigan | Evaluation of estimation methods for repeated time to event models: application to analysis of bone events during treatment of Gaucher Disease[END_REF], so the computing times by these new approaches are still reasonable compared with CTS. Another advantage of GQ and AGQ compared with CTS is that the computing time in GQ or AGQ does not depend on the number of subjects in the population design because of the properties of the population M F (equations 7 and 8). Consequently, a gain in computing time would be more obvious when designing studies in large populations. For instance, for a model with two random effects, simulation and fitting of 1000 datasets of 1000 subjects or 10 000 subjects will take about 30 hours and 800 hours respectively, while GQ will always take 6 hours and AGQ 20 minutes to evaluate the same design in all subjects.

GQ and AGQ would more easily allow evaluation of a larger number of designs because of the property (8), while CTS is certainly more time consuming because one has to simulate and analyse again new datasets whenever adding or removing an elementary design. One perspective of this work is to evaluate the performance of alternative sampling techniques in the stochastic integration part of our AGQ approach, for example we would try using Latin Hypercube sampling, which might be less time consuming than standard random sampling, as suggested [START_REF] Ueckert | Application of Quasi-Newton algorithms in optimal design[END_REF].

Here we considered a diagonal variance-covariance matrix Ω, but one may want to take into account the correlation between random effects in the model, for example, between those on E max and ED 50 as in [START_REF] Plan | Performance comparison of various maximum likelihood nonlinear mixed-effects estimation methods for dose-response models[END_REF]. Therefore, it would be interesting to include the full Ω in GQ and AGQ approaches as well, based on calculations proposed by [START_REF] Dumont | Influence of covariance between random effects in design for nonlinear mixed effect models with an illustration in paediatric pharmacokinetics[END_REF]. Moreover, we should also evaluate the GQ and AGQ with different levels of between-subject variability or with increasing variance as in growth data. It would also be necessary to evaluate these approaches not only for dose-response models but also for discrete data mixed models and survival frailty models where the linearisation approach to design would work poorly.

Other possible alternatives to linearisation for evaluation of M F when designing studies could be spherical radial integration [START_REF] Monahan | Spherical radial integration rules for Bayesian computations[END_REF] or stochastic approach and importance sampling. Stochastic approach to evaluate M F can be inspired by the SAEM algorithm [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF], using Louis's [START_REF] Louis | Finding the observed information matrix when using the EM algorithm[END_REF] or Oakes's [START_REF] Oakes | Direct calculation of the information matrix via the EM algorithm[END_REF] formulas for M F , with simulation of individual parameters and stochastic approximation of expectations. Importance sampling is more time consuming than AGQ in estimation, but has the advantage of being versatile in handling distributions other than the normal for both random effects and residual errors [START_REF] Pinheiro | Approximations to the log-likelihood function in the nonlinear mixed-effects model[END_REF]. However, it has not yet been assessed for design evaluation.

Once an appropriate approach to evaluating designs and to computing the D-criterion is chosen, the determinant of the Fisher information matrix can be maximised using the iterative Fedorov-Wynn algorithm [START_REF] Fedorov | Theory of Optimal Experiments[END_REF][START_REF] Wynn | Results in the theory and construction of D-optimum designs[END_REF] within a finite set of possible designs, based on the Kiefer-Wolfowitz equivalence theorem [START_REF] Kiefer | The equivalence of two extremum problems[END_REF], as implemented in design software such as PFIM [START_REF] Retout | Design in nonlinear mixed effects models: Optimisation using the Fedorov-Wynn algorithm and power of the wald test for binary covariates[END_REF], PkStaMP and PopDes. As suggested by [START_REF] Fedorov | Optimal Design for Nonlinear Response Models[END_REF], we would compute the normalised predicted variance for all possible combinations of candidate points, which hits its maximum value of dim(Ψ) at the support points of the D-optimal design. However, a limitation of the current approaches proposed here is that they require a priori knowledge of the model and its parameters, which leads to designs that are only locally optimal.

Sensitivity analyses with respect to the model and the parameter values would be necessary to quantify how the results vary. Alternatives, such as adaptive designs [START_REF] Leonov | An adaptive optimal design for the Emax model and its application in clinical trials[END_REF][START_REF] Foo | Adaptive optimal design for bridging studies with an application to population pharmacokinetic studies[END_REF] or robust designs based on Bayesian criteria [START_REF] Han | Bayesian experimental design for nonlinear mixed-effects models with application to HIV dynamics[END_REF][START_REF] Dette | Efficient experimental designs for sigmoidal growth models[END_REF][START_REF] Gotwalt | Fast computation of designs robust to parameter uncertainty for nonlinear settings[END_REF][START_REF] Abebe | Bayesian D-optimal designs for the two parameter logistic mixed effects model[END_REF] would be interesting to explore.

To conclude, when computing the Fisher information matrix in NLMEM for design evaluation and optimisation, the linearisation approach is accurate for models with mild nonlinearity, as has already been demonstrated for several PK [START_REF] Retout | Further developments of the Fisher information matrix in nonlinear mixed effects models with evaluation in population pharmacokinetics[END_REF][START_REF] Nguyen | Design evaluation and optimisation in crossover pharmacokinetic studies analysed by nonlinear mixed effects models[END_REF], PK/PD [START_REF] Bazzoli | Fisher information matrix for nonlinear mixed effects multiple response models: evaluation of the appropriateness of the first order linearization using a pharmacokinetic/pharmacodynamic model[END_REF] and viral dynamic models [START_REF] Retout | Design in nonlinear mixed effects models: Optimisation using the Fedorov-Wynn algorithm and power of the wald test for binary covariates[END_REF][START_REF] Guedj | Design evaluation and optimization for models of hepatitis C viral dynamics[END_REF].

This procedure is very fast and relatively simple, is available in several design software programs, and is a very useful tool for designing longitudinal studies while avoiding extensive simulations.

However, when the models are complex and have never been evaluated in design approaches, nonlinearity measures [START_REF] Bates | Relative curvature measures of nonlinearity[END_REF][START_REF] Cook | Curvatures for parameter subsets in nonlinear regression[END_REF]) of the studied model should be investigated for the evaluated or the optimised design before drawing final conclusions. If the nonlinearity is strong and therefore the linearisation approximation is to be avoided, then GQ and AGQ are relevant alternatives: GQ is very time consuming and is only useful for a small number of random effects, while AGQ is much less time consuming while providing adequate results. An appropriate choice for the number of quadrature nodes is to be defined, depending on the level of nonlinearity. Further evaluations for other types of models are needed before implementing AGQ in PFIM software for efficiently designing trials while reducing the number of samples per subject, which can be very important both ethically and practically when performing studies in patients.

Figure 1 :

 1 Figure1: Individual response versus dose profiles (in grey) for one dataset simulated using linear model (M1R), E max model with sigmoidicity coefficient γ = 1 (M4R) and E max model with γ = 3 (M5R), with 4 doses per subject (0, 100, 300, 1000). The black curves represent the response predicted by the model for the fixed effect parameters.

3 Figure 2 :

 32 Figure 2: First order (panel A) and second order (panel B) derivatives of the function f of the E max model with respect to the structural parameters E max and ED 50 for two sigmoidicity coefficients γ = 1 and γ = 3.

Table 1 :

 1 Different models and designs evaluated

	Structural	Number of	Number of doses	Doses	Name
	model	parameters P	per subject n
			4	(0, 100, 300, 1000) M1R
	Linear	2
			2	(100, 300)	M1S
	E max	γ = 1	4	(0, 100, 300, 1000) M2R
	(E 0 not estimated,	2	2	(300, 1000)	M2S
	fixed at same value		4	(0, 100, 300, 1000) M3R
		γ = 3
	in all subjects)		2	(300, 1000)	M3S
	E max	γ = 1	4	(0, 100, 300, 1000) M4R
	(E 0 estimated,	3	3	(0, 300, 1000)	M4S
	varying		4	(0, 100, 300, 1000) M5R
		γ = 3
	between subjects)		3	(0, 300, 1000)	M5S

  Relative standard error RSE (%) for the linear dose-response model obtained by clinical trial simulation (white bar) versus those predicted by linearisation (light grey bar), Gaussian quadrature (dark grey bar) and adaptive Gaussian quadrature (black bar), for the fixed effects (µ), the standard deviation of the random effects (ω), and of the random error (σ ), with rich or sparse designs: M1R (4 doses per subject), M1S (2 doses per subject). The 95% confidence intervals of the RSE obtained by clinical trial simulation are displayed on top of the white bars.

	M1R							M1S
		40
	RSE (%)	10 20 30	CTS LIN GQ AGQ				RSE (%)	10 20 30 40
		0						0
				µ E 0	µ S	ω E 0	ω S	σ	µ E 0	µ S	ω E 0	ω S	σ
	Figure 3: RSE (%) M2R	10 20 30 40	CTS LIN GQ AGQ				RSE (%) M2S	10 20 30 40
			0						0
				µ E max	µ ED 50	ω E max	ω ED 50	σ	µ E max	µ ED 50	ω E max	ω ED 50	σ
	M3R
	RSE (%)	10 20 30 40					RSE (%)	10 20 30 40
			0						0
				µ E max	µ ED 50	ω E max	ω ED 50	σ	µ E max	µ ED 50	ω E max	ω ED 50	σ

Acknowledgements

The authors thank IFR02 -Paris Diderot University and Hervé Le Nagard for the use of the "Centre de Biomodélisation".