Emmanuelle Comets
email: emmanuelle.comets@inserm.fr

Christelle Rodrigues

Vincent Jullien

Moreno Ursino

Conditional non-parametric bootstrap for non-linear mixed effect models

Keywords: Bootstrap, Conditional distribution, Non-linear mixed effect models, Uncertainty of parameter estimates Conditional

Purpose: Non-linear mixed effect models are widely used and increasingly integrated into decision-making processes. Propagating uncertainty is an important element of this process, and while standard errors (SE) on parameters are most often computed using asymptotic approaches, alternative methods such as the bootstrap are also available. In this article, we propose a modified residual parametric bootstrap taking into account the different levels of variability involved in these models.

Methods: The proposed approach uses samples from the individual conditional distribution, and was implemented in R using the saemix algorithm. We performed a simulation study to assess its performance in different scenarios, comparing it to the asymptotic approximation and to standard bootstraps in terms of coverage, also looking at bias in the parameters and their SE.

Results: Simulations with an Emax model with different designs and sigmoidicity factors showed a similar coverage rate to the parametric bootstrap, while requiring less hypotheses. Bootstrap improved coverage in several scenarios compared to the asymptotic method especially for the variance param-

Introduction

Nonlinear mixed effect models (NLMEM) are widely used namely in pharmacokinetics to characterise the evolution of drug concentrations after dose administration, and in pharmacodynamics where they are used to integrate disease models and pathophysiology to better understand the response of biological systems to treatment [START_REF] Lalonde | Model-based drug development[END_REF]. Modelling and simulation has become a key component in the analysis and design of clinical trials to establish the most promising therapeutic candidates, as well as the dosage regimens more likely to be effective and safe [START_REF] Wang | Model-informed drug development: current US regulatory practice and future considerations[END_REF], and is also increasingly used to individualise diagnostics and treatments [START_REF] Desmée | Nonlinear joint models for individual dynamic prediction of risk of death using Hamiltonian Monte Carlo: application to metastatic prostate cancer[END_REF]. As the models become increasingly more comprehensive, and the extrapolations increasingly more sophisticated, it becomes also important to report the confidence we have in them to propagate this uncertainty into decision making processes and realistic evaluation. In the present article, we focus on the issue of standard errors of estimation (SE), which measure the precision of estimation for the parameters in the model. In maximum-likelihod estimation, they are often computed using an asymptotic approximation as the inverse of the Fisher information matrix [START_REF] Retout | Fisher information matrix for nonlinear mixed-effects models: evaluation and application for optimal design of enoxaparin population pharmacokinetics[END_REF], while Bayesian methods consider the entire posterior distributions and use it to assess and propagate uncertainty [START_REF] Gelman | Bayesian Data Analysis[END_REF]. How well the asymptotic approximation holds in practice depends on the sample size, both in terms of number of subjects and in terms of number of samples per subject, and on the degree of non-linearity in the model [START_REF] Dartois | Evaluation of uncertainty parameters estimated by different population PK software and methods[END_REF]. Other approaches to derive the standard errors include log-likelihood profiling [START_REF] Lindstrom | Nonlinear mixed effects models for repeated measures data[END_REF], Sampling Importance Resampling (SIR) [START_REF] Dosne | dOFV distributions: a new diagnostic for the adequacy of parameter uncertainty in nonlinear mixedeffects models applied to the bootstrap[END_REF] and bootstrap methods [START_REF] Thai | A comparison of bootstrap approaches for estimating uncertainty of parameters in linear mixed-effects models[END_REF].

Bootstrap methods mimic the process of repeating a study several times to assess the distribution of an estimator, using only the data available in the original study through resampling. For independent and identically distributed observations, Efron [START_REF] Efron | Bootstrap methods: another look at the jacknife[END_REF] proposed resampling with replacement from the pool of available measurements, creating bootstrap datasets on which the estimator of interest is applied in order to approximate its distribution. In more complex settings, like linear regression or non-homoscedastic error models, issues such as within-subject correlation and heteroscedastic variances have to be considered. Several bootstrap approaches have been proposed in these situations, with the most intuitive being the case bootstrap, also called paired bootstrap [START_REF] Kümmel | Confidence and Prediction Intervals forPharmacometric Models[END_REF], where the resampling is performed at the level of the subject.

