Search tips
Search criteria 


Logo of aapsjspringer.comThis journalToc AlertsSubmit OnlineOpen Choice
AAPS J. 2009 September; 11(3): 615.
Published online 2009 August 28. doi:  10.1208/s12248-009-9138-8
PMCID: PMC2758131

Evaluation of an Extended Grid Method for Estimation Using Nonparametric Distributions


A nonparametric population method with support points from the empirical Bayes estimates (EBE) has recently been introduced (default method). However, EBE distribution may, with sparse and small datasets, not provide a suitable range of support points. This study aims to develop a method based on a prior parametric analysis capable of providing a nonparametric grid with adequate support points range. A new method extends the nonparametric grid with additional support points generated by simulation from the parametric distribution, hence the name extended-grid method. The joint probability density function is estimated at the extended grid. The performance of the new method was evaluated and compared to the default method via Monte Carlo simulations using simple IV bolus model and sparse (200 subject, two samples per subject) or small (30 subjects, three samples per subjects) datasets and two scenarios based on real case studies. Parameter distributions estimated by the default and the extended-grid method were compared to the true distributions; bias and precision were assessed at different percentiles. With small datasets, the bias was similar between methods (<10%); however, precision was markedly improved with the new method (by 43%). With sparse datasets, both bias (from 5.9% to 3%) and precision (by 60%) were improved. For simulated scenarios based on real study designs, extended-grid predictions were in a good agreement with true values. A new approach to obtain support points for the nonparametric method has been developed, and it displayed good estimation properties. The extended-grid method is automated, using the program PsN, for implementation into the NONMEM.

Key words: empirical Bayes estimates, extended grid method, NONMEM, nonparametric method, parameter distribution


A major objective and challenge of any nonlinear mixed-effect (“population”) analysis is to accurately quantify and separate different sources of variability in model parameters describing drug pharmacokinetics and pharmacodynamics. This is not an easy task and pharmacometricians often experience difficulties quantifying all variability terms and/or parameter correlations, which may lead to an inefficient model building process and the selection of a suboptimal model. Moreover, certain assumptions are involved in the most commonly used parametric estimation process, such as: (1) the parameters in the populations are lognormally distributed (or symmetrically distributed on log-scale) or (2) the assumption that parameter correlations are linear. Although sometimes these assumptions may seem to be reasonable, often they do not hold true and there is no method to prove otherwise. This could lead into numerous problems: biased results, difficulty in detecting model misspecification, failure of models to perform realistic simulations, or the interference of inadequate assumptions with the estimation procedure itself (1,2). In such a case, a more general methodology, known as the nonparametric nonlinear mixed effects modeling may be an alternative.

Nonparametric population methods are attractive and increasingly used tools in pharmacometric analysis as they relax the parameter distribution assumption. These methods estimate entire parameter distributions as a collection of parameter values (support points) and associated population probability. However, from a numerical point of view, these methods are computer intensive for models with many parameters and for datasets with a large number of subjects (1,35). Long runtimes can be efficiently minimized due to the high parallelization efficiency of these nonparametric algorithms. Recently, a new nonparametric approach has been introduced in NONMEM VI, which is a maximum likelihood method and, like any other nonparametric method, involves two steps: (1) estimation of support points for a nonparametric distribution and (2) estimation of population probabilities associated with each support point. What differentiates this method from other nonparametric methods is the first step, which involves empirical Bayes estimation of points of support. Although the empirical Bayes estimation is still a parametric method of describing individual parameter estimates, this innovative shortcut has proven to be beneficial for the subsequent nonparametric estimation from several points of view: (1) the empirical Bayes estimates (EBEs) is a best point estimate of individual parameters which can be obtained from the parametric estimation process (6), and the EBEs contain useful information from the preceding parametric analysis (7); thus, they also serve as a valuable starting point for the subsequent nonparametric analysis; (2) use of EBEs helps to shorten the usually long analysis time of the first step in nonparametric estimation; (3) this approach has been recently implemented into the software NONMEM, which is known for its data format and model specification flexibility. Following the release of this method, the idea of using EBEs as starting support points for a subsequent nonparametric motivated other method developers, and a new algorithm based on this principle has been suggested (8). It is important to underline that the nonparametric method in NONMEM evaluates only one set of support points, the EBEs, and therefore it is a stationary method, as opposed to other nonparametric methods which are iterative procedures allowing points of support to move to other locations. Sophisticated algorithms have been suggested which allow points of support to spread out in space by a mathematically optimal method (4).

