117  NONMEM

Introduction to NONMEM Version VI

1 I.2.1 Estimation

1. With the Hybrid Estimation method and an intraindividual mean-variance model (see NONMEM Users Guide VII), the η, ε-interaction is now allowed. This interaction applies to the elements of η for which conditional estimates are obtained. Simply include the INTERACTION option on the $ESTIMATION record.

(See the $ESTIMATION Help Item.)

2. With the LAPLACIAN Estimation method and a mean-variance model (see NONMEM Users Guide VII), the η, ε-interaction is now allowed. In this case the numerical method for obtaining second partial η-derivatives must be used, and therefore too, a “slow” minimization method will be used. Simply include the INTERACTION option on the $ESTIMATION record.

(See the Help Item.)

3. Background: With the First-Order Estimation method, if the model is not already linear in ε, it is first approximated by a first-order Taylor expansion in ε about ε = 0. This approximate model is then approximated further by a first-order Taylor expansion in η about η = 0, followed by the η, ε-interaction term being dropped. Finally, the Extended Least Squares Estimation method (Beal 1984) is used.

Now the interaction term may be retained, resulting in the First-Order Estimation method with η, ε-interaction. Simply include the INTERACTION option on the $ESTIMATION record. The result is generally not very different from that obtained when the interaction term is dropped. Generally, a slightly less biased estimate of Σ is obtained, and when there are systematic interindividual differences in the magnitude of interindividual random variability itself, a slightly less biased estimate of Ω is obtained.

(See the $ESTIMATION Help Item.)

1.1 Stieltjes Method is fully implemented

A new population estimation method - called the Stieltjes Estimation method - is available in NONMEM VI 2.0. It is a two-stage method. A population conditional estimation method should be used to “quickly” obtain fairly reliable parameter estimates; this comprises the first stage. Include the METHOD=CONDITIONAL option on the $ESTIMATION record. The first-stage estimates are used as initial estimates in the second stage, wherein another objective function is minimized to obtain a final set of parameter estimates. The second stage takes much more time than the first stage, and the resulting parameter estimates are at least as good as the first-stage estimates, perhaps better. In the second stage, conditional estimates of η are not used, but the functions that are minimized to obtain conditional estimates are still evaluated at many values of η. Generally though, sufficiently good estimates are obtained by using conditional estimation methods, or even the First-Order Estimation method, and one needs to use the Stieltjes Estimation method only with somewhat pathological data designs. Simply also include the STIELTJES option (along with an optional GRID option) on the $ESTIMATION record.

1.2 4a. Further notes on Stieltjes for NONMEM VI 2.0

A set of tentative population estimates are first obtained using some 1st- or 2nd-order method specified on the $ESTIM record in the usual way. (However, the Hybrid method cannot be used). At this first stage, for each individual, a tentative value for the integral (i.e. an area) is obtained. Then numerical integration is used to obtain second-stage estimates. The integral will be exact in the case where the intraindividual model is of mean-variance kind and is linear in eta. A default grid on eta-space is used for numerical integration.

This is a description of how one might set one’s own grid. A grid is specified in terms of ellipsoidal coordinates, with origin at eta=0. With these coordinates, the length of a point P is given in fractile units - the first-stage estimate of the fraction of area enclosing all points with same or less length as that of P. A direction is specified by m=n-1 angles, where n is the number of (active) etas. At the second stage, the fraction of area enclosing all points of length less than r0 is taken to be identical to that computed at the first-stage. To obtain the fraction of area enclosing points of length greater than r0, numerical integration is used. A grid is obtained by first taking the interval [r0,r1] of the length axis and dividing this interval into nr equal subintervals. For each such subinterval, there is a region enclosing all points of length lying in the subinterval. There is yet another region (the tail region) enclosing all points of length between r1 and .9999. Within each region, the integrand is evaluated at a number of points, distributed fairly uniformly throughout the region, but with approximately the same length. The number of points is taken to be (ns**m)*(2**n), where ns is some integer. Think of ns as the number of points in a single quadrant of a 2-dimensional ellipse in n-space, and ns**m as the number of points in a orthant of m-dimensional space. 2**n is the total number of orthants. To use a different grid, include the option GRID=(nr,ns,r0,r1) on the $ESTIM record. Constraints on these numbers are nr<=100, ns<=9999, 0<r0<r1<1. nr and ns should be integers. If r1>.9999, there is no tail region. The present default values are: GRID=(3,1,.6,.9).