Alternatively, the residual bootstrap resamples residuals in the empirical distribution (non-parametric residual bootstrap) or in an external distribution with suitable characteristics (parametric boostrap) [START_REF] Efron | An introduction to the bootstrap[END_REF]. Bootstraps were extended to NLMEM by considering the different levels of variability involved at the individual and population level. Das and Krishen [START_REF] Das | Some bootstrap methods in non-linear mixed-effects models[END_REF] proposed a residual bootstrap involving two levels of resampling. In small sample studies, however, the distribution of the residuals does not reflect the variability in the population because of a phenomenon called shrinkage, causing a regression to the mean due to the limited amount of information to jointly estimate a large number of parameters and/or the residual variability. Carpenter et al. [START_REF] Carpenter | A novel bootstrap procedure for assessing the relationship between class size and achievement[END_REF] proposed to correct the distribution of residuals before resampling, using a ratio between the empirical variance-covariance matrix and the estimated variance matrix of the random effects. Different bootstrap methods were summarised and evaluated both in linear and non-linear mixed effect models by Thai et al. [START_REF] Thai | A comparison of bootstrap approaches for estimating uncertainty of parameters in linear mixed-effects models[END_REF][START_REF] Thai | Evaluation of bootstrap methods for estimating uncertainty of parameters in nonlinear mixed-effects models: a simulation study in population pharmacokinetics[END_REF]. For linear mixed effect models with homoscedastic variance, the case, parametric and non-parametric bootstraps exhibited good performances in terms of absence of bias and coverage rates, while the asymptotic method proved reliable for the estimation of the parameters and their SE [START_REF] Thai | A comparison of bootstrap approaches for estimating uncertainty of parameters in linear mixed-effects models[END_REF]. These bootstraps were then evaluated in NLMEM, where the results were more contrasted. In this more challenging settings, all bootstraps and the asymptotic method exhibited some measure of bias, especially on a model involving non-linear pharmacokinetics where some parameters have strong non-linearity [START_REF] Thai | Evaluation of bootstrap methods for estimating uncertainty of parameters in nonlinear mixed-effects models: a simulation study in population pharmacokinetics[END_REF]).

The objective of the present manuscript was to develop a new bootstrap approach and evaluate it by comparison to usual bootstraps in NLMEM. Our proposed approach, the conditional non-parametric bootstrap, is a residual bootstrap using samples from the the conditional distribution of the individual parameters, which have been shown to correct the shrinkage in the empirical Bayes estimates in non-linear models [START_REF] Lavielle | Enhanced method for diagnosing pharmacometric models: random sampling from conditional distributions[END_REF]. We use simulations inspired by the study performed by Plan et al. [START_REF] Plan | Performance comparison of various maximum likelihood nonlinear mixed-effects estimation methods for dose-response models[END_REF], using an Emax and a Hill model to investigate the influence of non-linearity with both sparse and rich designs.

Methods

Statistical models

Let the random variable y ij denote the observation of the longitudinal data, which we will assume in this work to be continuous, at time t ij for subject i = 1, ..., N (j = 1, ..., n i ). Following [START_REF] Lavielle | Mixed effects models for the population approach: models, tasks, methods and tools[END_REF], the NLMEM can be formalised as a hierarchical model involving a parametric structural model f defining the longitudinal evolution, depending on fixed parameters µ and subject-specific random effects η i with a parametric distribution, coupled with a probability distribution on the measurements, here a Gaussian with a standard deviation defined by a parametric function g:

         y ij = f (x ij , θ i ) + g(x ij , θ i , σ) ǫ ij θ i = h(µ, η i ) η i ∼ N (0, Ω) ǫ ij ∼ N (0, 1) (1)