However, EBEs do have shortcomings as they suffer from the “shrinkage phenomenon” in certain circumstances (9), and therefore they may not always provide a suitable range of support points. Also, it is well known that both parametric and nonparametric methods have drawbacks when the number of subjects in the datasets is low. For parametric methods, a small sample size prevents precise estimation of variance and proper assessment of the distribution assumption, while, for nonparametric analysis, the methods available estimate nonparametric distribution at the maximal number of support points equal to number of individuals in the datasets. In this circumstance, the discrete nature of nonparametric methods is often questioned for simulation purposes of novel scenarios, as simulations from the model only can represent as many unique individuals as were present in the original study.

The aim of this paper is to: (1) develop a novel nonparametric method that has improved estimation properties when data are sparse and sample size is small, (2) evaluate this novel method via Monte Carlo (MC) simulation studies followed by re-estimation using a range of different scenarios, and (3) make the method available for users by implementation into the software PsN (10).


The extended-grid method improves the density of the nonparametric grid by addition of new support points. New support points are simulated from the parametric distribution with the mean and variance estimated in the preceding parametric analysis. Similar alternative methods of generating points of support include using inflated/deflated variances, simulations from other distributions (e.g., uniform), sampling from the parameter space after dividing it into equiprobable subspaces, etc. These extensions will however not be investigated in this work apart from simulations using inflated variance. The extended-grid method does not have a required upper limit for the number of support points, and it also may or may not keep the EBEs from the preceding parametric step in the nonparametric grid. The former will be the method used in the present work. The entire nonparametric joint probability density function is estimated at the new extended grid. This method leads to a nonparametric distribution estimated at a number of support points larger than the number of individuals in the datasets, and thus it may address both imprecision due to small sample size and imprecise and/or shrunk EBEs.

The performance of the extended-grid method was evaluated and compared to the default method available in NONMEM VI via simulations and subsequent analysis of pharmacokinetic data using an intravenous linear one-compartment model with various distribution shapes for the parameters CL and V and with various numbers of observations per subjects. The extended-grid method was further evaluated via simulations using two more complex models with more realistic parameter estimates which were based on real case studies (sparse and small datasets), showing poor performance with the default method. Additionally, the method was also evaluated using a real data. Details of all studied cases follow.

Studied Scenarios

Simulation Exercise Based on a Simple Model

Simulation Settings

An MC simulation study was carried out in order to investigate estimation properties of the new methods. The major interest was to explore if this method can accurately and precisely estimate entire distribution shapes (specifically those other than the normal) and if it can provide an improvement compared to the default method available in NONMEM VI. Thus, two simulation conditions when the default nonparametric method is likely to perform poorly were chosen to be studied: small sample scenario (low number of subjects) resulting in low numbers of support points and sparse data scenario (low number of observation per subject and high residual variability) leading to uninformative data and poor estimation of support points.

Small Sample Scenario

Pharmacokinetic datasets (100 replicates) were simulated from a one-compartment intravenous bolus model with first-order elimination. Each dataset comprised of 30 individuals with two observations per individual placed at 1 and 5 h postdosing. The parameter distribution of V was simulated to be heavy-tailed. This was done by simulating a mixture of two normal distributions with the same mean, which was 100 L, but different variance, which was 0.09 (80% of the subjects) and 1 (20% of the subjects). The mean value of CL used for simulation was 10 L/h with variance of 0.09, giving rise to a half-life for a typical patient of 6.93 h (2.91 and 17.49 h for the fifth and 95th percentile).

A unit dose was assumed, and the homoscedastic residual error model was used, shown in Eq. 1, where yij stands for the simulated observations, equation M1 for individual prediction, and εij for a normally distributed random variable with a mean of 0 and default variance of 0.01, giving rise to residual variability of 10% (CV) (11).

equation M2

Sparse Data Scenario

This scenario was similar to the previous one, but each dataset included 200 subjects. Also, the residual variability used for simulation purposes was higher (30% CV corresponding to a variance of 0.09) in order to generate data with considerable stochastic noise, thus with a low information content per subject.

Analysis of Simulated Data