(See the $ESTIMATION Help Item.)

5. With minimizations that terminate with the message “MINIMIZATION SUCCESSFUL” in the report file, a warning message will now sometimes also be output:

HOWEVER, PROBLEMS OCCURRED WITH THE MINIMIZATION.
REGARD THE RESULTS OF THE ESTIMATION STEP CAREFULLY, AND ACCEPT THEM ONLY
AFTER CHECKING THAT THE COVARIANCE STEP PRODUCES REASONABLE OUTPUT.

This may happen when e.g. the final gradient vector is too large.

6. Along with the ETABAR statistic, its estimated standard error is now also output. As with Version V, the corresponding P-value that is output relates to the hypothesis that the true value of ETABAR is 0. With Version VI, this is now a default P-value, and one may obtain an alternative P-value. This is done by implementing a subsequent NONMEM problem which simulates a new data set under the model used with the current problem and obtains parameter estimates for these new data along with an ETABAR statistic, and by requesting that the P-value corresponding to this statistic relate to the hypothesis that the true value of ETABAR for the subsequent problem is the same as that with the original problem. This provides a type of check on the adequacy of the data-analytic model (but the necessity to satisfy this check, like all such checks, is subjective). To get this alternative P-value, generate a model specification file with the original problem, input this file with the subsequent problem, and include the new option ETABARCHECK on the $ESTIMATION record of the subsequent problem.

(See the $ESTIMATION and ETABAR Help Items.)

7. With Version V, there are situations where the computation of ETABAR with a mixture model was not sensible. Now multiple ETABAR vectors are output, one for each subpopulation. For each subpopulation, ETABAR is computed over all individuals who have been classified (via a standard Bayesian type computation) to be in the subpopulation. Due to possible misclassification, the default P-value may not be reliable, and it is not output. One may, however, use the alternative P-value described in point 6.

8. A simultaneous analysis of both regular-continuous and odd-type data can be implemented more easily than can be done with Version V. A flag F_FLAG can now be set in a NONMEM-PRED common NMPR17 (via NM-TRAN abbreviated code), which directs NONMEM to interpret the value set in Y (in NM-TRAN abbreviated code) or in F (in a user-supplied PRED or ERROR routine) as being either a “prediction” or a likelihood, on an observation-specific basis.

(See the NMPR17 Help Item.)

9. A frequency-prior (Gisleskog et al. 2002) may be specified for data-analytic purposes, or to simulate parameter values (in a Simulation Step). Two NONMEM utility routines help in this regard. NWPRI allows the THETA vector to be constrained by a “normal shaped prior”, and/or the OMEGA matrix to be constrained by an “inverse-Wishart shaped prior”, where both priors are specified by the user. TNPRI allows all parameters to be constrained by a prior that arises automatically from a NONMEM analysis of a prior data set.

(See the PRIOR and NWPRI and TNPRI Help Items and the NWPRI and TNPRI examples.)

1.3 9a. Implementing PRIOR routines

Implementing PRIOR routines has been made easier by a new NM-TRAN control record, $PRIOR, which may be used instead of a user-written PRIOR routine.

(See the $PRIOR help item.)

10. With an intraindividual mean-variance model (see Guide VII), the conditional likelihood for an observation is “conditioned” on values of the population parameters. With Version VI, it may easily be conditioned further; the condition being that the observation is inside or outside a given one- or two-sided interval. This is useful if e.g. blood concentrations below a known quantification limit have been eliminated from the data set. In this case one may use a one-sided interval with lower boundary equal to the quantification limit. It is assumed that the observations of the mean-variance model are normally distributed. With an observation record simply set the boundaries of the interval, YLO and YUP, in the NONMEM-PRED common NMPR12 (via NM-TRAN abbreviated code).

(See the NMPR12 and ROCM44 Help Items, and the YLO example.)

