J Antic
email: j.antic@envt.fr

C M Laffont

D Chafaï

D Concordet

Comparison of nonparametric methods in nonlinear mixed effects models

Keywords: Nonlinear mixed effects models, Nonparametric, NPML, NPEM, NONMEM, SNP, EBEs

During the drug development, nonlinear mixed effects models are routinely used to study the drug's pharmacokinetics and pharmacodynamics. The distribution of random effects is of special interest because it allows to describe the heterogeneity of the drug's kinetics or dynamics in the population of individuals studied. Parametric models are widely used, but they rely on a normality assumption which may be too restrictive. In practice, this assumption is often checked using the empirical distribution of random effects' empirical Bayes estimates. Unfortunately, when data are sparse (like in patients phase III clinical trials), this method is unreliable. In this context, nonparametric estimators of the random effects distribution are attractive. Several nonparametric methods (estimators and their associated computation algorithms) have been proposed but their use is limited. Indeed, their practical and theoretical properties are unclear and they have a reputation for being computationally expensive. The article evaluates four nonparametric methods in comparison with the usual parametric method. Statistical and computational features are reviewed and practical performances are compared in simulation studies mimicking real pharmacokinetic analyses. The nonparametric methods seemed very useful when data are sparse. On a simple pharmacokinetic model, all the nonparametric methods performed roughly equivalently. On a more challenging pharmacokinetic model, differences between the methods were clearer.

Introduction

Nonlinear mixed effects (NLME) models are increasingly used in several biomedical applications. Especially, they have become a routine tool for the analysis of pharmacokinetic (PK) and/or pharmacodynamic (PD) data obtained in clinical trials during drug development. Furthermore, the US Food and Drug Administration strongly recommends their use for regulatory review (FDA, 1999). For each individual, PK data include the dosing history (amount of drug administered, route and time of administration), measured drug concentrations in blood or in plasma, measurement times and values of clinically relevant covariates (e.g. sex, age, weight, renal function...). The drug's systemic concentration-time profile depends on fixed effects and on individual-specific random effects quantifying the drug's absorption, distribution and elimination. NLME models allow to estimate the distribution P * of random effects, even when within-individual information is sparse. They are particularly interesting in phase II or phase III clinical trials that involve a large group of individuals representative of the target population, with often few information per individual. The objective of these analyses is to give a rationale for dosing recommendations and identify subpopulations of patients that require a dosing adjustment.

The most widely used model is parametric: it is usually assumed that P * is a (log-) normal distribution. The model is said parametric because the distribution of random effects P * depends only on a finite number of parameters: its mean and variance. In this framework, the maximum likelihood estimator (MLE) is generally favored because it is consistent and efficient. The nonlinear relationship between the observations and the random effects makes the likelihood not explicit. Several algorithms have been proposed to compute approximations of the MLE. The most popular are First Order Conditional Estimates with Interaction (FOCEI), First Order Conditional Estimates (FOCE) and First Order (FO) algorithms. They are implemented in the gold standard software NONMEM [START_REF] Boeckmann | NONMEM Users Guide[END_REF].

However, the parametric assumption may be too restrictive to describe very heterogeneous populations. In practice, this can be checked by looking at the empirical distribution of individual predictions of random effects, known as empirical Bayes estimates (EBEs) [START_REF] Pinheiro | Mixed-Effects Models in S and S-PLUS[END_REF], EMEA (2007)). The problem is that, when the number of measurements per individual is low (the dataset is qualified as sparse), the EBEs are strongly influenced by the parametric assumption and their empirical distribution is unreliable [START_REF] Savic | Un)informativeness of empirical Bayes estimate-based diagnostics[END_REF].