Each simulated dataset was analyzed with: (1) the default nonparametric estimation method and with (2) the extended-grid estimation method. In addition, each simulated data was analyzed with a true (mixture) model using FO and FOCE method to generate a reference. For small sample scenario, parameter estimates in parametric analysis were obtained using the FOCE method, while, for sparse data scenario, the FO method was used, although there was no major difference with respect to parameter estimates between FO and FOCE in this particular case. We also performed estimation with the default nonparametric estimation preceded with LAPLACE for a reference. When the extended grid method was used, the default nonparametric grid (EBEs) was extended with 200 additional support points generated via simulation using mean and variance obtained in preceding parametric analysis, for both scenarios. Thus, for the extended-grid method, the nonparametric distribution was estimated at 230 and 400 support points, for small sample and sparse data scenario, respectively. For the sparse data scenario, addition of 400, 800, and 1,600 support points in the nonparametric grid was also explored for possible improvements over the 200 originally used additional support points. Generation of additional support points via simulations using an inflated variance, compared to that estimated from the parametric step, was also investigated for the extended-grid method. This was studied in cases when 400, 800, and 1,600 additional support points were generated to enrich the nonparametric grid. Variances were inflated 1.5 and two times the original value obtained in preceding parametric analysis. In this case, change in the nonparametric objective function value (NPOFV) was also observed for models with the same number of support points but different variance inflation factor.

Assessment of the Results

Cumulative density functions (cdf) of estimated distributions were compared visually with the cdf of the true (simulation) distribution, for randomly chosen samples. Also, prediction intervals based in 100 estimated nonparametric distributions were constructed and compared with the true distribution. Numerically, estimated distributions were evaluated at different percentiles (5th, 10th, 20th, 25th, 30th, 50th, 70th, 75th, 80th, 90th, and 95th). Parameter values at these percentiles for the true distribution were obtained through large sample simulations (n = 5,000). Comparison between the estimated and true distributions of CL and V were made at the evaluated percentiles using the relative percentile error (RPE) according to Eq. 2 where Pest stands for the estimated parameter value and Ptrue for the true parameter value at the given percentile.

equation M3

The bias and precisions in parameter estimates at given distribution percentiles were presented by plotting the distribution of RPE as box plots, where the median of the RPEs at each percentile characterizes the bias of the method (Eq. 3) while the precision of the method is given by the variance of the sample of RPE at evaluated percentiles.

equation M4

Simulation Exercise Based on the Design and Model from Real Case Studies

Sparse Pediatric Example

The dataset described previously by a parametric NONMEM model was used in this study to further evaluate the extended-grid method. This was a study of desmopressin pharmacokinetics conducted in 84 children with nocturnal enuresis aged 6 to 12 years (12). It was a dose-escalation study where six doses of desmopressin ranging from 30 to 480 µg as well as placebo were administered sublingually. Each child was sampled twice, at one early (0–2 h postdosing) and one late (3–24 h postdosing) time point, and one third of the children had an additional sample at an intermediate time point. In total, 139 samples were available from 72 children. A one-compartment model with first-order absorption and transit compartment chain for description of the absorption delay described the data best. This model contained random effects in five parameters (CL, V, Ka, mean transit time (MTT), and number of transit compartments (NN)). Details of this study are published previously (12). The final parameter estimates of population means and variances from this (parametric) model were considered as the true parameters. One hundred datasets were simulated from this model using original sparse study design in children.

Parallel Bioavailability Example

The data comprised of 17 individuals, of whom nine received the drug via intravenous administration route while the remaining eight subjects received the orally administrated drug. In total, 122 observations were available for the analysis. The data were kindly provided by Dr. Rik Schoemaker. The pharmacokinetics of the drug was successfully described by a two-compartment model with first-order absorption and a lag time. Random effects were introduced in five parameters, elimination rate constant (k), distribution rate constant (k23), absorption rate constant (ka), bioavailability (F), and lag time parameter (alag). The final parameter estimates were considered as the true parameters (Table I) and, using this model and the study design, 100 datasets were simulated.

Table I
Parameter Values Used for Simulations of the Real Data Example 2

Analysis of the Datasets

For both examples, nonparametric distributions of all five parameters were estimated using the default nonparametric method in NONMEM as well as the extended-grid method under three different conditions: using 200, 500, and 1,000 additional support points to enrich the nonparametric grid. Results were evaluated in a similar manner as for the simulation study described above. Additionally, for each parameter distribution, 90% prediction intervals based on 100 nonparametric distributions estimated on the extended grid were constructed and compared both to the true distribution and to the nonparametric distribution from the original real data.