11. With an intraindividual mean-variance model (see Guide VII), the possible values of an observation may be categorized into one- and/or two-sided intervals, and the conditional likelihood for an observation may be taken to be the probability that the observation is in the category into which it is observed to fall. Categorization may be preceded by conditioning as described in point 10 above. Interval boundaries may be random variables. It is assumed that the observations of the mean-variance model are normally distributed, so that the probabilities for the categories are essentially specified by probits. So this feature allows a probit analysis of data that might otherwise be handled by a logit analysis. With an observation record, simply set the lower boundary, CTLO, and upper boundary, CTUP, (and their η-derivatives) of the observed interval in the NONMEM-PRED commons NMPR13 and NMPR14 (via NM-TRAN abbreviated code).

(See the NMPR13 and NMPR14 and ROCM45 Help Items.)

12. The individual contributions to the objective function value are now stored in ROCM50 along with the individual ID values.

(See the ROCM50 Help Item.)

1.4 12a. Elements of ROCM50 may be used in abbreviated code

See ROCM50 below.

(See the ROCM50 Help item and ROCM50 example.)

13. There is a new option on the $ESTIMATION record. The objective function value and the gradients can be calculated after the individual contributions have been sorted from the lowest to highest absolute values. In some cases this may decrease rounding errors. Simply include the SORT option on the $ESTIMATION record.

(See the $ESTIMATION Help Item.) The prototype code for sorting the individual components of the objective function value and the gradients was provided by Dr. Mark Sale, Next Level Solutions, LLC.

1.5 Default Boundary Test is User-selectable

With NONMEM VI 1.0, a new test (the Default Boundary Test) was introduced that can give rise to the error message

PARAMETER ESTIMATE IS NEAR ITS BOUNDARY THIS MUST BE ADDRESSED BEFORE
THE COVARIANCE STEP CAN BE IMPLEMENTED

With NONMEM VI 2.0, this Default Boundary Test is optional for OMEGA, SIGMA, and/or THETA.

(See the $ESTIMATION Help Item.)

2 I.2.2 General

1. There are a number of new NONMEM read-only commons, some containing particularly useful variables.

• ROCM32: the numbers of the “current” individual record (NIREC) and data record (NDREC) within the individual record
• ROCM46: the number of individual records containing at least one observation (NINDR), and the numbers of the first (INDR1) and last (INDR2) such individual records
• ROCM48: the total number of data records in the “current” individual record (LIREC)
• ROCM49: PRED_, RES_ and WRES_, the PRED, RES and WRES data items, available at problem-finalization (ICALL=3) during any pass through the data set

All these variables are available via NM-TRAN abbreviated code.

(See the ROCM32 and ROCM46 and ROCM48 and ROCM49 Help Items.)

2. Data records are passed sequentially to the PRED routine; heretofore, always from the first data record of an individual record to the last record of the individual record. With Version VI, variables (RPTI, RPTO, RPTON and PRDFL) in NONMEM-PRED common NMPR10 can be set and/or tested (via NM-TRAN abbreviated code), thus allowing a subsequence of data records to be repeatedly passed to PRED multiple times before the next data record following the last record of the subsequence is passed. Subsequences can be nested. This “repetition feature” allows e.g. kinetics to be computed by a convolution integral. The user can exercise control over which pass through a sequence it is, during which the output from PRED with a given data record will be used by NONMEM. This “PRED control” allows e.g. PRED output with a given record to involve computations over subsequent, as well as prior, records.

(See the NMPR10 Help Item and the Repetition 1 and Repetition 2 examples.)

3. The maximum number of subpopulations allowed with a mixture model is now an installation parameter in SIZES. The default value is 4, the maximum number allowed with Version V, and still, a fairly good number to use. With most data sets it is hard to discern more than 4 subpopulations. On the other hand, a recent paper (Beal 2005) describes a technique for dealing with discrete covariates with missing values, which could greatly increase the required number of mixture subpopulations. (See the SIZES Help Item.)