In these equations, θ i denotes the individual parameters, which we assume are related through a function h to the p-dimensional vector containing the fixed effects µ, and the q-dimensional vector containing the random effects η i . ǫ ij is a random variable assumed to be normally distributed. σ denotes variance parameters entering the function g, which expresses the standard deviation of the measurement error and is generally either constant (homoscedastic variance) or a function of f . The random effects η i and the residual errors ǫ ij are assumed to be independent for different subjects and to be independent of each other for the same subject. We will denote by ξ i = {x i1 , ..., x ini } the design variables, for example the n i sampling times for subject i.

The parameters of the model Ψ = (µ, Ω, σ) can be estimated by maximum likelihood. In NLMEM, the likelihood associated with equation ( 1) is intractable as individual likelihoods need to integrate out the unknown parameters θ i over their distribution D θ . In this study, we used the SAEM algorithm to obtain parameter estimates [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF][START_REF] Lavielle | Mixed effects models for the population approach: models, tasks, methods and tools[END_REF]. The maximum likelihood estimate (MLE) θ of θ is asymptotically normally distributed with mean θ and asymptotic covariance matrix given by the inverse of the Fisher information matrix (FIM), denoted M F , usually obtained through a first-order approximation of the model [START_REF] Retout | Fisher information matrix for nonlinear mixed-effects models: evaluation and application for optimal design of enoxaparin population pharmacokinetics[END_REF]. During the estimation process, the individual parameters θ i act as nuisance parameters, but their conditional distribution is used during the stochastic approximation phase to sample the unknown parameters η i and obtain a complete dataset from which the conditional log-likelihood is derived [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF]. The estimated individual parameters, called empirical Bayes estimates (EBE), can be defined as the mode or the median of the conditional distribution obtained at the end of the estimation process.

Conditional bootstrap

Bootstrap consists in repeatedly generating pseudo-samples distributed according to the same distribution as that the original sample. In NLMEM, bootstrap methods include the case bootstrap (Case), resampling the entire vector of individual observations, the parametric bootstrap resampling residuals from a theoretical distribution (Par) and the non-parametric bootstrap resampling from estimated residuals (NP) [START_REF] Thai | Evaluation of bootstrap methods for estimating uncertainty of parameters in nonlinear mixed-effects models: a simulation study in population pharmacokinetics[END_REF]. We propose a new bootstrap method, which we call conditional non-parametric residual bootstrap (cNP), using samples from the conditional distributions instead of shrinkage-corrected empirical Bayes estimates.

The algorithm is as follows: after fitting the model to the data, we obtain K samples for each vector η i , η ik (k = 1, ..., K) from the conditional distri-butions, center them and resample vectors of individual parameters η * i with replacement. In this work we implement this step by sampling from all conditional samples over i = 1...N, k = 1...K but we could also resample for each individual within his or her conditional samples. Error residuals are sampled from the centered and shrinkage-corrected residuals as in NP bootstrap. More details on bootstrap algorithms and implementation are given in supplementary material.

Simulation study

Models:

The model used for this simulation is a Hill (or sigmoid E max ) model, a standard model in dose-response studies where the effect of a drug in response to a dose d, E(d) follows:

E(d) = E 0 + E max d γ i d γ i + ED γ 50 (2)
It involves 4 parameters, the initial effect E 0 , the maximum effect E max , the concentration at which half the maximum effect is achieved EC 50 and the sigmoidicity factor γ which controls the non-linearity of the model through the curvature. Interindividual variability was modelled through a log-normal distribution for all parameters, except for γ which was assumed to be the same for all subjects (no IIV). A correlation was simulated between E max and EC 50 , and we used a proportional error model with coefficient of variation σ. The parameters used in the simulation were the same as those used by Plan et al. [START_REF] Plan | Performance comparison of various maximum likelihood nonlinear mixed-effects estimation methods for dose-response models[END_REF], and are given in Table 1.

Table 1 Parameters used in the simulation for the Emax model (see [START_REF] Plan | Performance comparison of various maximum likelihood nonlinear mixed-effects estimation methods for dose-response models[END_REF]).

