NONMEM® (NONlinear Mixed Effects Model) is the most widely used software for the population pharmacokinetic/pharmacodynamic (PK/PD) modeling and simulation, also is the standard tool for drug development as well as regulatory approval. NONMEM® uses inter-subject variability parameter η to describe the individual differences in PK/PD parameters. However the standard error (SE) of the empirical Bayes estimate (EBE) of the inter-subject variability parameter η is not provided in the NONMEM® output file. The SE of EBE is an enhanced model diagnosis tool to assess the precision of the estimated inter-subject variability parameter η. Also the SE of EBE can be used for simulation considering the uncertainty in the inter-subject variability [1].

The SE of EBE of the inter-subject variability parameter η is not available in the output file of NONMEM®. Therefore postprocessing after a modeling NONMEM run is necessary to obtain the SE of EBE. Writing one's own code requires careful thinking and programming expertise. A Perl based automated utility to calculate the SE of EBE is available in PsN® version 2.2.6 or later [2]. However, a more direct method to obtain the SE of EBE is possible with NONMEM® VI using one of the internal variables, POSTV. Unfortunately, this feature is not available in the previous versions of NONMEM®. The POSTV matrix is not available in the new version of NONMEM® 7, either. However, the elements of the POSTV matrix are available as ETC columns in the .phi file, one of the new output files from NONMEM® 7.

The objectives of this study are to illustrate how to directly obtain the SE of EBE by the use of NONMEM® VI internal matrix POSTV and to evaluate the performance by comparing the results with other methods. The contents of this paper are as follows. In the methods section, theoretical backgrounds for calculation of the EBE and the standard errors of EBE are overviewed along with the descriptions on the four models and the two datasets used to assess the different ways of calculating the standard of errors of EBE. Then detailed explanations on the three methods, including the suggested simple method of obtaining the covariance matrix of EBE in NONMEM® VI using verbatim code, are provided. In the results section, the numerical results are summarized and compared. Lastly, the discussion section concludes with a few remarks and conclusions.

The following notations are used for PK/PD modeling using NONMEM®:

θ, Ω, and Σ are constants while η and ε are random variables. The function f1 represents a structural model describing the relationship between the PK/PD observations and the model parameter θ, while η and ε represent the stochastic model components describing the randomness unexplained by the structural model.

The typical PK/PD of the population is summarized by θ, and individual PK/PD parameters are expressed as a combination of θ and η. For example, systemic clearance is expressed as

CL=θ_CL·exp(η_CL), Equation 1

where θ_CL is the typical value and η_CL is the corresponding inter-subject variability following normal distribution with mean zero and covariance, e.g., Ω₁₁. In other words, the typical value of a PK/PD parameter is the value of PK/PD parameter with η=0. F is called the typical predication (PRED) when η=0, and the individual prediction (IPRE or IPRED) when η=η̂ the EBE of individual eta.

In addition to the inter- and intra-subject variability, there can be another level of random variability called inter-occasion variability which describes the random variability among the different treatment occasions, for example, different visits. Implementation of inter-subject and intrasubject variability can be accomplished using the corresponding NONMEM® variance parameters η and ε. However, implementation of inter-occasion variability requires an elaborate coding using the inter-subject variability parameter η, because η now has to match the new variability introduced to model the inter-occasion differences. The developed method of computing standard error (SE) of EBE using NONMEM® VI internal matrix POSTV can be applied to η for inter-subject variability as well as η for inter-occasion variability without further modification of NONMEM® control code.