In that context, nonparametric (NP) estimators, which allow P * to live in an infinite dimensional space, are attractive. Several NP methods have been proposed (by method, we mean an estimator plus its computation algorithm). [START_REF] Laird | Nonparametric Maximum Likelihood Estimation of a Mixing Distribution[END_REF], [START_REF] Pfanzagl | Consistency of maximum likelihood estimators for certain nonparametric families, in particular: mixtures[END_REF] and [START_REF] Lindsay | The geometry of mixture likelihoods : a general theory[END_REF] studied the MLE when no assumption is made on P * . Two algorithms have been introduced to compute it. [START_REF] Mclachlan | A maximum likelihood method for random coefficient regression models[END_REF] exploited a duality with the D-optimal design theory to develop the Nonparametric Maximum Likelihood (NPML) algorithm. [START_REF] Schumitzky | Nonparametric EM algorithms for estimating prior distributions[END_REF] proposed a generalization of the EM algorithm (Dempster, [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]: the Nonparametric EM (NPEM) algorithm. An advanced version of the NPEM algorithm, named NPAG, has been developed by the USC Laboratory of Applied Pharmacokinetics (USC Pack, 2002). [START_REF] Gallant | Semi-nonparametric maximum likelihood estimation[END_REF] and Davidian & Gallant (1993) investigated a Smooth (or Semi) Nonparametric (SNP) estimator assuming that P * admits a smooth probability density function π * with respect to the Lebesgue measure. [START_REF] Kuhn | Estimation par maximum de vraisemblance dans des problèmes inverses non linéaires[END_REF] proposed a consistent logspline estimator of π * . Finally, since 2006, NONMEM has introduced a discrete estimator whose support points are the EBEs [START_REF] Boeckmann | NONMEM Users Guide[END_REF]. In the present article, we will refer to it as NPNM which stands for NP in NONMEM.

Several articles evaluated one of these NP methods on simulated data: [START_REF] Roe | Comparison of population pharmacokinetic modeling using simulated data: results from the population modeling workgroup[END_REF] studied NPML; [START_REF] Bustad | Parametric and Non parametric Population Methods : Their Comparative Performance in Analysing Clinical Dataset and Two MonteCarlo Simulation Studies[END_REF] evaluated NPEM; [START_REF] Savic | Evaluation of the nonparametric method in NONMEM V Iβ[END_REF] focused on NPNM. However, there is a lack of a comparison between the different NP methods. Besides these studies do not include comparison with the empirical distribution of EBEs which remains widely used, despite its weaknesses. This article offers to carry out a detailed review and comparison of four NP methods among the most widely used or documented ones: NPML, NPEM, SNP and NPNM, with focus on PK applications. The empirical distribution of EBEs is also included in the study. The logspline estimator [START_REF] Kuhn | Estimation par maximum de vraisemblance dans des problèmes inverses non linéaires[END_REF] is not further investigated here, because it is only defined for unidimensional random effects whereas the simplest PK model deals with bidimensional random effects.

The article is organized as follows. The model and the notations are presented in section 2. Section 3 describes the different methods studied. Section 4 is dedicated to the review of their asymptotic properties. Section 5 presents two PK simulation studies that empirically compare the performance of the methods. The article ends with a conclusion in section 6.

Notations and framework

We consider the following general NLME model:

Y i = f i (d i , θ * , X i ) + Γ i (d i , θ * , X i )ε i (1)
where

-Y i = (Y ij , j = 1..n i ) is the n i -vector of observations on individual i.
-d i is the known experimental design for individual i.

-θ * is an unknown p-vector of fixed parameters (or fixed effects).

-X i = (X 1i , ..., X qi ) is the q-vector of real random effects associated with individual i. The X i s are assumed independent and identically distributed (iid) from the probability measure P * :

X i iid ∼ P * i = 1..N.
We denote Π * the cumulative distribution function (cdf) of X i , and, if it exists, π * its probability density function (pdf) with respect to the Lebesgue measure. -f i is a known real function depending on i, d i , θ * and X i . An important feature of NLME models is the nonlinearity of f i with respect to the random effects

X i . -ε i = (ε ij , j = 1..n i ) is a white noise: ε i iid ∼ N (0, I n i ) i = 1..N. ε i is independent of X i . -Γ i (d i , θ * , X i ) is a known positive low triangular n i × n i -matrix depending on d i , θ * and X i .
The challenge of the statistical inference is the estimation of θ * and P * from the observation of Y i , i = 1..N while the X i s are unobserved.

In phase II or phase III clinical trials, the number of individual measurements (n i ) is usually not large enough to be under asymptotic conditions and can even be very low. Thus, the only asymptotic examined in the following will be the number N of individuals in the sample, whereas n i will remain finite and possibly small (with possible differences between individuals).

It is important to note that, in our NP context, the first moment of P * can be infinite. Consequently, it cannot be considered as a fixed effect as usually done in parametric models. One goal of the PK and/or PD analyses is to explain the mean variations of X i with covariates (age, body weight, renal function...). Obviously in this case, it is necessary to assume that the first moment of P * is finite. The ability of the NP methods to apply in this context will be discussed in the following.

3 Studied methods

The EBEs' empirical distribution (EBEsED)

In population PK and/or PD analyses, it is common to assume the (log-) normality of X i . In this parametric model, P * only depends on its mean and variance. These parameters and the fixed effects θ * are generally estimated by MLE: ( θparametric

N , Pparametric N ) = argsup θ,P∈Nq L N (θ, P),
where N q denotes the set of the q-dimensional (log-)normal distributions and L N (θ, P) the likelihood of θ and P.

After this parametric estimation step, some predictions of the random effects X i , known as empirical Bayes estimates (EBEs), can be computed:

EBE i = argsup x P[X i = x|Y i , Pparametric N , θparametric N ],
where

P[X i |Y i , Pparametric N , θparametric N ] is the pdf of X i conditionally to Y i , θparametric N and Pparametric N .
In practice, the (log-) normality of X i is often checked by looking at the EBEs' empirical distribution (EBEsED). In that context, EBEsED can be considered as an NP estimator of the random effects distribution P * : it is discrete with EBEs as support points and equal frequencies for all individuals (equal to 1/N ). Its cdf ΠEBEsED N is:

ΠEBEsED N (x) = 1 N N i=1 1 [EBE i x] ,
where 1 [EBE i x] equals 1 when, for all j = 1..q, EBE i (j) x(j), and 0 otherwise.

The drawback of EBEsED is that it gives the same weight (1/N ) to all EBEs even if the accuracy of each EBE i can be very different depending on the number of measurements per individual (n i ) and on the quality of the experimental design (d i ) (i.e. the choice of the measurement times).

Nonparametric Maximum Likelihood (NPML)

When no assumption is made on P * , the maximum likelihood estimator (MLE) P N P N is defined as:

PNP N = argsup P∈Pq L N (θ, P) = argsup P∈Pq N i=1 P(Y i |X i = x, θ)P(dx)
where P q denotes the space of all probability measures on R q and P(Y i |X i = x, θ) the density of Y i conditionally to X i = x with θ for fixed effects. The following theorem gives an handy way to compute PNP N .

Theorem 1 (Discreteness of PNP N [START_REF] Lindsay | The geometry of mixture likelihoods : a general theory[END_REF])) Assuming that

{(P[Y i |X i = x, θ], i = 1..N ), x ∈ R q } is
a compact set, P N P N can be expressed as a discrete measure with at most N support points.

The proof is elementary. It deals with convexity and the theorem of Caratheodory. From a computational point of view, this discreteness has two useful consequences. First, for a given N , the likelihood maximization can be restricted to the finite dimensional set of all discrete measures with at most N support points. Secondly, the integral with respect to PNP N is a finite sum, and thus the likelihood L N (θ, PNP N ) has an explicit form:

L N (θ, PNP N ) = N i=1 P[Y i |X i = x, θ]d PNP N (x) = N i=1 N l=1 P[Y i |X i = X l , θ] × p l ,
where X l , l = 1..N denote the support points of PNP N and p l , l = 1..N their probabilities.

Despite the explicit writing of the likelihood, its maximum is still not explicit. Therefore, iterative algorithms have been developed to compute PNP N . [START_REF] Mclachlan | A maximum likelihood method for random coefficient regression models[END_REF] used an optimal design analogy to build the NPML algorithm.

NPML is derived from the Fedorov algorithm [START_REF] Fedorov | Theory of Optimal Experiments[END_REF]

L+1 is defined as:

X (k+1) L+1 = argsup X ∈R q ∂L N (θ, P(k) N ) ∂X ,
where

∂L N (θ, P(k) N ) ∂X
is the directional derivative of the likelihood in the direction X . If argmax , .., p

X ∈R q ∂L N (θ, P(k) N ) ∂X = 0,
(k+1) L+1 ) = argsup L N (θ, P(k+1) N ),
where the supremum is over all p l , l = 1..L + 1 between 0 and 1 such that L+1 l=1 p l = 1. 3. Index limitation step. If P(k+1)

N has N + 1 support points, a point can be removed from the support without a decreasing of the likelihood. Such a point X can be determined by the resolution of the linear system:

P Y |X × X = b
where P Y |X is an (N + 1) × (N + 1)-matrix and b an (N + 1)-vector:

     P Y |X ij = P[Y i |X i = X (k+1) j ] i = 1..N, j = 1..N + 1, P Y |X ij = 1 i = N + 1, j = 1..N + 1. If P Y |X is singular then b i = 0 ∀i = 1..N + 1, else b i = 1 P[Y i ] i = 1..N and b N +1 = 0.
Theorem 2 (convergence of NPML, [START_REF] Mclachlan | A maximum likelihood method for random coefficient regression models[END_REF]) Each iteration of NPML increases the sample likelihood. Besides, if the maximum likelihood estimator PNP N is unique:

P(k) N -→ k→+∞ PNP N
where P(k) N is the NPML estimator after k iterations.

Nonparametric Expectation Maximization (NPEM)

As [START_REF] Mclachlan | A maximum likelihood method for random coefficient regression models[END_REF], [START_REF] Schumitzky | Nonparametric EM algorithms for estimating prior distributions[END_REF] proposed an algorithm to compute the MLE P N P N defined in the preceding section.

The Expectation Maximisation (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] is acknowledged to be very efficient for likelihood maximization with unobserved data. However, despite a wide area of applications, this algorithm does not permit to compute PNP N . Actually, the EM algorithm is restricted to absolutely continuous estimators of P * . [START_REF] Schumitzky | Nonparametric EM algorithms for estimating prior distributions[END_REF] extended the EM algorithm to the specific computation of P N P N . The Nonparametric EM (NPEM) algorithm consists of 2 steps. At iteration (k + 1), given

P(k) N with L support points (L N ), P (k+1)
N is defined in two steps:

1. Support optimization step. For l = 1..L, the support points X (k+1) l of P(k+1)

N are computed as:

X (k+1) l = argsup X ∈R q N i=1 P[X (k) l |Y i , P(k) N , θ] × ln P[Y i |X i = X , θ],
where X

l , l = 1..L are the support points of P(k) N . 2. Frequencies optimization step. For l = 1..L, the probabilities p (k+1) l of X (k+1) l have an explicit form:

p (k+1) l = 1 N N i=1 P[X (k) l |Y i , P(k) N , θ]
Theorem 3 (Convergence of NPEM, [START_REF] Schumitzky | Nonparametric EM algorithms for estimating prior distributions[END_REF]) Each iteration of NPEM increases the likelihood.

NPEM often leads to measures with a low number of support points. Actually, NPEM does not allow to increase the number of support points, but some points can be removed if their probability becomes zero.

Smooth Nonparametric (SNP)

In most applications, the random effects X i have absolutely continuous distributions. Davidian & Gallant (1993) proposed a smooth estimation of the density π * of X i . More precisely, they considered the set SN P K of all density functions with finite mean m of form:

[P K (U -1 (x -m))] 2 × N q (x -m; 0, U U T ) + π 0 (x -m),
where -K is a fixed truncation parameter, -N q (x; M, V ) is the density of a q-dimensional normal random variable, with mean M and variance matrix V , -P K (x) is a polynomial in x = (x 1 , ..., x q ) of degree K, -U is a non singular upper triangular matrix, U T is the transposed matrix of U, -is small positive real number, -π 0 is a strictly positive probability density function with zero mean, and such that π 0 S is bounded, where . S is a weighted Sobolev norm (see [START_REF] Gallant | Semi-nonparametric maximum likelihood estimation[END_REF] for details).

The term π 0 (x -m) ensures that all the densities considered are strictly positive. In practice, this lower bound can generally be disregarded ( = 0) without affecting the behavior of the log-likelihood.

The SNP estimators θSNP N and πSNP N of θ * and π * are defined as the MLE over the set SN P K :

( θSNP N , πSNP N ) = argsup θ,π∈SN P K L N (θ, π(x)dx)
The likelihood is not explicit, but it can be approximated by Gauss-Hermite quadrature. The maximization of the approximated likelihood with respect to θ, the coefficients of the polynomial P K and the terms of U under the constraints πSNP N (x)dx = 1 and xπ SN P N (x)dx = m can be performed thanks to standard algorithms. Davidian & Gallant (1993) suggested for instance the Sequential Quadratic Programming algorithm.

Nonparametric in NONMEM (NPNM)

The software NONMEM (version VI) has recently introduced a nonparametric subroutine. As EBEsED, this estimator, referred to as NPNM, is a discrete measure with the EBEs used as support points. However, in contrast to EBEsED, the frequencies p 1 , .., p N associated with support points EBE 1 , .., EBE N are not necessarily the same for all individuals. They maximize the log-likelihood:

(p 1 , .., p N ) = argsup ln L N ( θparametric N , PNPNM N ) = argsup N i=1 ln( N l=1 p l P[Y i |X i = EBE l , θparametric N ]),
under the constraint 0 p l 1 for all l = 1..N and N l=1 p l = 1. It shall be noticed that, thanks to the discreteness of NPNM, the log-likelihood is explicit. Besides, this optimization is a convex problem.

Asymptotic behaviors

The EBEs' empirical distribution (EBEsED)

In this section, the consistency of EBEsED, considered as an NP estimator of P * , is investigated. If all X i (i = 1..N ) are observed, their empirical cdf is an unbiased consistent estimator of the cdf Π * . In general mixed effects models, X i are not observed but EBE i , defined in section 3, can give an accurate prediction of X i . Besides, if θ * is known and all n i are large, the accuracy of this prediction is little affected by a misspecification of P * (see [START_REF] Van Der Vaart | Asymptotic Statistics[END_REF] for instance): conditionally on

X i = x i , EBE i -x i a.s. -→ n i →+∞ 0 i = 1..N even if P * is misspecified.
Consequently, when θ * is known and all n i are large, EBEsED is consistent in the sense that, if for all i, n i ∼ +∞,

ΠEBEsED N (x) = 1 N N i=1 1 EBE i x a.s. -→ N →+∞ Π * (x), ∀x ∈ R q .
However, when data are sparse, the properties of EBEsED can be questioned.

In a general NLME model, when n i remains fixed, the rigorous study of the asymptotic (in N ) behavior of EBEsED appears to be a hard matter. However, for a specific linear homoscedastic mixed effects model, the following theorem shows that EBEsED is not consistent.

Theorem 4 (non consistency of EBEsED) Consider the model:

Y i = X i + θ * ε i , i = 1..N,
with the same notations as for model 1. Besides, we suppose that n i is the same for all individuals i. Under the assumptions:

1. |E[X i ]| < ∞ and var[X i ] < ∞, 2. ∃M < ∞, lim N →+∞ sup i=1..N |EBE i | M .
If -→ means weak convergence, we have:

PEBEsED N -→ N →+∞ P * .
Proof In this linear mixed effects model, the EBE i , PEBEsED N , its mean and its variance are explicit. They depend on the parametric estimators of θ * , E[X i ] and var[X i ], which, in this linear mixed effects model, are consistent even if P * is misspecified. Then, by the law of large numbers, assumption 1, the continuous mapping theorem, one can easily show that:

E[ PEBEsED N ] -→ N →+∞ E[P * ] and var[ PEBEsED N ] -→ N →+∞ 1 1 + 1 n i θ * 2 var[P * ] var[P * ] = var[P * ]
By assumption 2, it implies the non weak consistency of PEBEsED N (van der Vaart, 1998, example 2.21).

The shrinkage phenomenon, described by [START_REF] Savic | Un)informativeness of empirical Bayes estimate-based diagnostics[END_REF], is an illustration of this result: when the data are sparse (that is the number of measurement per individual is small) or 'uninformative' (the variance of the residual error is high in comparison with the variance of X i ), the variance of EBEsED shrinks, i.e. it is smaller than the variance of P * . As expected, when individual information is important (n i is large), the bias of var[ PEBEsED N ] is insignificant. To our knowledge, the (non) consistency of EBEsED has not been established in sparse general NLME models.

NPML and NPEM

NPML and NPEM aim at computing the same estimator, PNP N , defined as:

PNP N = argsup P∈Pq L N (θ, P).
Here P is unrestricted since P q denotes the set of all probability measures on R q .

For a specific NLME model, [START_REF] Pfanzagl | Consistency of maximum likelihood estimators for certain nonparametric families, in particular: mixtures[END_REF] gave the asymptotic properties of PNP N . His result has been extended to a non-equally distributed (n i and d i are not the same for all individuals) and heteroscedastic NLME model by [START_REF] Chafaï | On Nonparametric Maximum Likelihood for a class of stochastic inverse problems[END_REF].

Theorem 5 (Consistency of PNP N [START_REF] Chafaï | On Nonparametric Maximum Likelihood for a class of stochastic inverse problems[END_REF])) Under the assumptions 1. ∀(P 1 , P 2 ) ∈ P q × P q , P 1 = P 2 ⇒ L N (θ, P 1 )dy 1 ..dy N = L N (θ, P 2 )dy 1 ..dy N in P N × P N , 2. the fixed effects θ * are known.

The maximum likelihood estimator PNP N is weakly consistent:

PNP N -→ N →+∞ P * .
Assumption 1 ensures identifiability of the model. Assumption 2 is very restrictive because θ * is generally unknown. [START_REF] Pfanzagl | Large deviation probabilities for certain nonparametric maximum likelihood estimators[END_REF] showed that assumption 2 can be relaxed if θ * is consistently estimated whatever P * . But, to the best of our knowledge, such an estimation of θ * is not currently available for a general NLME model. [START_REF] Lai | Nonparametric estimation in nonlinear mixed effects models[END_REF] proposed to jointly estimate θ * and P * by likelihood maximization. This procedure is consistent in an Ibragimov-Has'minski framework:

∃I, 1 I N, I -→ N →+∞ +∞ such that ∀i = 1..I, n i ∼ +∞.
As noticed in section 2, it is usual in PK and/or PD analyses to constrain the mean of X i (E[X i ] = m < ∞) in order to explain between subjects variations with the available covariates. Unfortunately, to the best of our knowledge, generalization of theorem 5 when the mean of X i is constrained is not trivial. [START_REF] Mallet | Handling Covariates in Population Pharmacokinetics, with an Application to Gentamicin[END_REF] circumvent this difficulty by considering all covariates as random effects. In this case, the distribution of interest is the distribution of X i conditionally to covariates. The joint distribution of X i and the covariates is first estimated by unconstrained likelihood maximization (with NPML or NPEM). The distribution of X i conditionally to covariates is then easily obtained by the Bayes formula. The method provides a consistent estimator of the distribution of X i conditionally to covariates. Besides, the specification of a regression model between X i and the covariates, and estimation of the regression parameters, are no longer necessary. On the other hand, this method increases the dimension of the distribution to be estimated. This highly increases the computation time and requires very large samples because the estimation of a joint distribution is much more demanding than the estimation of the distribution of X i .

Smooth Nonparametric (SNP)

The estimators θSNP

( θSNP N , πSNP N ) = argmax (θ,π)∈Θ×SN P K L N (θ, π(x)dx),
where Θ ⊂ R p is compact and SN P K is the set of smooth density functions defined in section 3.4. Under the assumptions:

1. E[X i ] = m is finite. The actual density π * of X i can be written as:

π * (x) = h 2 (x -m) + h 0 (x -m),
where is a positive real number, h and h 0 are probability density functions with zero mean. h S and h 0 S are bounded, where . S denotes a weighted Sobolev norm, see [START_REF] Gallant | Semi-nonparametric maximum likelihood estimation[END_REF]. h 0 is strictly positive. The actual fixed effect θ * lives in Θ.

2. K -→ N →+∞ +∞ 3. there is a function L(θ, π, θ * , π * ) that is continuous in (θ, π) with respect to (θ, π) = ( θ 2 + π 2 S∞ ) 1/2 such that, with probability 1, lim N →∞ sup θ×π |L(Y i , i = 1..N |π(x)dx, θ) -L(θ, π, θ * , π * )| = 0
where . S∞ is an infinite weighted Sobolev norm (see [START_REF] Gallant | Semi-nonparametric maximum likelihood estimation[END_REF] for details). 4. the model is identifiable.

The following convergences hold:

πSNP N -π * S∞ a.s. -→ N →+∞ 0 and θSNP N -θ * a.s. -→ N →+∞ 0
Consistency with respect to . S∞ implies L1-consistency of πSNP N and of its derivatives up to a certain order, and consistency of the first moments of X i .