Parameter

Scenarios:

We investigated different scenarios to assess the influence of nonlinearity, with γ = 1 and 3, and design (Table 2). For each value of γ we considered two designs, a rich design with N=100 subjects given 4 doses each (0, 100, 300, 1000), and a sparse design with N=200 subjects divided in 4 groups of 50 subjects given 2 doses each among the following 4 combinations (0,1000),(100,1000),(0,300) and (100,300). To evaluate the performance of bootstraps with a smaller sample size, the two scenarios with a Hill model were also evaluated with N=20 subjects (rich design) and 20 subjects per group (sparse design). For each scenario, we simulated K=200 datasets. Bootstrap distributions were obtained with B=100 bootstrap samples for each bootstrap. Additional simulations with up to B=500 bootstrap samples were used to assess the number of bootstrap samples necessary.

Bootstrap evaluation: The bootstraps were evaluated in terms of coverage rate, and we also computed the bootstrap mean θB,k , the bootstrap standard error (SE) and the bootstrap confidence interval (CI) for each estimated parameter from the bootstrap distributions. The bootstrap parameters estimates and their SE were defined as unbiased with relative error within ±5%, moderately biased (from ±5% to ±10%) and strongly biased (> ±10%). The coverage rate of the 95% CI was considered to be good (from 90% to 100%), low (from 80% to 90%) and poor (< 80%). A good bootstrap was defined as a method without bias for parameter estimates and their corresponding SE, and ensuring a good coverage rate of the 95% CI.

All computations and simulations were performed in R [START_REF] Development | R: A Language and Environment for Statistical Computing[END_REF] using the development version of the saemix library [START_REF] Comets | Parameter estimation in nonlinear mixed effect models using saemix, an R implementation of the SAEM algorithm[END_REF]. Details of the computations are given in Supplementary Material.

Results

Bootstrap distributions

As an illustration, Figure 1 shows the distribution of the parameters obtained for one randomly selected simulated dataset obtained with the different bootstraps in scenario S Hill,R (Appendix B shows example of simulated data). For each parameter, we plot in red the true value (used in the simulation) and in blue the value estimated by saemix (for E 0 and ω Emax the two values are almost identical so the lines are overlayed). In this example the distributions are for the most part similar, but we notice some discrepancy for ω Emax where the NP bootstrap is shifted to the left, and in some cases the distribution appears more peaked, usually with Par or with Case.

Bootstrap performances

Figure 2 shows the total relative bias on parameter estimates and their SE for the 6 scenarios simulated in this study. Dashed lines delineate absolute relative biases within 5% and dotted lines absolue relative biases within 10%. To make the graph easier to read, each line of the plot, representing one scenario, has its own scale on the Y-axis because of the difference in information between the different designs. The results from the estimation by saemix on each dataset are shown in grey, and the SE in this case are obtained using the asymptotic FIM, and the bootstrap distributions in different colours. The plots in the left column, representing the relative bias on the fixed and random parameters in the model, show a similar pattern across the different bootstraps: all the fixed parameters are estimated without bias across all scenarios except S Hill,smallS which had the least information. Concerning random effects, the covariance term seems most sensitive to a sparse design. For the rich designs (scenarios S Emax,R and S Hill,R ), the asymptotic method underestimated the SE for the random effects in S Emax,R and for nearly all the parameters in S Hill,R with biases as high as -30%. With the rich design however, the bootstraps manage to compensate some of the bias and reduce it to the 10% range. Additional results suggest that biased bootstrap estimates are associated with biased estimates on the original dataset (Appendix C).

With sparse designs (S Emax,S and S Hill,S ), although the bootstraps remain biased for the estimation of the SE, they correct in part the high bias of the asymptotic method. Finally, the last two lines show the scenarios with a small number of subjects, where all methods fail to estimate properly the parameters and their SE, although there again most of the bias on the parameters is in fact due to poor estimates by saemix.