The random variable η cannot be directly observed but can be estimated using subject level observations (Y, i.e., dependent variable) and the final estimates of θ, Ω, and Σ. The estimates of each subject's η (η̂) has several namesrealized η, post hoc η, empirical Bayes estimate (EBE) of η, or maximum a posteriori (MAP) estimate of η. There are several numerical methods to estimate the aforementioned PK/PD parameters. The available estimation methods from NONMEM® VI and their characteristics are as follows [3]. FO method employs "First Order" linear approximation of F with respect to η. Laplacian method (L) uses "Laplacian" integral approximation of the objective function using up to the second order partial derivatives. FOCE method stands for "First Order Conditional Estimation," for which individual η (η̂) is estimated during the minimization process [3]. Individual η estimates (η̂) are not calculated during the minimization process for the FO method, however they can be estimated with the final θ, Ω, and Σ with the POSTHOC option. In the FOCE or Laplacian method, computation of the individual η estimates (η̂) is necessary during the minimization process. The difference of FOCE objective function with respect to the Laplacian objective function is the use of the first order partial derivatives for the linear approximation whereas the Laplacian method uses up to the second order partial derivatives [4]. The computation time increases in the order of FO, FOCE, and Laplacian. The INTERACTION estimation option considers interaction of η and ε, and uses η̂ instead of 0 for η during the calculation of variance of Y. INTERACTION option does not affect the output for the additive residual error model, but produces different estimates which are regarded as more accurate for other residual error models (proportional error, combined additive and proportional error, power function error, etc.).

When f₁ is linear with respect to η, joint objective function is used for estimation of θ, Ω, Σ, and [5]. When f₁ is nonlinear with respect to η, NONMEM constructs a separate objective function to estimate EBE in addition to the objective function to estimate θ, Ω, and Σ. Objective function for EBE in NONMEM® as well as ADAPT 5 [6], a similar PK/PD modeling and simulation software, is as follows:

where ij denotes j-th observation of the i-th individual, and LL() is the log-likelihood.

The uncertainty of the estimated θ, Ω, Σ, and η̂ can be obtained by computing the corresponding standard error which is the standard deviation of estimates. However, special consideration is required for η̂ because there is an upper limit for the uncertainty in η̂. Even if there are no data for a particular individual, ω̂ (estimated inter-subject variability for the population data at hand) can be an estimate of the uncertainty about η̂ of this individual, because of the assumption that the individual belongs to the population characterized by ω̂. For this case, individual η shrinkage, the standard error of divided by the corresponding ω̂, could be a better measure of uncertainty [7].

In NONMEM® VI, the standard error of EBE of η̂ can be computed using the following asymptotic expression for covariance of an estimate, assuming that the estimated θ, Ω, and Σ are the true values [5].

The usefulness of EBE of η is manifold. Individual PK/PD parameters and accordingly individual drug concentrations or effects are described using EBE [8-10]. EBE can be used for screening covariates for the structural model development [11]. If relationships between covariates and EBEs do not exist, EBE theoretically should have no trend with the covariates. Apparent relationship between EBE and covariate indicates the necessity of further refining the structural model, for example, including that covariate in the structural model. However, such a method is not recommended when there is difference between the empirical distribution of η̂ and the estimated Ω, i.e., when the shrinkage is greater than about 20%. The empirical distribution of EBE can be used to test the normal distribution assumption on Ω, or the validity of the mixture model, or the necessity of using some other kind of nonparametric distribution. Since EBEs are used in the objective function calculation in FOCE or Laplacian method, both methods rely on accurate EBEs. In NONMEM® VI, EBE can be used as support points to calculate the joint density for NONPAR (nonparametric) estimation option [12].

Five estimation methods were tested for each model: FO (First Order), FOCE (First Order Conditional Estimation), FOCEI (First Order Conditional Estimation with Interaction), L (Laplacian), and LI (Laplacian with Interaction). FOI (First Order with Interaction) was explored using Methods 1 & 2, however the result was not used for comparison because POSTV was not available with METHOD= ZERO INTERACTION option. Other estimation options were set at default values except that the maximum number of function evaluation (MAXEVALS) was set at a generous number of 9999. As a measure of comparing the results from the three methods, the frequency of discrepancy in the results greater than 5% compared with the R code result (Method 1) is shown in Table 3. Method 1 is used as reference to compare the degrees of discrepancy of computation results from the other methods, because the R code was straightforwardly implemented equation 3. On the other hand, the detailed computation algorithms of the other two methods are either unavailable to the public or very difficult to track. EBEs and standard errors computed by method 1 are in appendix tables.