Evaluation of the Method with the Real Data

The dataset described previously by a parametric NONMEM model was used in this study to further evaluate the extended-grid method. This was a study of desmopressin pharmacokinetics conducted in 28 healthy volunteers (12). A single dose of 240 μg was studied, and rich sampling was conducted (16 samples per subject). A one-compartment model with first-order absorption and transit compartment chain for description of the absorption delay described the data best. This model contained random effects in five parameters (CL, V, Ka, MTT, and NN). Details of this study are published previously (12).

First, EBE distribution and nonparametric distribution using default nonparametric methods were estimated using a data-rich study design. Shrinkage was assessed for each estimated distribution. Then, data were censored at two different levels, inducing low and high shrinkage. With these censored data, estimation of nonparametric distribution was repeated using the default nonparametric method and extended-grid method.

Visual assessment of results was done for two parameters (CL and V) which the original distribution estimated using full data showing negligible shrinkage (<1%) and therefore were considered as reliable and close to the true values.


Simulation Exercise Based on a Simple Model

When visually inspected by the comparison of their cdf, the parameter distributions estimated using the extended-grid method showed good agreement with the true parameter distributions, for both small sample scenarios and sparse data scenario as shown in Fig. 1 where a randomly chosen example from both scenarios is shown.

Fig. 1
A visual comparison of distributions for volume of distribution estimated using a default nonparametric method and extended grid method with the true distribution in small and sparse dataset illustrations. A randomly chosen example is shown for both scenarios ...

When evaluated numerically, in the small sample example, bias was slightly larger for the extended-grid method but still low for both methods (less than 10%); however, imprecision was decreased as illustrated in Table II. The average value of RPE variance at each evaluated percentile was 586 for the default method (range 89–4,296) and 344 (range 52–2,453) for the extended-grid method (Table II and Fig. 2). This can also be seen in Fig. 3 where prediction intervals created based on 100 estimated distributions were somewhat tighter for the extended-grid method than for the default method, which suggests improved precision of the new method.

Table II
Performance of the Default Nonparametric Estimation Method and the Extended Grid Method with 200 Additional Points of Support for the Small Datasets Evaluated as Bias and Precision at Different Parameter Distribution Percentiles
Fig. 2
RPEs (n = 100) at 11 evaluated distribution percentiles (5th, 10th, 20th, 25th, 30th, 50th, 70th, 75th, 80th, 90th, and 95th) using default nonparametric estimation method in NONMEM (default), extended grid method (extend), and the true ...
Fig. 3
Prediction intervals created based on 100 estimated nonparametric distributions for the small sample scenario. The true distribution is shown with the thick line

For the sparse data illustration, bias was slightly lower for the extended-grid method compared to the default method, 3.3% (range 1.3–8.4%) versus 5.9% (range 1.3–9.9%), and imprecision was substantially decreased from the average variance across all percentiles of 204 (range 44–563) for the default method down to 80 (range 20–238) for the extended-grid method (Fig. 4). Addition of 400, 800, and 1,600 support points in the extended grid improved results even further with decreased bias down to 2.4% (range 0.0–6.6%), 2.4% (range 0.0–5.3%), and 2.9% (range 0.0–5.9%), respectively. Imprecision was similar when 200 additional support points were added: 87 (range 20–337), 85 (range 15–337), and 83 (range 16–319) for 400, 800, and 1,600 support points, respectively. The summary for bias and imprecision at different percentiles for each evaluated case is shown in Table III. Using inflated variances to generate additional points of support did not further improve results for any of the studied cases. This was confirmed both by increase of biases (Table IV) as well as by the increase in NPOFV (Fig. 5). A slight decrease of imprecision was observed at the tail distribution percentiles.

Fig. 4
RPEs (n = 100) at 11 evaluated distribution percentiles (5th, 10th, 20th, 25th, 30th, 50th, 70th, 75th, 80th, 90th, and 95th) using default nonparametric estimation method in NONMEM (default), extended grid method (extend), and the true ...
Table III
Performance of the Default Nonparametric Estimation Method Preceded by FO and LAPLACE and the Extended Grid Method with 200, 400, 800, and 1,600 Additional Points of Support for the Sparse Datasets Evaluated as Bias and Precision at Different Parameter ...
Table IV
Influence of Different Inflation Factors (1.5 and 2) Used to Determine the Distribution Variance used to Generate Additional (400, 800, and 1,600) Points of Support for Sparse Data Scenario
Fig. 5
Change in nonparametric OFV between models with default and inflated variance was used to generate the support points. Three cases are shown: generation of 400, 800, and 1,600 support points