Comparing the conditional non-parametric bootstrap (cNP) with the nonparametric bootstrap (NP), we note that cNP behaves generally similarly to NP in the first three scenarios in terms of parameter estimates, while it performs slightly better in S Hill,S (sparse design, higher model non-linearity) and in the scenarios with less subjects except for the estimation of the covariance term in scenario S Hill,smallR . The estimates of the SEs are also on the whole slightly more accurate with cNP. In Figure 2 we used the standard deviation of the bootstrap distribution to estimate the SE, however the bootstrap is a non-parametric approach designed to produce the distribution of the estimator, and this distribution may not be symmetrically distributed. Therefore in Figure 3 we plot the coverage for each parameter in the 6 scenarios (rich designs in the left column and sparse designs in the right column), using the 4 bootstrap distributions and the asymptotic confidence intervals. This plot shows good coverage for the case, parametric, and conditional non-parametric bootstraps in the first four scenarios, while the asymptotic method struggles with the non-linearity in the Hill model (see also Appendix D and E). cNP performs better than NP in all scenarios investigated, with coverage rates closer to the nominal value (a comparison between the coverage rates for these two approaches is given in Figure S10, Appendix I). Similar bias and coverage rate were found for σ in both bootstrap approaches, which is consistent with using the same correction in the cNP as in the NP bootstrap. To complement the coverage itself, Appendix F shows a figure showing the percentage of cases in which the true parameter is either lower or higher than the confidence interval. This figure indicates that apart from very few cases, lower coverage is due to the confidence interval being too low with respect to the true value of the parameter. We also note in S Hill,S and S Hill,smallS that except for E 0 , the bootstraps for fixed effects tend to yield coverages close to 100%, suggesting overly wide confidence intervals. To investigate this further we also computed 50% CI (see Appendix G) and found that indeed, coverages could go as high as 70% in particular in S Hill,S and S Hill,smallS for these parameters, with none of the bootstraps appearing to have a clear advantage.

Concerning runtimes, to obtain 100 bootstrap samples for one dataset from the SHill,R scenario, parametric and standard residual bootstraps took 5% longer than the case bootstrap, while cNP took 40% longer. Most of the additional burden for the cNP bootstrap is however due to the step involving the estimation of the conditional distribution in the current, non optimised, version of the code, and the conditional distribution is normally estimated to perform diagnostic graphs after an saemix run. Finally, additional results concerning the influence of the number of bootstrap samples on the performances are given in Appendix H.

Discussion

The aim of the present article was to propose a new bootstrap approach for non-linear mixed effect models, and evaluate it by comparison to three bootstrap methods previously recommended [START_REF] Thai | Evaluation of bootstrap methods for estimating uncertainty of parameters in nonlinear mixed-effects models: a simulation study in population pharmacokinetics[END_REF]. We use samples from the conditional distributions proposed by Lavielle and Ribba [START_REF] Lavielle | Enhanced method for diagnosing pharmacometric models: random sampling from conditional distributions[END_REF] to bypass the issue of shrinkage for non-parametric bootstrap.