There was no discernable difference in the point estimates of EBE of η among all of the three methods, for the various estimation methods and the datasets. The degree of differences in the results from the three computation methods for the SE of EBE of η depended on the employed estimation method. All the three methods produced almost the same results when Laplacian estimation method (L or LI) was used. For the conditional estimation methods (FOCE or FOCEI), there was no noticeable difference between the computed SE estimates from R code (Method 1) and PsN (Method 2). When FOCEI estimation method was used for the sparse sampling dataset of phenobarbital, the percentage occurrence of POSTV method results (Method 3) differing more than 5% compared with R code result (Method 1) was less than 1%. For the extensive sampling dataset of theophylline, greater than 5% difference in EBEs occurs 3% to 17% with FOCEI. The degree of difference got larger for the simple FO estimation. Greater than 5% difference between Method 1 and Method 3 occurred more frequently (21% to 86%) with the FO method.

Estimation results are usually provided as predicted mean or typical values for the unknown parameters of the system along with measures of inter-correlation and precision associated with the parameter estimates. Standard errors of the parameters are easy to explain and understand about the precision of the parameters. In population analysis where hierarchical structure (ie, population level and individual level) is involved, investigation of the SE of EBE for individual level can be used as an enhanced model diagnosis tool to assess the precision of the estimated inter-subject variability parameter as well as for simulation considering the uncertainty in the inter-subject variability, in addition to the SE of population parameters.With the availability of SE of EBE during simulation, the discussion about the central tendency or outlying values for individual predicted values can be possible, which can be a useful information for developing indivualized therapy.

Three methods of computing the standard error of empirical Bayes estimate of the inter-subject variability parameter η for NONMEM were explored and the performance was evaluated. The SE of EBE of the inter-subject variability parameter η can be computed in a direct and robust way by using the POSTV matrix in NONMEM® VI without additional postprocessing efforts. The POSTV matrix is not available in other versions of NONMEM®. However the elements of the POSTV matrix are available as ETC columns in the .phi file from the new version of NONMEM® 7. Such SE of EBE from NONMEM® 7 was almost identical with the POSTV from NONMEM® VI. The proposed simple and robust method of computing SE of EBE using POSTV is not reported in the standard NONMEM® output file. The reason might be that there are other methods for computing SE of EBE using uncertainties associated with θ, Ω and Σ together. The computation differences with respect to the result from R code may be due to the fact that R function "optim" uses numerical calculation of gradient and hessian while NONMEM® uses analytical partial derivatives partially for the PK models using ADVAN2 routine. However, when LAPLACIAN is used, NONMEM® also uses numerical methods for gradient and hessian calculation so that the result showed no difference with respect to PsN and R.

There are a few limitations of using POSTV for obtaining SE of EBE. The method assumes no uncertainties for the estimated θ, Ω and Σ, and cannot be used with the METHOD=ZERO INTER estimation option. The discrepancies of the computation results with respect to the other two computation methods were observed, and the differences might be attributed to the analytical versus numerical calculation of derivatives.

PsN is a convenient tool to obtain SE of EBE, however PsN fits modified model to the modified individual datasets. In addition, PsN may produce quite different point estimates of EBE when the EBE values are close to zero. Therefore, computing the point estimate of EBE using individual data should be used with caution. The source of this difference might be due to the different minimization algorithm for EBE computation to the algorithm for the estimation of θ, Ω and Σ. Not much is known about what algorithm is used for computing EBE, whereas quasi-Newton type derivative free algorithm (a variant or modification of derivative free Davidon-Fletcher-Powell algorithm in ZXMIN1. FOR) is used for estimating θ, Ω and Σ.

In conclusion, POSTV matrix can be used for a simple and robust computation of SE of individual ETA values in NONMEM® VI. The results were comparable with other computation methods. The degree of discrepancy was smallest in LI method, smaller in FOCEI than FO.