Assumption 1 imposes regularity conditions on π: the bound on h S restricts violent oscillatory behavior and tail's thickness, conditions on h 0 bound the tails from below. This assumption seems realistic for most of the applications because numerous type of densities, like skewed, leptokurtic, platykurtic, or multi-modal, are allowed. Assumption 2 ensures that the potential complexity of πSNP N increases with the sample size. The truncation parameter K thus plays a key role in the SNP estimation. In practice, its determination for a given sample size N is not obvious. [START_REF] Fenton | Convergence rates of SNP density estimators[END_REF] investigated deterministic rules. When X i is observed and univariate, they found that the optimal convergence rate (with respect to the L1-norm) was achieved for K = N 1 5 . To our knowledge, there has been no extension of this result to NLME models. In these models, Davidian & Gallant (1993) recommend to 'inspect plots [...] for all models between those selected by Schwarz and Akaike criteria inclusively and make a visual selection'. An automatic rule can be necessary when an a priori knowledge or the subjectivity of the analyst could distort the figures interpretation. In that case, Davidian & Gallant (1993) advise to use the Hannan-Quinn criterion which gives good results in others Hermite series expansion studies.

It is noteworthy that theorem 6 proves the consistency of the fixed effects' estimation jointly with the distribution of random effects' estimation. Such a result is not established for the others NP methods.

It shall be noticed that theorem 6 ensures that πSNP N and θSNP N remain consistent when the mean of X i is constrained: this allows to use covariates for explaining the mean variations of X i .

Given that the (log-) normal distribution is a special case of SNP estimation (for K = 0), the (log)-normality assumption can be tested with, for instance, the likelihood ratio test.

Nonparametric in NONMEM (NPNM)

To the best of our knowledge, asymptotic properties of NPNM have not been documented.

Comparison on simulation studies

In order to compare the various NP methods in practical situations, we performed two simulation studies. This section is organized in 3 parts: part 1 describes the design of the simulations, part 2 presents the various criterion used to compare the methods' performance and part 3 is devoted the results.

Design of simulations

General strategy

The simulations were conducted in a PK framework. In order to simulate realistic datasets, the model's settings were chosen according to real population PK studies of phenobarbital in infants. The phenobarbital is a barbiturate generally used to treat neonatal seizures. Two simulations were performed: simulation 1 illustrated a simple PK analysis, simulation 2 illustrates a more challenging PK analysis.

Simulation 1 was inspired by the phenobarbital PK analysis of [START_REF] Grasela | Neonatal population pharmacokinetics of phenobarbital derived from routine clinical data[END_REF]. The PK of phenobarbital after intravenous administration is described with a one compartment model with first order elimination, as follows: for

j = 1..n i , [f i (d i , θ * , X i )] j = D X i2 × exp -X i1 ×t ij /X i2 .
The known experimental design d i consists of the administered dose D at time 0 and the measurement times t ij . The 2 random effects are: the clearance (X i1 ) and the volume of distribution (X i2 ). The measurement error is assumed to be proportional to the expected concentration:

Γ i (d i , θ * , X i ) = θ * × diag(f i (d i , θ * , X i )),
where θ * is the only fixed effect. The clearance was simulated with a bimodal distribution. Clearance and volume of distribution were considered independent. In order to investigate the influence of data sparseness, a rich and a sparse experimental design were simulated (in PK, a dataset is usually considered rich if, for each individual i, n i q, q being the number of random effects).

Simulation 2 was inspired by the phenobarbital PK analysis of [START_REF] Yukawa | Population pharmacokinetic investigation of phenobarbital by mixed effect modelling using routine clinical pharmacokinetic data in Japanese neonates and infants[END_REF]. The PK of phenobarbital after oral administration is described with a one compartment model with first-order elimination and first-order absorption, as follows: for j = 1..n i

[f i (d i , θ * , X i )] j = DX i3 X i2 (X i3 -X i1 /X i2 ) (exp -X i1 ×t ij /X i2 -exp -X i3 ×t ij ).
The known experimental design d i consists of the administered dose D at time 0 and the measurement times t ij . The 3 random effects are: the apparent clearance (X i1 ), the apparent volume of distribution (X i2 ) and the absorption rate constant (X i3 ). The measurement error is assumed to be proportional to the expected concentration. The apparent clearance was simulated with a bimodal distribution. Apparent clearance and volume of distribution were considered correlated.

For each simulation, sample of different sizes were simulated (N = 50, 100, 200, 300, 400) in order to assess the consistency of the methods.

For each simulation and sample size, 100 independent datasets were simulated and analyzed with the various methods studied (EBEsED, NPML, NPEM, SNP, NPNM).

More details on simulations 1, on simulation 2 and on computational settings are provided in the following sections.

Details on simulation 1

In order to choose the values of θ * and of the first moments of X i , we analyzed the real dataset [START_REF] Grasela | Neonatal population pharmacokinetics of phenobarbital derived from routine clinical data[END_REF], without taking into account the covariates and assuming that X i is log-normal. The estimation was performed using NONMEM software (subroutine ADVAN1 TRANS2, method FOCEI). The estimated values were:

-12.8% for θ * , -0.00613 liters per hour (L/h) for the mean clearance and 34.3% × 0.00613 for its standard deviation, -1.60 liters (L) for the mean volume of distribution and 47.3% × 1.60 for its standard deviation.

We simulated X i with the following distribution with a minor (30%) and a major (70%) sub-population (see figure 1):

X i1 iid ∼     
ln N (0.00735, (14.84% × 0.00735) 2 ) with probability 70% ln N (0.00330, (20.27% × 0.00330) 2 ) with probability 30% and X i2

iid

∼ ln N (1.60, (47.32% × 1.60) 2 ),
where ln N (m, v) denotes the log-normal distribution with mean m and variance v, X i1 and X i2 were independently simulated. The values of modes and proportions were arbitrary chosen in such a way that the overall mean and variance were similar to those estimated from the real dataset used by [START_REF] Grasela | Neonatal population pharmacokinetics of phenobarbital derived from routine clinical data[END_REF]. In the simulations, we supposed that a single dose D was administered at time t = 0 to all individuals. The value of D was set to the mean of the initial doses of the real dataset used by [START_REF] Grasela | Neonatal population pharmacokinetics of phenobarbital derived from routine clinical data[END_REF]: D = 25.7µg. We simulated some rich and some sparse datasets. For rich datasets, we simulated:

-2 observations per individual for 85% of the individuals: one early and one late, -3 observations per individual for 15% of the individuals: one early, one intermediate and one late.

For sparse datasets, we simulated:

-1 single observation per individual for 85% of the individuals. This observation was early for half of the individuals, late for others individuals, -3 observations per individual for 15% of the individuals: one early, one intermediate and one late.

The observation times were: 12 minutes (min) for the early observation, 6 hours (h) for the intermediate one and 246 h for the late one. These times respectively correspond to the minimum, the median and the maximum of the observation times following intravenous administration in the real dataset.

Details on simulation 2