In the present article, we use a simple simulation example based on a sigmoid E max model as a proof of concept to evaluate this new bootstrap, on the assumption that bootstrap methods should provide good results in simple cases first if they are to be trusted with more complex situations. In our simulations, the conditional non-parametric bootstrap showed good coverage in different scenarios, and performed better than the standard non-parametric bootstrap. Its performance was comparable to the case and parametric bootstraps, which showed good performance in a simulation setting with a homogenous and informative experimental design, in line with previous studies [START_REF] Thai | A comparison of bootstrap approaches for estimating uncertainty of parameters in linear mixed-effects models[END_REF][START_REF] Thai | Evaluation of bootstrap methods for estimating uncertainty of parameters in nonlinear mixed-effects models: a simulation study in population pharmacokinetics[END_REF]. A question arising at the time was whether it would be possible to improve the correction of the residuals before resampling, to benefit from the non-parametric bootstrap's ability to maintain the structure of the original dataset, while being less dependent on model assumptions than the parametric bootstrap. That question echoed similar interrogations around the usefulness of graphical and statistical diagnostics for NLMEM [START_REF] Karlsson | Diagnosing model diagnostics[END_REF][START_REF] Combes | Powers of the likelihood ratio test and the correlation test using empirical bayes estimates for various shrinkages in population pharmacokinetics[END_REF] in designs with low information where the empirical Bayes estimates exhibit high shrinkage. Recently, Lavielle and Ribba [START_REF] Lavielle | Enhanced method for diagnosing pharmacometric models: random sampling from conditional distributions[END_REF] proposed new diagnostic tools which use the conditional distribution of individual parameters, obtained as a product of stochastic estimation methods that have been gaining increasing popularity in the field of pharmacometrics and beyond. These conditional distributions allow to bypass the issue of shrinkage, and it was therefore a natural step to consider using these distributions to sample from. The individual conditional distributions are a by-product of the estimation in the SAEM algorithm, and the mode or mean of which are used as estimates of the individual parameters. With linearisationbased algorithms, these individual parameters are obtained by minimising a Bayesian criterion, but additional computations could be used to derive the same distributions. Sampling individual parameters from the individual conditional distributions p(ψ i |y i ) generates a sample of parameters representative of the population distribution p(ψ i ) [START_REF] Lavielle | Enhanced method for diagnosing pharmacometric models: random sampling from conditional distributions[END_REF]. Diagnostic graphs using these samples, for example of parameters versus covariates, do not suffer from shrinkage. When we use these resampled parameters to build the bootstrapped datasets in the cNP bootstrap, we therefore expect to generate datasets representative of the true interindividual variability. In unbalanced designs, it might also be worth considering resampling the random effects within each subjects conditional distributions; indeed, the posterior distributions capture the uncertainty on individual parameter estimates, reflecting the informativeness of the individual sampling design which differs according to the number and position of the samples. As such, our cNP bootstrap retains the attractive feature of the NP bootstrap of preserving the structure of the original dataset with an unchanged number of samples and subjects, and combines it with a more appropriate distribution for the random effects. An extension for future work will therefore be to investigate unbalanced designs, where some subjects may have more information than others, such as the design mixing rich and sparse data in [START_REF] Thai | Evaluation of bootstrap methods for estimating uncertainty of parameters in nonlinear mixed-effects models: a simulation study in population pharmacokinetics[END_REF], where the non-parametric bootstrap greatly overestimated the SE for some variance parameters but a stratified version led to a much lower coverage. By using conditional distributions instead of theoretical distributions as in the parametric bootstrap, we could also expect the cNP to be more robust to misspecifications in the random effects. This will however need to be tested in a further study, as we show here that bootstraps cannot recover from very poor asymptotic estimates. Also, conditional distributions involve model assumptions on the distribution of the random effects, and their properties have not been studied fully yet, especially in small samples.

For the residual errors ǫ , we used the same correction in cNP as in the NP bootstrap, inflating the residuals so that their overall variability matches the population estimate of the residual error. However, this is not entirely satisfactory as the ǫij are also affected by shrinkage, called ǫ-shrinkage by Karlsson and Savic [START_REF] Karlsson | Diagnosing model diagnostics[END_REF], and a global level correction incurs similar penalties as for η-shrinkage. This could be the reason why cNP and the NP bootstraps gave similar results for σ. In the present study, the designs were homogenous in terms of informativeness, but it would be interesting to evaluate the difference between cNP and NP, in unbalanced designs. A further extension to the cNP will be to compute the residual errors for the conditional samples drawn from the individual conditional distribution of the parameters, and resampling from these to reflect the individual residual variability.