4. With a mixture model, mixture probabilities can be individual-specific, and so the MIX routine is called repeatedly - with each individual. Useful individual-specific information for setting the probabilities may be made available to MIX by means of the $CONTR record. With Version VI, the data record in NONMEM read-only common ROCM31 serves to provide an additional way for individual-specific information to be made available (and so now serves an additional purpose to that of providing a template data record). With each individual, when MIX is called, the first data record of the individual record may be found in the common. (See Section 4.1.6 and the ROCM31 and $CONTR and $MIX Help Items.)

5. With Version V, when a model specification file (MSF) is input, and the Simulation Step is implemented along with other NONMEM steps, the parameter values used in the simulation (the “true values”), as well as the initial parameter values used in conjunction with the other steps, are the final parameter estimates included in the MSF. With Version VI, more flexibility is available with a new option - TRUE - on the $SIMULATION record, which allows the old behavior, but this behavior is not the default behavior.

(See the $SIMULATION Help Item.)

6. At problem-finalization, a flag (SKIP_) can now be set in NONMEM-PRED common NMPR15 (via NM-TRAN abbreviated code), which permits a problem (with subproblems), or a superproblem or superproblem iteration to be ended and for NONMEM to continue with the next problem, superproblem or superproblem iteration, respectively.

(See the NMPR15 Help Item.)

7. There are a number of new NONMEM-PREDPP commons, some of which are described above. Of particular additional interest:

• NMPR16: parameter values produced during the Simulation Step when a frequency prior (see above) is used during that step

8. There is a new step, NONPARAMETRIC. NONMEM obtains either marginal cumulatives or conditional (non-parametric) estimates of etas.

(See the $NONPARAMETRIC and $SIMULATION and ROCM18 Help items, and the example.)

2.1 Increase in the number of records

Increase in the number of records that may be specified.

If the number of records in the input data set is not specified in the NONMEM control file, NONMEM reads the data set to a NONMEM FINISH record or to end of file. With all previous versions of NONMEM, the maximum number of records that could be specified in the NONMEM control file was 9999. When the NONMEM data set contained more than 9999 records, NM-TRAN put a FINISH record at the end (unless the WIDE option of $DATA was used). A FINISH record can cause problems, e.g., if another program is later used to process the NONMEM data set FDATA. With NONMEM VI 2.0, the maximum number of records that may be specified is 99999999. Only if the number of records is larger than 99999999 will a FINISH record will be present in FDATA.

(See the FINISH Help Item.)

2.2 Increase in the maximum number of function evaluations

Previously, the maximum number of function evaluations that could be specified for the the Estimation Step was 9999. Now, the maximum is 99999999. It may still be useful to specify a smaller number of function evaluations and use a Model Specification File to continue the Estimation step in a later run, but is it not necessary.

2.3 Table: PRED, RES, and WRES with NOAPPEND; more items

When the NOAPPEND option of a $TABLE record is used, PRED, RES, and WRES may be listed anywhere in a table and will appear where they are listed. Previously, they could listed anywhere in the table but would appear in the NONMEM table only as the final items (unless NOAPPEND was specified, in which case they would not appear at all).

When the NOAPPEND option of a $TABLE record is used, up to 54 items can be output in a table file. Otherwise, the limit is 50. Previously, these limits were 24 and 20, respectively.

2.4 Table and Scatter: more PRED-defined labels

Up to 50 distinct PRED-defined items may be displayed in all tables and scatters combined. Previously, the limit was 20.

Beal, Stuart L. 1984. “Population Pharmacokinetic Data and Parameter Estimation Based on Their First Two Statistical Moments.” Drug Metabolism Reviews 15 (1–2): 173–93. https://doi.org/10.3109/03602538409015064.

Beal, Stuart L. 2005. “Conditioning on Certain Random Events Associated with Statistical Variability in PK/PD.” Journal of Pharmacokinetics and Pharmacodynamics 32 (2): 213–43. https://doi.org/10.1007/s10928-005-0064-9.

Gisleskog, Paul Olof, Mats O. Karlsson, and Stuart L. Beal. 2002. “Use of Prior Information to Stabilize a Population Data Analysis.” Journal of Pharmacokinetics and Pharmacodynamics 29 (5–6): 473–505. https://doi.org/10.1023/A:1022972420004.