The value of θ * used for the simulations equals the one found by [START_REF] Yukawa | Population pharmacokinetic investigation of phenobarbital by mixed effect modelling using routine clinical pharmacokinetic data in Japanese neonates and infants[END_REF]: 25.2% for oral administration. [START_REF] Yukawa | Population pharmacokinetic investigation of phenobarbital by mixed effect modelling using routine clinical pharmacokinetic data in Japanese neonates and infants[END_REF] also found that, for an individual with a mean body weight of 2.9 kg and a mean postnatal age of 21 days: the apparent clearance has a mean of 0.0363 L/h and a standard deviation of 32%×0.0363 L/h, the apparent volume of distribution has a mean of 7.75 L and a standard deviation of 54% × 7.75 L, the absorption rate constant equals 50h -1 (it was considered as a known fixed effect).

We simulated X i with the following distribution (see figure 2):

X i1 iid ∼     
ln N (0.0465, (43.01% × 0.0465) 2 ) with probability 70% ln N (0.0124, (47.72% × 0.0124) 2 ) with probability 30%,

X i2
iid ∼ ln N (7.75, 17.5) and X i3 iid ∼ ln N (5.00, 2.25).

Within each sub-population, the apparent clearance and volume of distribution were correlated: correlation[ln(X i1 ), ln(X i2 )] = 0.6.

The absorption rate constant (X i3 ) was independent of apparent clearance (X i1 ) and volume of distribution (X i2 ). The values of modes and proportions were arbitrary chosen in such a way that: the mean of X i1 and X i2 and the standard deviation of X i2 were similar to the value found by [START_REF] Yukawa | Population pharmacokinetic investigation of phenobarbital by mixed effect modelling using routine clinical pharmacokinetic data in Japanese neonates and infants[END_REF]. The X i1 coefficient of variation was increased from 32% to 64% to obtain a larger between individuals variability, very usual in PK applications. The absorption rate constant was considered as a random effect in order to have a more complex model. Its mean value was decreased from 50h -1 to 5h -1 .

The value of the dose administered (D) was set to the mean of the real dataset's doses: D = 12, 400µg. In contrast to simulation 1, we only simulated sparse datasets, with: Figure 2. The distribution of random effects chosen in simulation 2. On the left, contour plot of the joint pdf of the apparent clearance and the apparent volume of distribution at quantiles 5%, 10%, 25%, 50%, 75%, 90% and 95%. The pdf of apparent clearance is bimodal. The apparent clearance and volume of distribution are correlated. On the right, the absorption rate constant's marginal pdf: it was independently simulated as a log-normal random variable.

-2 observations per individual for 85% of the individuals. One of these observations was around the maximum of the concentration-time profile. The other observation was early for half of the individuals, and late for others individuals. -4 observations per individual for 15% of the individuals: one early, one intermediate, one late and one around the maximum of the concentrationtime profile.

The early observation time was equal to 12 min, the intermediate one to 6 h and the late one to 246 h. The time around the maximum of the concentrationtime profile was set to 1h 24 min (this time achieved the maximum of the population mean concentration-time profile).

Computational settings

In order to compute EBEsED,NPNM and an estimation θparametric N of θ * , we assumed that P * was a log-normal distribution and used the software NONMEM VI with the $ESTIMATION (FOCEI method for simulation 1, FO method for simulation 2), option POSTHOC and $NONPARAMETRIC routines [START_REF] Boeckmann | NONMEM Users Guide[END_REF]). FO has been used in simulation 2 because FOCEI failed too often.

In order to perform NPML and NPEM, an estimation of θ * is needed: we chose θparametric N . NPML and NPEM algorithms were implemented in C ++ . NPNM was used as starting point for NPML and NPEM. A Fletcher-Reeves algorithm with 10 random different initializations was used to perform the support optimization step of NPML. The frequencies optimization step of NPML was performed with an EM algorithm as described by [START_REF] Mclachlan | A maximum likelihood method for random coefficient regression models[END_REF]. NPEM was stopped when:

ln L N ( θparametric N , P(k+1) ) -ln L N ( θparametric N , P(k) ) | ln L N ( θparametric N , P(k) )| < 10 -4 ,
where ln L N ( θparametric N , P(k) ) is the NPEM likelihood at iteration k. This stopping criteria was not convenient for NPML, because the likelihood improvement between 2 NPML iterations is always very small (each iteration of NPML modifies only 1 support point whereas NPEM modifies N support points). We considered the likelihood improvement between N iterations: NPML was stopped when:

ln L N ( θparametric N , P(k+N) ) -ln L N ( θparametric N , P(k) ) | ln L N ( θparametric N , P(k) )| < 10 -4 ,
where ln L N ( θparametric

N , P(k) ) is the NPML likelihood at iteration k.
A fortran 77 implementation of the SNP method, called nlmix [START_REF] Davidian | Nlmix: A Program for Maximum Likelihood Estimation of Nonlinear Mixed Effects Model with a Smooth Random Effects Density[END_REF], is available in the public domain (statlib, Carnegie-Mellon University). The SNP estimations of θ * , θSNP N , and of P * were performed with nlmix. The likelihood was computed by quadrature (with 15 quadrature points in simulation 1 and, to save computation time, only 10 in simulation 2). Following Davidian & Gallant (1993) and [START_REF] Davidian | Nlmix: A Program for Maximum Likelihood Estimation of Nonlinear Mixed Effects Model with a Smooth Random Effects Density[END_REF] recommendations, the matrix U was set diagonal for K > 0. Besides, for each dataset, several truncation parameters K (K=0,2,3 and 4) and, for each K, several initialization points (10 in simulation 1, 5 in simulation 2) were assessed. For a given K, we selected the initialization which performed the highest likelihood. The best K was then selected using the Hannan-Quinn criterion. Besides, when K was greater than 0, the matrix U was set diagonal.

Performance evaluation

Comparing the estimation of a distribution in inverse problems is not an straightforward task. Indeed, ways of comparing estimated and actual probability distribution of random effects are numerous, and in general not equivalent. Besides, in such inverse problems the objective in not always the estimation of the random effects' distribution. For instance in this study, the ability to detect the sub-population was of specific interest. Similarly, when the aim of the statistical study is more predictive than descriptive, the ability to predict the observed data Y i or the unobserved random effects X i could be a pertinent information. Therefore, it is necessary to use several criteria to compare the methods.

Estimation of P *