Simulation Exercise Based on the Design and Model from Real Case Studies

In the sparse pediatric example, the default method performed poorly, especially for the estimation of the absorption parameter distributions. The average absolute bias across all evaluated percentiles and all parameter distribution was 34% (range 0.2–172%), and imprecision was 1,548 (9–38,592). The extended-grid method improved results substantially for all parameter distributions both in terms of decreased bias with average value of absolute bias across all evaluated percentiles and all parameter distributions of 3.6% (range 0.1–11.8%) and imprecision value of 109 (range 9–1,305). Also, prediction intervals derived based on the distribution estimated on the extended grid showed a good agreement with the true parameter distribution. This was not the case with the default nonparametric distribution estimated using the true model and the original dataset (Fig. 6).

Fig. 6
Sparse pediatric example: 90% prediction intervals (blue line) and a median (black line) for each of the five parameter distributions created based on 100 estimated nonparametric distributions using extended grid method (1,000 additional support points) ...

In the parallel bioavailability example, the default method also showed poor behavior with estimated distributions substantially deviating from the true distribution (Fig. 7). An interesting trend was observed in absorption parameter (ka, F, and alag) default nonparametric distributions—the majority of the population probability was associated with the typical value. Conversely, the distributions estimated on the extended grid were in a good agreement with the true distribution (Fig. 7 and Table V).

Fig. 7
Parallel bioavailability example: 90% prediction intervals (blue line) and a median (black line) for each of five parameter distributions created based on 100 estimated nonparametric distributions using extended grid method (1,000 additional support points) ...
Table V
A Real Case Study—Parallel Bioavailability Example

Evaluation of the Method with the Real Data

In this example, nonparametric distribution estimated using extended grid offered an improvement compared to the default nonparametric method examined as visual agreement with the empirical Bayes distribution and nonparametric distribution estimated using a default method and rich design. This was the case for both CL and V and also for both levels of censoring, inducing low and high shrinkage (Fig. 8).