Although the bootstrap is often branded as an approach to gain a good understanding of the uncertainty in small samples, Broeker and Wicha have shown that this is mostly wishful thinking as the performance in small datasets is not good [START_REF] Broeker | Assessing parameter uncertainty in small-n pharmacometric analyses: value of the log-likelihood profiling-based sampling importance resampling (LLP-SIR) technique[END_REF]. We confirm this here in the two scenarios S Hill,smallR and S Hill,smallS . With even lower number of subjects, especially in the sparse design, the bias on the parameter estimates for the asymptotic method was even more apparent (data not shown) and the bootstrap approaches could not recover from these poor estimates. With these scenarios the coverage for all bootstraps was lower than 90%, and there was no clear trend in whether the bootstraps under-or over-estimated the parameters and their SE, which is a cautionary message for the usage of bootstrap in practice as from present and previous work [START_REF] Thai | Evaluation of bootstrap methods for estimating uncertainty of parameters in nonlinear mixed-effects models: a simulation study in population pharmacokinetics[END_REF][START_REF] Broeker | Assessing parameter uncertainty in small-n pharmacometric analyses: value of the log-likelihood profiling-based sampling importance resampling (LLP-SIR) technique[END_REF] the conditions in which bootstrap performances are adequate are not obvious. In our simulations, we found that the performance of the bootstraps degrade as the informativeness of the design decrease. Tools like PFIM [START_REF] Dumont | PFIM 4.0, an extended R program for design evaluation and optimization in nonlinear mixed-effect models[END_REF] can be used to evaluate the design before running a bootstrap approach, to assess whether a bootstrap procedure is appropriate.

A final point is the high computational effort incurred. Even the simplest version, the case bootstrap, requires a large number of runs; its theoretical shortcomings, such as the difference in information between bootstrapped datasets in unbalanced designs, or the difficulty in handling covariate distributions when they enter a model, have caused other even more computationally intensive bootstrap algorithms to be developed [START_REF] Thai | A comparison of bootstrap approaches for estimating uncertainty of parameters in linear mixed-effects models[END_REF]. These bootstraps require repeated model evaluations on top of the sampling step, making them somewhat slower. However, these approaches only need to be run for key models such as the final model in an analysis. On top of this, an additional source of uncertainty comes from the model itself, as a model is rarely an isolated object but an element in model building; methods like model averaging or sensitivity analyses should also be used to gain a full understanding of uncertainty [START_REF] Buatois | Comparison of model averaging and model selection in dose finding trials analyzed by nonlinear mixed effect models[END_REF].

Conclusion

Our new conditional residual bootstrap for NLMEM improved the coverage rate compared to the classical non-parametric residual bootstrap, proving nearly as efficient as the parametric bootstrap, which in these simulations where we used the true distributions for the parametric bootstrap can be considered as a gold standard. The case bootstrap also performed well, and all three bootstraps provided more appropriate estimates of the SE than the asymptotic method in sparse designs. The bootstrap estimates of the parameters and their SE showed bias in the presence of non-linearity and sparse information, highlighting the importance of considering the entire bootstrap distribution for uncertainty instead of using a normal approximation with mean and standard deviation derived from the distribution.

Fig. 1

 1 Fig. 1 Bootstrap distributions obtained for one randomly sampled distribution in scenario S Hill,R . The red line in each plot represents the true value of the parameter and the blue line the value estimated by saemix.

Fig. 2

 2 Fig.2Relative bias on parameters (left) and SE (right) for the 6 scenarios simulated in this study. Dotted lines delineate absolute relative biases within 5% and dashed lines absolue relative biases within 10%. Each line of plot has its own scale on the Y-axis.

Fig. 3

 3 Fig. 3 Coverage rates for the 6 scenarios simulated in this study. Dotted lines indicate a coverage of 90 and 100%.

Table 2

 2 Simulation scenarios.

	Simulation scenarios	Value of γ	Design Number of subjects
	S Emax,R	1	rich	100
	S Emax,S	1	sparse	200
	S Hill,R	3	rich	100
	S Hill,S	3	sparse	200
	S Hill,smallR	3	rich	20
	S Hill,smallS	3	sparse	80

Funding

This work was performed during the authors' research time and no other funding was used.

Data Availability

All data generated or analysed during this study are included in this published article and its supplementary information files (see Appendix A) and are available on Zenodo (https://doi.org/10.5281/zenodo.4059718).

Conflict of interest

None.

Author Contributions

Conceptualisation