To compare the estimated and the true distribution, the T1 distance (also known as Wasserstein or Kantorovich or transport distance) is convenient here because it can be computed for multidimensional discrete and continuous distributions. The T1-error made by PN (with cdf ΠN ) in estimating P * (with cdf Π * ) is defined as: T1-error( PN ) = | ΠN (x) -Π * (x)|dx, T1-error represents the absolute volume between the estimated and the true cdf.

Prediction of random effects

A prediction of unobserved random effects X i can be computed after an estimation PN of P * and θN of θ * :

EBE i ( θN , PN ) = arg sup X i P[X i |Y i , θN , PN ].
It is noteworthy that EBE i ( θparametric N , EBEsED) equals the parametric EBE i defined in section 3.1. For each dataset, the sum of absolute relative errors between EBE i ( θN , PN ) and the actual realization x i of X i , gives the average percentage error (denoted X-error) for the prediction of X i :

X-error p ( θN , PN ) = 1 N N i=1 | EBE i ( θN , PN ) p -x ip x ip |.

Prediction of observed data

In order to evaluate the errors made in predicting the individual concentration time profile of interest, we computed, for each dataset:

Y-error( θN , PN ) = 1 N N i=1 +∞ t=0 |f (t, EBE i ( θN , PN )) -f (t, x i )|dt,
where f (t, x) = f i (d, θ, x), where d represents the experimental design with dose D and measurement time t i = (t). The Y-error represents the average absolute area between the predicted and the real concentration-time profiles.

Results

Computation

The NP methods were only run when the NONMEM $ESTIMATION routine converged without warnings. Table 1 reports the number of these datasets. On these datasets, all NP methods converged (SNP was considered convergent if at least one initialization converged). for SNP) and of the time needed to compute the estimation of P * . For SNP, the computation time consists of the computation time needed to perform the MLE with K = 0, 2, 3 and 4, and for each K several initializations (10 for simulation 1, 5 for simulation 2). As expected, the computation time of all methods increased with the sample size and with the number of random effects. SNP appeared particularly sensitive to an increase of the random effects dimension: this can be explained by the use of quadratures to approximate the likelihood.

Estimation of P * For rich datasets (simulation 1 only), EBEsED performed as well as the NP methods. However, for sparse datasets (for simulations 1 and 2), it produced much more T1-error than the NP methods and its T1-error did not really decrease with the sample size: it seemed not consistent with respect to T1 metric.

On the contrary, in all situations, the T1-errors of all the NP methods de- creased closed to 0 with the sample size: all the NP methods seemed convergent with respect to T1 metric. For simulation 1, the NP methods performed roughly equivalently with respect to the T1-error. As expected, these errors were lower for the rich than for the sparse datasets. It is noteworthy that, in this simulation, NPNM performed nearly as well as the others methods, which is surprising as its support points suffer from shrinkage. For simulation 2, the difference between the NP methods was more perceptible. In particular, NPNM appeared less successful than NPML and NPEM: we suppose that, because of shrinkage, the EBEs did not provide a convenient support to describe P * in this 3-dimensional model. The performance of SNP was not as good as in the 2-dimensional model: we speculate that the 10 quadrature points were not enough to give an accurate approximation of the likelihood and/or that the 5 initializations assessed were not sufficient to reach the global maximum in this high dimensional optimization problem (e.g. for K = 4, the dimension equals 41). The difference between NPML and NPEM was very slight. q q q q q q q 50 100 200 300 400 0.00 0.02 0.04 0.06 0.08 0.10 sample size N q q 0.00 0.02 0.04 0.06 0.08 0.10 q q q 0.00 0.02 0.04 0.06 0.08 0.10 q q q q q q q 0.00 0.02 0.04 0.06 0.08 0.10 q q q 0.00 0.02 0.04 0.06 0.08 0.10 EBEsED NPML NPEM SNP NPNM q q q 50 100 200 300 400 0.00 0.02 0.04 0.06 0.08 0.10 sample size N q q q q q q q q 0.00 0.02 0.04 0.06 0.08 0.10 q q q q q q q q q q 0.00 0.02 0.04 0.06 0.08 0.10 q q q q q q q 0.00 0.02 0.04 0.06 0.08 0.10 q q q q q q 0.00 0.02 0.04 0.06 0.08 0.10 Figure 4. For simulation 1, boxplots of T1-errors for various methods with respect to the sample size N : on the left for the rich datasets, on the right for the sparse datasets.

q q q q 50 100 200 300 400 0 1 2 3 4 sample size N q q q q q q q q q q 0 1 2 3 4 q q q q q q q 0 1 2 3 4 q q q q 0 1 2 3 4 q q q q q q q 0 1 2

Maximization of the likelihood

In order to evaluate the performances of the algorithms for likelihood maximization, we looked at the likelihood obtained with the various methods. The likelihood of EBEsED considered here is not the approximated likelihood given by NONMEM software after the parametric estimation, but the (explicit) likelihood of the NP discrete measure with EBEs as support points and equal frequencies. In that sense, it is comparable with the likelihood of NPNM, NPML and NPEM. The likelihood of SNP is not considered here since it is not comparable with the others because the estimation of θ * is dif-ferent. Table 2 reports the average log-likelihood of the various methods for samples of size N = 400. Results for the other sample sizes are not presented here but led to the same conclusions. In theory, for a given estimation of θ * , the likelihood of NPML and NPEM should be greater than the likelihood of EBEsED and NPNM, since these algorithms maximize the likelihood over the set of all probability measures. Besides by definition, the likelihood of NPNM cannot be less than the EBEsED's one. Our results are in agreement with the theory for all the datasets. For simulation 1, the likelihood improvement due to the optimization of support points (difference between NPML/NPEM and NPNM) is slight in comparison with the likelihood increasing due to the frequencies optimization (the difference between NPNM and EBEsED), in both sparse and rich datasets. On the contrary, for simulation 2, the optimization of the support points provided a real improvement of likelihood. Overall, these results are coherent with the T1-error results presented above. In all situations, in average, NPEM reached a higher likelihood than NPML, with a lower computational cost. Therefore, it can be concluded that the NPEM algorithm was, in average, more efficient than the NPML algorithm, under the studied conditions.

Prediction abilities

Figure 6 displays the boxplots of prediction errors for the clearance (X-error 1 ), for simulation 1, the sparse datasets and samples of size N = 400. In the overall population, EBEsED produced roughly as much X-error 1 as the NP methods, despite its poor estimation of P * . Actually, when looking at the 2 sub-populations separately, we can see that: in the major sub-population, EBEsED seemed better than the discrete NP methods, whereas in the minor