Fig. 8
Example with the real dataset. Comparison of the nonshrunk empirical Bayes distribution (green line) and nonparametric distributions estimated using (1) default NONMEM method and full rich design (red line), default NONMEM method and sparse design (blue ...


A new method for the estimation of the nonparametric distribution has been developed and evaluated in this work. One of the aims with this work was to develop a method with improved estimation properties for small sample sizes. This is also a challenging problem for parametric methods: due to small sample size, it is hard to estimate variance precisely and to assess the distributional assumption. In general, it is preferable for any method if the interindividual variability model is based on a large number of subjects. The default behavior of existing nonparametric methods is that the parameter distributions is defined at the maximal number of support points equal to the number of individuals in the dataset. This is often seen as a disadvantage of the nonparametric methods in general, when the estimated distribution is to be used for simulation of new future subjects—the parameter values that can be sampled are limited. Also, tails of the distribution are usually not well defined, which is the case with few subjects under all circumstances. Therefore, from nonparametric distribution, new subjects are sometimes simulated using a smoothed distribution. Such smoothing is generally arbitrary, and the final smoothed distribution may contain parameters that have no support in the studied population. The extended-grid parameter distribution will have support in real data for all parameters in the distribution. With the extended-grid method, the nonparametric distribution is estimated at a larger number of support points which also allows the uncertainty in the parameter distribution to be incorporated. This uncertainty relates to uncertainty in individual parameter estimates but not the uncertainty coming from sample size (number of individuals in the dataset).

The extended-grid method estimates parameter distributions without increasing bias in the parameter estimates and with improved precision over the default method. Precise estimation of parameter variability is challenging when individual data are sparse and noisy, not only for nonparametric methods but also for less numerically demanding parametric methods. An additional challenge for nonparametric methods is to precisely estimate the entire probability density function from such data. The default nonparametric method in NONMEM, similar to the parametric method, performs not as well when this type of data is analyzed. This is because the nonparametric estimation method in NONMEM uses EBEs as support points, and EBEs are known to be unreliable when individual data are uninformative (13). In such cases, the range of support points does not include the full parameter range, and the subsequent nonparametric estimations result in biased and imprecise parameter estimates. Thus, another important aim with development of this new method was to improve estimation properties when individual data are uninformative. Indeed, the approach of obtaining support points implemented in extended-grid method seems to be a better alternative than using only EBEs themselves. The extended-grid method, as implemented here, keeps the EBEs in its grid but also fills the parameter space with additional support points.

In this work, as a default method, we added 200 additional points of support in the nonparametric grid for both studied cases. This number seemed to be sufficient in order to obtain satisfactory estimation properties of the method for these particular cases. We also tested addition of 400, 800, and 1,600 support points which improved results even further. Also, the estimation procedure itself was not substantially longer with respect to time. Run times for nonparametric estimation on the grid containing 200, 400, and 800 additional support points were 12, 15, and 22 s, respectively, for rather simple cases which we tested (PK models with two to five parameters). However, with each new parameter to be estimated, the parameter space that need to be filled up with support points will have an extra dimension, which will also require addition of new support points. Thus, theoretically, if one would like to keep the same density of the nonparametric grid, increase in number of support points to be added would follow the power function with the power equal to the number of dimensions. If a modeler would like to preserve a fine grid with a large number of steps (support points) for a complex model, it is clear that number of support points needed would increase rapidly (number of support points = number of steps(parameters)). In our example, 200 additional support points were sufficient to fill two-dimensional space and to produce satisfactory results. If the parametric method is expected to have performed well, a “low” number is expected to do well. With low number, at least as many support points as in the original study are considered, as this number usually indicates how well the parameter distributions is aimed to be characterized in the population. If the parametric fit is shrunk, more support points are needed. Another criterion for choosing number of additional support points can be obtained by comparison of parametric and nonparametric estimated parameter distributions. If these are similar between two methods, it seems sufficient to simulate low numbers of additional support points. If the nonparametric distribution differs substantially from the parametric in the presence of shrinkage, this is an indication that more support points need to be generated from the parametric fit. If nonparametric variances are higher than the parametric ones, generation of support points by simulations from inflated parametric variances can be considered.

In order to generate additional support points, we used simulation from the normal distributions using a mean and variance estimated in preceding parametric analysis. Again, for the studied conditions, this seemed to be satisfactory. The majority of parametric algorithms including FOCE or stochastic algorithms (SAEM, MC-PAM) are rather good at quantifying the variance of parameter distribution (14). However, this method is disputable if the parametric estimate is expected to be biased towards too low variances or correlations too close to ±1. This was also observed sporadically in our analysis. If that is the case, variances or covariances to simulate support points need to be changed. In our case, when parametric parameter variances and covariances were reasonably adequate, inflation of variances did not help to improve results further. Also, the NPOFV, which is the measurement of the exact likelihood, increased. There is little experience of using NPOFV for model building decisions. Rival nonparametric models are seldom nested so the likelihood ratio test, based on objective function value (OFV) cannot be used. This is because, with any new addition of a structural model parameter, the number of additional parameters to be estimated in the model will be equal to or less than the number of subjects. Moreover, implementation of the likelihood ratio test procedure for inference and model selection is based on the assumption that the intraindividual and interindividual errors are normally distributed, which is not the case within nonparametric framework (1). For nonnested models with different number of parameters, where each support point represents a parameter, criteria such as the Akaike, Schwarz, or Bayesian information may possibly be useful (1), but limited published information is available on the use of such criteria. However, for two models which have the same number of parameters, a comparison based on OFV seems appropriate with a lower OFV indicating a better fit. Generation of support points using simulations from other distribution than normal (e.g., uniform, using parameter estimates, mean, and variance, from preceding parametric fit) was also tested, and this change did not affect analysis to a greater extent (data not shown).

One of the characteristics of majority of nonparametric methods is that these methods are not able to estimate simultaneously interindividual and residual variability. An exception to this is NPAG implemented in USC*PACK (3,4), for which applications (15,16) have reported estimated residual error. The underlying methodology for this estimation has to our knowledge not been published. Therefore, it is a common procedure to assess the residual variability a priori via parametric analysis, and subsequently the nonparametric estimation is conditioned on this estimate. Therefore, it is essential to use both appropriate models for residual variability as well as estimation methods which are able to estimate components of residual variability well. It has been shown that, in some circumstances, FO is not able to separate variability components well, and all variability is expressed as a residual variability, which further affects the subsequent nonparametric estimation to a great extent. The extended-grid method also suffers from this phenomenon, which was noticed in our analysis (data not shown). Therefore, more sophisticated algorithm, like FOCE or stochastic algorithms, such as SAEM, MC-PAM, and MCMC would be appropriate to employ for separation of variability components.

The method was also tested using simulations based on two real case studies. In the first pediatric example, data available to develop a model were rather sparse. Indeed η-shrinkage in the final model ranged from 20% to 70%, indicating rather unreliable individual parameter estimates. This will have negative consequences for generation of support points in the default nonparametric method, but the extended-grid method showed rather good behavior even in this case of sparse data. The second real data example was based on the small study and it was motivated by discussion on NONMEM user list from July 4th 2007 ( where certain problems with the default version of the method were pointed out. Firstly, the default method in NONMEM would give a different nonparametric outcome depending on the method used to run the parametric step (FOCE or FOCE with INTERACTION) regardless the fact that these methods produce exactly the same parametric fit, which was treated as a software bug. Another more worrying problem is that the default method assigns parameter values also to subjects that do not require that particular parameter in their model, for example, absorption parameters are assigned to subjects receiving an intravenous drug, where truly, for these subjects, absorption parameters should be treated as missing. This is indeed a consequence of EBE estimation which, in the case of no information on particular parameters from an individual, assigns the most probable value to that individual, which is the population mean. This causes a problem for subsequent nonparametric estimation, which treats these support points (ETA = 0) as active support points and, as a consequence of lack of information, a large part of nonparametric probability is associated with the population value. This affects the entire joint density estimation, which results in a distorted distribution shape, not only for the absorption parameters but also for the other parameters common to all subjects.

The extended-grid method offered a solution to both problems. It gave the same result no matter which method was used to estimate EBEs (FOCE or FOCE with INTERACTION) and provided estimates of distribution shapes which were in a good agreement with the true ones, despite the small dataset and the parallel study design with some subjects not having information about some parameters.

Another method for generation of support points has been also reported earlier—the inflation method (17). With respect to EBE shrinkage, this method provides a better range of support points by estimating EBEs from the parametric distribution with inflated variances. We also compared this method with the extended-grid method (data not shown), and it showed adequate performance when variances were inflated with twice the value of the estimated standard deviations; however, as it still uses EBEs as support points, this method did not perform that well when data were even sparser. Indeed, the inflation method was never better than the extended-grid method. However, it is easier to implement as it only requires users to change the variance values prior to EBE estimation. The inflation method may be useful for large datasets with not too large EBE shrinkage.

The evaluation of method was also exemplified using a single real data example. The real data studies often present problems not encountered in simulation studies. The main reason is that the model assumptions are often wrong for real data. Also, the true distribution is never known. In our example, we have used nonshrunk EBE distribution and the default nonparametric distribution estimated using rich design as a reference distribution because these were not exhibiting shrinkage and therefore were considered reliable. However, showing benefit of one approach over another on a single dataset is seldom convincing, and therefore evaluation with the real data would require more extensive evaluation with the larger number of the real datasets/models in order to overcome such limitations for using real data. However, such an extensive investigation definitely requires a separate presentation and is beyond the scope of present manuscript.

The extended-grid method is more complex to implement and practically requires several steps which are not that straightforward. Because of that, we have implemented an automated version of the extended-grid method in the software PsN (freely available at (10). In this automated version, the user has an option of choosing the number of support points to generate, as well as to decide on the inflation factor to use for simulations.

The availability of nonparametric estimation method in NONMEM enabled a wider use of nonparametric estimation methods in general. In addition to recommended use of these methods whenever normality assumption of random effects does not hold, these methods are also recommended to be used as an aid in parametric model building, as a replacement of the final parametric model and as a bridging tool between NONMEM and other nonparametric software. However, the relative merit of different nonparametric methods possible to implement in NONMEM VI depends on data and model structure: for data with many subjects and rich individual data, the new methods offer no great advantage, as the default method performed already well in these scenarios in earlier studies (7,17). For sparse individual data and for small datasets, the extended-grid method was preferable. The inflation method was never better than the extended-grid method but easier to implement and may be a viable alternative for large datasets with sparse individual data.


In conclusion, a method for nonparametric estimation has been developed and explored which resulted in improved estimation properties compared to the default method in NONMEM for small and/or sparse data. The extended-grid method is automated, using the program PsN, for implementation into the NONMEM program.


We thank Dr. Rik Schoemaker for providing an example for parallel bioavailability example, Pontus Pihlgren for technical assistance, and Dr. Shasha Jumbe for valuable comments on the manuscript. We would also like to thank two anonymous reviewers for their valuable comments which greatly helped to improve the manuscript.


1. Davidian M, Giltinan DM. Nonlinear models for repeated measurement data. Boca Raton: Chapman & Hall/CRC; 1995. pp. 191–217.
2. Ette EI, Wiliams PJ. Pharmacometrics: the science of quantitative pharmacology. Hoboken: Wiley; 2007. pp. 223–245.
3. Bustad A, et al. Parametric and nonparametric population methods: their comparative performance in analysing a clinical dataset and two Monte Carlo simulation studies. Clin Pharmacokinet. 2006;45(4):365–383. doi: 10.2165/00003088-200645040-00003. [PubMed] [Cross Ref]
4. Leary R, Jelliffe R, Schumitzky A, Guilder MV. An adaptive grid non-parametric approach to pharmacokinetic and dynamic (PK/PD) population models. In 14th IEEE Symposium on Computer-Based Medical Systems (CMBS'01) p. 0389; 2001.
5. Mallet A, et al. Nonparametric maximum likelihood estimation for population pharmacokinetics, with application to cyclosporine. J Pharmacokinet Biopharm. 1988;16(3):311–327. doi: 10.1007/BF01062140. [PubMed] [Cross Ref]
6. Savic RM, Kjellsson MC, Karlsson MO. Evaluation of the nonparametric estimation method in NONMEM VI. Eur J Pharm Sci. 2009;37(1):27–35. doi: 10.1016/j.ejps.2008.12.014. [PubMed] [Cross Ref]
7. Savic RM, Karlsson MO. Sequential parametric-nonparametric modeling - advantages, disadvantages and usefulness. In ACoP, Arizona, US: Tucson,; 2007.
8. Leary R. An evolutionary nonparametric NLME algorithm. in PAGE. 2007. Copenhagen, PAGE 16 (2007) Abstr 1100 [].
9. Karlsson MO, Savic RM. Diagnosing model diagnostics. Clin Pharmacol Ther. 2007;82(1):17–20. doi: 10.1038/sj.clpt.6100241. [PubMed] [Cross Ref]
10. Lindbom L, Pihlgren P, Jonsson N. PsN-Toolkit—a collection of computer intensive statistical methods for non-linear mixed effect modeling using NONMEM. Comput Methods Programs Biomed. 2005;79(3):241–257. doi: 10.1016/j.cmpb.2005.04.005. [PubMed] [Cross Ref]
11. Carroll RJ, Ruppert D. Transformation and weighting in regression. London: Chapman & Hall; 1988.
12. Osterberg O, et al. Pharmacokinetics of desmopressin administrated as an oral lyophilisate dosage form in children with primary nocturnal enuresis and healthy adults. J Clin Pharmacol. 2006;46(10):1204–1211. doi: 10.1177/0091270006291838. [PubMed] [Cross Ref]
13. Savic, RM, Karlsson MO. Importance of shrinkage in empirical Bayes estimates for diagnostics: problems and solutions. AAPS J2009, Aug 1st (in Press). [PMC free article] [PubMed]
14. Girard P, Mentre F. A comparison of estimation methods in nonlinear mixed effects models using a blind analysis. In PAGE 14 (2005) Abstr 834. Pamplona, Spain,; 2005.
15. de Wildt SN, et al. Population pharmacokinetics and metabolism of midazolam in pediatric intensive care patients. Crit Care Med. 2003;31(7):1952–1958. doi: 10.1097/01.ccm.0000084806.15352.da. [PubMed] [Cross Ref]
16. Drusano GL, et al. Use of preclinical data for selection of a phase II/III dose for evernimicin and identification of a preclinical MIC breakpoint. Antimicrob Agents Chemother. 2001;45(1):13–22. doi: 10.1128/AAC.45.1.13-22.2001. [PMC free article] [PubMed] [Cross Ref]
17. Savic RM, Kjellsson MC, Karlsson MO. Evaluation of the nonparametric estimation method in NONMEM VI beta. PAGE 15 (2006) Abstr 937,; 2006.

Articles from The AAPS Journal are provided here courtesy of American Association of Pharmaceutical Scientists