Conclusion

The goal of this article was to compare several NP methods for estimation of the random effects' distribution P * in NLME models, with applications to PK analyses. We studied four NP methods: NPNM, NPML, NPEM and SNP.

EBEsED which is often used to check the parametric assumption, has also been included in the study.

As expected, EBEsED appeared to be an efficient estimator of P * when the experimental design was rich but was unreliable when the experimental design was sparse. In that case, the need for NP methods was confirmed.

The NPML and NPEM algorithms aim at computing the MLE when no assumption is made on P * . This MLE is a discrete measure with at most N support points. This discreteness has the advantage of making the likelihood explicit. On the other hand, it can make the interpretation difficult since P * is generally continuous. The consistency of this MLE is established, but under restrictive assumptions: all fixed effects are known and no covariates are used to explain mean variations of X i . In our simulation studies, NPML and NPEM gave satisfactory results even if the fixed effect (relative to residual error) was previously estimated using a parametric model. In general, NPEM reached a higher likelihood with a lower computational cost. Both NPML and NPEM present the risk of being trapped in local maxima since they are deterministic algorithms.

In contrast to these discrete NP methods, the SNP method provides a continuous estimator of P * . Its statistical properties are established under mild assumptions. In particular, the fixed effects can be estimated and covariates can be used to explain the mean variations of X i . Besides, the (log-)normality can be tested. SNP have two main weaknesses: the selection of the truncation parameter and the likelihood optimization. To select the truncation parameter, no criterion has well established properties in NLME models. As for the likelihood optimization, this is difficult for two reasons. First, the likelihood is not explicit: here, it was approximated with quadratures which is computationally very expensive when the number of random effects is important. Second, the likelihood certainly has several local maxima: the algorithm used here was deterministic and so several initializations were needed to avoid local maxima. Investigations on stochastic algorithms would certainly be fruitful. Lastly, NPNM is an easy-to-compute discrete estimator. On a simple PK model, it performed nearly as good as the others NP estimators, but was less satisfactory in a more complex PK model. Since its support points can be affected by shrinkage, its statistical properties are questionable. Finite mixtures of (log-) normal distributions could also be an alternative to the parametric assumption (Lemenuel-Diot, [START_REF] Lemenuel-Diot | Estimating Heterogeneity in Random Effects Models for Longitudinal Data[END_REF], [START_REF] Wang | Nonlinear random effects mixture models: Maximum likelihood estimation via the EM algorithm[END_REF]). If P * is really a finite mixture, the MLE is consistent and some computational algorithms are available. If not, one could think that increasing the number of mixture components would allow to approximate any kind of distribution, but, to our knowledge, such a result is not established.

In conclusion, all the studied NP methods seemed promising but would probably require more investigations on statistical and/or computational issues in order to be more widely used in population PK and/or PD analyses.

)

  be the MLE defined as:

Figure 1 .

 1 Figure1. The distribution of random effects chosen in simulation 1: contour plots of the joint pdf of clearance and volume of distribution at quantiles 5%, 10%, 25%, 50%, 75%, 90% and 95%. The pdf of clearance is bimodal. Clearance and volume of distribution are independent.

Figure 3

 3 Figure 3 displays the average computation times. The computation time of each method consists of the computation time needed to compute the fixed effects estimation ( θparametric N

Figure 4

 4 Figure4and 5 display boxplots of the T1-errors. For rich datasets (simulation 1 only), EBEsED performed as well as the NP methods. However, for sparse datasets (for simulations 1 and 2), it produced much more T1-error than the NP methods and its T1-error did not really decrease with the sample size: it seemed not consistent with respect to T1 metric. On the contrary, in all situations, the T1-errors of all the NP methods de-

Figure 3 .

 3 Figure3. Average computation times with respect to the sample size N : on top, for simulation 1, on the left, for rich datasets, on the right for sparse datasets, at the bottom, for simulation 2 (for sparse datasets only).

Figure 5 .

 5 Figure 5. For simulation 2, boxplots of T1-errors with respect to the sample size N (the datasets are sparse).

Table 2

 2 Average log-likelihood of the various discrete NP estimators, for samples of size N = 400. * the likelihood of EBEsED computed here is not the parametric likelihood but the NP likelihood considering EBEsED as a discrete NP estimator with EBEs as support points and equals frequencies.

	method	Simulation 1	Simulation 1	Simulation 2
		Rich datasets Sparse datasets Sparse datasets
	EBEsED *	-2403.3	-1518.6	-7260.7
	NPNM	-2377.9	-1450.5	-7235.9
	NPML	-2376.7	-1450.1	-7183.4
	NPEM	-2376.3	-1449.1	-7166.9

sub-population, the discrete NP methods seemed better than EBEsED. For all methods, the X-error 1 was more important in the minor than in the major sub-population. Surprisingly, SNP performed equivalently to EBEsED with respect to X-error 1 , even in the minor sub-population. q 0.0 0.2 0.4 0.6 0.8 1.0 q 0.0 0.2 0.4 0.6 0.8 1.0 q 0.0 0.2 0.4 0.6 0.8 1.0 q q q 0.0 0.2 0.4 0.6 0.8 1.0 q 0.0 0.2 0.4 0.6 0.8 1.0 EBEsED NPML NPEM SNP NPNM q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Figure 6. Boxplots of prediction errors for the clearance (X-error 1 ), for simulation 1, sparse datasets and samples of size N = 400: on the left for the overall population, in the middle within the major sub-population (around 70% of the individuals), on the right within the minor sub-population.

Figure 7 displays boxplots of prediction errors (Y-error), for simulation 1, the sparse datasets and samples of size N = 400. In the overall population, EBEsED made less Y-error than the discrete NP methods. Within the major sub-population, EBEsED was better than all the NP methods. However, within the minor sub-population, all the NP methods appeared better than EBEsED. It is noteworthy that SNP was the best in this sub-population, despite its poor prediction of clearance (X-error 1 ). q 0 1000 3000 5000 q q q 0 1000 3000 5000 q q 0 1000 3000 5000 q q 0 1000 3000 5000 q q q 0 1000 3000 5000 EBEsED NPML NPEM SNP NPNM 0 1000 3000 5000 q q q 0 1000 3000 5000 q 0 1000 3000 5000 q q 0 1000 3000 5000 q 0 1000 3000 5000 q 0 1000 3000 5000 0 1000 3000 5000 q q 0 1000 3000 5000 q q q q q 0 1000 3000 5000 q q 0 1000 3000 5000