64  $TABLE: Additional Statistical Diagnostics, Associated Parameters, and Output Format

Introduction to NONMEM 7.6.0

1 Requesting a Range of Etas to be Outputted: Etas(x:y) (NM73)

Instead of requesting each ETA specifically in a $TABLE item list, a range of etas may be requested:

ETAS(2:4)

is equivalent to requesting ETA2, ETA3, and ETA4.

ETAS(5)

or

ETAS(5:LAST)

is equivalent to requesting ETA(5), ETA(6), … to ETA(NETAS).

As of nm74, if NPOPETAS is specified on a $MSFI record, then the value of NPOPETAS is used as LAST.

As of NM74, more flexible syntax is available. The word TO may be used in place of the semi-colon :. The BY expression may be used:

ETAS(1 TO 10 by 3) prints out etas 1, 4, 7, 10. ETAS(LAST TO 1 by -3) prints out etas 10, 7, 4, 1 (assuming LAST=10).

A number list may be given:

ETAS(1,5,12,4) prints out etas 1, 5, 12, 4. ETAS(4:1) prints etas 4, 3, 2, 1. ETAS(4:1 by -2) prints etas 4, 2. ETAS(1:4 by -1) prints etas 4, 3, 2, 1 (the by value sets the direction).

The $SCAT will also interpret this syntax, for example,

$SCAT ETAS(1:2) VS ETA3

is equivalent to

$SCAT ETA1 ETA2 VS ETA3

However, unlike $TABLE, $SCAT will ignore implied endings, such as

$SCAT ETAS(1:LAST) VS ETA3

And just interpret it as

$SCAT ETA1 VS ETA3

2 Excluding Records from Being Outputted: The EXCLUDE_BY option (NM74)

A data item or defined variable may be identified on a $TABLE record as an EXCLUDE_BY variable, which if not 0, will exclude the record. For example:

$PK
...
EXCL=0
IF(ID.GE.45.AND.ID.LE.53) EXCL=1
...
$TABLE ID TIME DV IPRED CL V1 Q V2 ETAS(1:LAST) EXCLUDE_BY EXCL NOAPPEND FILE=exctable.par NOPRINT

The table exctable.par will not list records from subjects 45 to 53. If more than one exclusion variable is listed, then if any of these have a non-zero value, the record will be excluded.

To further assist in selecting records, the following variables are available, against which you may compare with NDREC, the present record number for an individual:

• FIRSTREC: first record of subject.
• LASTREC: last record of subject.
• FIRSTOBS: first observation record of subject (for which MDV=0 or 100).
• LASTOBS: last observation record of subject (for which MDV=0 or 100).
• FIRSTDOS: First record of subject with EVID=1 or EVID=4. FIRSTDOS=-1 when there are no dose records, or PREDPP is not used.
• LASTDOS: Last record of subject with EVID=1 or EVID=4. LASTDOS=-1 when there are no dose records, or PREDPP is not used.
• EFIRSTREC: first record of subject during estimation (so, among records for which MDV=0 or MDV=1).
• ELASTREC: last record of subject during estimation (so, among records for which MDV=0 or MDV=1).
• EFIRSTOBS: first observation record of subject during estimation (so, among records for which MDV=0).
• ELASTOBS: last observation record of subject during estimation (so, among records for which MDV=0).
• EFIRSTDOS: First record of subject with EVID=1 or EVID=4 during estimation (so, among records for which MDV=1). EFIRSTDOS=-1 if no dose record PREDPP is not used.
• ELASTDOS: Last record of subject with EVID=1 or EVID=4 during estimation (so, among records for which MDV=1). ELASTDOS=-1 if no dose record PREDPP is not used.
• IRECIDX: IRECIDX+1 is starting absolute record number in the data set for the present subject (so while NDREC always sets to 1 for the first record of each subject, IRECIDX accumulates record number positions from previous subjects. For subject 1, IRECIDX+1=1, for subject 2, IRECIDX+1=LASTREC of previous subject +1).

During ICALL=1 (INFN initialization), ICALL=3 (INFN finalization), or ICALL=4 (simulation) only FIRSTREC and LASTREC are available. The other values will be set to -1.

To refer to these variables in $PK, for example, insert include the nonmem_reserved_general file statement, which contains their location information:

$PK
include nonmem_reserved_general
...
EXCL=1
IF(NDREC==LASTOBS) EXCL=0
IF(NDREC==FIRSTOBS) EXCL=0
...

$TABLE ID TIME IPRED EXCLUDE_BY EXCL NOPRINT NOAPPEND file=mytable.TAB

Then, only the record containing first and last observations will be printed to the table.

3 Selecting LASTONLY and FIRSTLASTONLY Records (NM74)

In addition to FIRSTONLY, you may select LASTONLY to just report the last record of each subject, or FIRSTLASTONLY, to report just first and last records. The FIRSTONLY, LASTONLY, and FIRSTLASTONLY options are mutually exclusive, and may not be used in combination.

4 Requesting Standard Errors to User-Defined and PREDPP Variables (NM74)

With the variance-covariance of the etas, and evaluation of the partial derivative of the particular table item with respect to eta by central finite difference method, and by principle of propagation of errors, the standard errors of the table items are evaluated. To report standard errors associated with etas (individual variables) in the tables for user defined variables, set

$TABLE ...VARCALC=1.

See setest.ctl in the examples directory. If using RFORMAT formatting system, make sure to allow enough format fields to include reported standard errors. In addition, full variances-covariances among all user-defined and PREDPP variables will be outputted to root.vpd (the FORMAT used for this file is that defined in the $EST statement). To only print to the vpd file, and not report SE’s to the table, set VARCALC=2. The variances are produced by using the variance-covariance of the etas/phis reported in the root.phi file, centered about the values reported in the table files. For classical NONMEM methods, this means variables are centered at the eta-mode assessed variable values, and approximate variance-covariance from the MAP estimation step. For EM methods:

FNLETA=1,3: variance-covariance of etas from last iteration of last estimation method is used, and centered about FNLETA based evaluations of variables (MAP, eta-mode).

FNLETA=0: variance-covariance of etas from last iteration of last estimation method is used, and centered about method based evaluations of variables (conditional means for EM methods, conditional modes for classical and ITS methods, MCMC posterior means for BAYES).

FNLETA=2: variance-covariance, and etas used for assessing variables are obtained from external source.

As of nm75, table root.vpt is also created, which incorporates variance-covariances associated with THETAs (or OMEGAs and SIGMAs, listed in the .cov file) as well as those associated with etas/individual variables (phc() listed in the .phi file). To present the comparable total standard errors in a user defined table, set $TABLE ...VARCALC=3. The standard errors in the .vpt table are most accurately assessed if the problem is MU modeled. If it is not MU modeled, then only standard errors of variables that do not vary with etas (variables derived exclusively from THETAs, for example), are assured to be correctly calculated. It also takes into account that the .phi tables contain variances to phi’s and etas under fixed theta conditions for all methods except BAYES, where the variances to phis or etas are sampled under varying theta conditions.

Note that to obtain the PRED and its error, you should define a variable that is based on IPRED, but evaluated only when COMACT=1 (NONMEM passes through the model system with ETA=0):

IPRED=F
If(COMACT==1) PREDU=IPRED

So, PREDU is IPRED evaluated with ETA=0. Then, list this in a table:

$TABLE ... PREDU ... VARCALC=3

Keep in mind that PREDU has only statistical errors associated with THETAS (and SIGMAS, depending on how PREDU is evaluated, based on the variance-covariance in the .cov file), and not etas (the phc() variances in the .phi file do not contribute to PREDU, since etas=0 when COMACT=1), so VARCALC=3 is required. No statistical errors are evaluated for weighted residual diagnostics (WRES, NPD’s, etc).

The vpd and vpt files contain entries for each COMACT the NONMEM used to calculate the table items. For the vpd file, variance-covariance among user-defined variables and PREDPP variables are placed in this file, associated with individual variances of the .phi file. The COMACT column value is usually 2, or 3 if $NONPARAMETRIC ETAS was implemented. The subsequent columns contain the values of user-defined table items, followed by variance-covariance of these items (labeled VAR1_VAR2, where VAR1 and VAR2 are the particular table items).

For the vpt file, the total variance-covariance among user-defined variables and PREDPP variables are placed in this file, associated with individual variances of the .phi file, as well as the variance-covariance of the fixed effects variables (THETAS, SIGMAS, OMEGAS) of the .cov file. The COMACT column value is 1 to indicate the record contains the values during a COMACT=1 pass, that is with ETAs set to 0. Thus, an individual predicted value (IPRED) on this record would in fact by the population predicted value (PRED). The subsequent columns contain the values of user-defined table items, followed by variance-covariance of these items (labled VAR1_VAR2, wher VAR1 and VAR2 are the particular table items). The variance IPRED_IPRED for a record with COMACT=1 would then have essentially the statistical error associated with the variance-covariance of the appropriate THETAs. If you are satisfied to pick up the PRED values and their errors from the vpt table, then it is not necessary to define a COMACT=1 evaluated variable in the control stream, as you can get it from IPRED on the COMACT=1 rows in the .vpt table. But if you want the appropriate PRED variable with standard error presented in the user-requested $TABLE table, then the PREDU must be defined as above. For non-parametric estimated etas, when $NONPARAMETRIC ETAS is used, then a COMACT 3 pass will be made, and their values will be presented in the vpt table with the COMACT set to 3. See Guide VIII for more on COMACT.

See the example ..\examples\predu.ctl, which prints out IPRED, PRED, and PREDU in the various tables (predu.tab, predu.phi, predu.cov, predu.vpt, predu.vpd), and you can compare their values and associated standard errors.

5 New diagnostic items

Additional types of pred, res, and wres values may be requested than the usual set available in NONMEM VI. They may be specified at any $TABLE command or $SCATTER command, as one would request PRED, RES, or WRES items. If $TABLE statements succeed multiple $EST statements within a run, the table results (as well as scatter plots if requested via $SCATTER) will pertain to the last analysis.

5.1 OBJI

These are objective function values for each individual. The sum of the individual objective function values is equal to the total objective function.

5.2 NPRED, NRES, NWRES

These are non-conditional, no eta-epsilon interaction, pred, res, and wres values. These are identical to those issued by NONMEM V as PRED, RES, and WRES.

5.3 PREDI, RESI, WRESI

These are non-conditional, with eta-epsilon interaction, pred, res, and wres values. These are identical to those issued by NONMEM VI as PRED, RES, and WRES. The WRESI will not differ from NWRES if INTERACTION was not selected in the previous $EST command.

5.4 CPRED, CRES, CWRES

These are conditional, no eta-epsilon interaction, pred, res, and wres values as described in (Hooker et al. 2007). The conditional mode etas (from FOCE or ITS, also known as conditional parametric etas (CPE), empirical bayes estimates (EBE), posthoc estimates of etas, or mode a posteriori (MAP) estimates) or conditional mean etas (from Monte Carlo EM methods) will be referred to as η̂ (eta hat), must be available from a previous $EST MAXEVAL>0 command. The conditional weighted residuals are estimated based on a linear Taylor series approximation that is extrapolated from the conditional mean or mode (or posthoc) eta estimates, rather than about eta=0:

\[ \text{CPRED}_{ij} = f_{ij}(\hat{\eta}) - g_{ij}(\hat{\eta})\hat{\eta} \]

using the nomenclature of Guide I, Section E2. Then

\[ \text{CRES}_{ij} = y_{ij} - \text{CPRED}_{ij} \]

The population variance covariance of observed data described in Guide I, E.2 is also evaluated at eta_hat: \(C_i(\hat{\eta})\):

\[ \text{CWRES}_i = C(\hat{\eta})_i^{-1/2} (y_i - \text{CPRED}_i(\hat{\eta})) \]

Because of the linear back extrapolation, it is possible for some CPRED values to be negative. Users may prefer to request NPRED CRES CWRES, or NPRED RES CWRES. The conditional weighted residual will not differ from the non-conditional weighted residual if FO was selected in the previous $EST command.

In NM72, if $EST INTERACTION was not specified prior to requesting $TABLE CWRES, then the population variance-covariance is evaluated at eta=0: \(C_i(\eta=0)\). In NONMEM 7.1.0 and 7.1.2, regardless of INTERACTION setting in a previous $EST statement, \(C_i(\hat{\eta})\) is used.

5.5 CPREDI, CRESI, CWRESI

These are conditional, with eta-epsilon interaction, pred, res, and wres values. The conditional mode or conditional mean etas must be available from a previous $EST MAXEVAL>0 command.

5.6 EPRED, ERES, EWRES

The EPRED, ERES, EWRES are Monte-Carlo generated (expected, or exact) pred, res, and wres values, and are not linearized approximations like the other diagnostic types.

The expected diagnostic items are evaluated using predicted function and residual variances evaluated over a Monte Carlo sampled range of etas with population variance Omega, and assuming asymptotic values for averaging over the y domain. Define

\[ \text{EPRED}_{ij} = \int_{-\infty}^{\infty} f_{ij}(\eta) p(\eta | 0, \Omega) d\eta \]

is the expected predicted value for data point j of subject i for a given subject, evaluated by Monte Carlo sampling, overall possible eta. The probability density of eta:

\[ p(\eta | 0, \Omega) d\eta \]

is a multivariate normal distribution with eta variance Ω. The 1×n_i vector of EPRED for a given subject, where n_i is the number of data points to that subject, is then:

\[ \text{EPRED}_i = \int_{-\infty}^{\infty} f_i(\eta) p(\eta | 0, \Omega) d\eta \]

Then the corresponding residual vector for observed values y_i is

\[ \text{ERES}_i = y_i - \text{EPRED}_i \]

The residual (epsilon) variance matrix using the nomenclature in Guide I, Sections E.2 may be

\[ V_i(\eta) = \text{diag}(h_i(\eta) \Sigma h_i(\eta)') \]

or it may be the more complicated form described in section of E.4 in the case of L2 data items. Then, the expected residual (epsilon) variance (assessed by Monte Carlo sampling) is

\[ EV_i = \int_{-\infty}^{\infty} V_i(\eta) p(\eta | 0, \Omega) d\eta \]

The full variance-covariance matrix of size n_i×n_i, that includes residual error (epsilon) and inter-subject (eta) variance contributions is:

\[ EC_i = EV_i + \int_{-\infty}^{\infty} (f_i(\eta) - \text{EPRED}_i)(f_i(\eta) - \text{EPRED}_i)' p(\eta | 0, \Omega) d\eta \]

And is the expected population variance, Monte Carlo averaged over all possible eta. Then, following the Guide I, section E nomenclature, the population weighted residual vector for subject i is:

\[ \text{EWRES}_i = EC_i^{-1/2} \text{ERES}_i \]

where the square root of a matrix is defined here by default as evaluated by diagonalizing the matrix, and multiplying its eigenvector matrices by the square roots of the eigenvalues. Selecting the WRESCHOL option obtains the square root of the matrix by Cholesky decomposition.

5.7 ECWRES

ECWRES is a Monte Carlo assessed expected weighted residual evaluated with only the predicted function evaluated over a Monte Carlo sampled range of etas with population variance Omega, while residual variance V is always evaluated at conditional mode (from the most recent FOCE/ITS estimation) or conditional mean (from the most recent IMP/IMPMAP/SAEM analysis) eta (η̂), so that

\[ ECC_i = V_i(\hat{\eta}) + \int_{-\infty}^{\infty} (f_i(\eta) - \text{EPRED}_i)(f_i(\eta) - \text{EPRED}_i)' p(\eta | 0, \Omega) d\eta \]

and

\[ \text{ECWRES}_i = ECC_i^{-1/2} \text{ERES}_i \]

As with CWRES, the eta_hat (conditional mode or mean) values must be available from a previous $EST MAXEVAL>0 command.

Thus, ECWRES is the Monte Carlo version of CWRES, while EWRES is the Monte Carlo version of CWRESI.

In NM72, if $EST INTERACTION was not specified prior to requesting $TABLE CWRES, then the residual variance is evaluated at eta=0: \(V_i(\eta=0)\). In NONMEM 7.1.0 and 7.1.2, regardless of INTERACTION setting in a previous $EST statement, \(V_i(\hat{\eta})\) is used.

5.8 NPDE

The NPDE is the normalized prediction distribution error (reference (Comets et al. 2008): takes into account within-subject correlations), also a Monte Carlo assessed diagnostic item. For each vector of simulated data y_ki:

\[ \text{ESRES}_{ki} = y_{ki} - \text{EPRED}_i \]

its decorrelated residual vector is calculated:

\[ \text{ESWRES}_{ki} = EC_i^{-1/2} \text{ESRES}_{ki} \]

and compared against the decorrelated residual vector of observed values EWRES_i such that

\[ \text{pde}_i = \frac{1}{K} \sum_{k=1}^{K} \delta(\text{EWRES}_i - \text{ESWRES}_{ki}) \]

For K random samples, where

\[ \delta(x) = 1 \text{ for } x \geq 0 \] \[ \delta(x) = 0 \text{ for } x < 0 \]

For each element in the vector. Then, an inverse normal distribution transformation is performed:

\[ \text{npde}_i = \Phi^{-1}(\text{pde}_i) \]

5.9 NPD

The NPD is the correlated normalized prediction distribution error (reference (Brendel et al. 2006): does not take into account within-subject correlations), also a Monte Carlo assessed diagnostic item. For each vector of observed data y_i of subject i, and vector of simulated etas η_k:

\[ \text{IWRES}_{ki} = V(\eta_k)_i^{-1/2} (y_i - f_i(\eta_k)) \]

These are then averaged over all the random samples:

\[ \text{pd}_i = \frac{1}{K} \sum_{k=1}^{K} \delta(\text{IWRES}_{ki}) \]

Then, an inverse normal distribution transformation is performed:

\[ \text{npd}_i = \Phi^{-1}(\text{pd}_i) \]

The default PRED, RES, and WRES will be given the same values as PREDI, RESI, and WRESI, when INTERACTION in $EST is specified, or NPRED, NRES, and NWRES when INTERACTION in $EST is not specified.

As the PRED, RES, and WRES, may be referenced in a user-supplied $INFN routine, or in $PK or $PRED (when ICALL=3) as PRED_, RES_, WRES_, so the additional variables may be referenced by their names followed by _ (for example EWRES_).

5.10 NPDTYPE=0 (default)

If NPDTYPE=1 is set as a $TABLE option, then the strict stochastic (Monte Carlo) method over the data y domain as well as etas, as referenced in (Comets et al. 2008; Nguyen et al. 2012), is implemented for NPD diagnostics:

\[ \text{pd}_i = \frac{1}{K} \sum_{k=1}^{K} \delta(y_{ki} - y_i) \]

\[ \delta(x) = 1 \text{ for } x \geq 0 \] \[ \delta(x) = 0 \text{ for } x < 0 \]

Where y_ik is the kth simulated vector of data, and y_i is the observed data. Similarly for npde data:

\[ \text{EPRED}_i = \frac{1}{K} \sum_{k=1}^{K} y_{ki} \]

\[ \text{ERES}_i = y_i - \text{EPRED}_i \]

\[ \text{ESRES}_{ki} = y_{ki} - \text{EPRED}_i \]

\[ EV_i = \frac{1}{K} \sum_{k=1}^{K} (y_{ik} - \text{EPRED}_i)(y_{ik} - \text{EPRED}_i)' \]

\[ \text{EWRES}_i = EV_i^{-1/2} \text{ERES}_i \]

\[ \text{ESWRES}_i = EV_{ki}^{-1/2} \text{ESRES}_i \]

\[ \text{pde}_i = \frac{1}{K} \sum_{k=1}^{K} \delta(\text{EWRES}_i - \text{EWRES}_{ik}) \]

5.11 INTERPTYPE=0 (default) (nm751)

By default the interpolation for NPD/NPDE evaluations type is 0, that is, the nearest value in the quantile area is chosen, without interpolation. If INTERPTYPE=1, then an interpolation is performed using an inverse NORM linearized interpolator. In addition, for NPDE it will extrapolate quantile positions that are outside the Monte Carlo sampled range, allowing severe outliers to be less truncated. For example, by default, if the EXAMPLE=300, and a quantile position is 0.001, then the 1/300 value is returned, the lowest value from 300 samples. With INTERPTYPE=1, an extrapolation is performed to better represent the 0.001 value, and will be less than the lowest value among 300 samples simulated. The INTERPTYPE has less impact with increasing ESAMPLE.

5.12 CIWRES, CIPRED, CIRES, CIWRESI (NM72)

The CIWRES is the conditional individual weighted residual as evaluated during the estimation, equivalent to (DV-F)/(F*SQRT(SIGMA(1,1))) for simple problems with proportional residual error. With L2 data or CORRL2 data, the individual weighted residuals are in their decorrelated forms:

\[ \text{CIWRES}_i = V(\hat{\eta})_i^{-1/2} (y_i - f_i(\hat{\eta})) \]

when INTERACTION in the previous $EST record is set, and a conditional analysis (non-FO) was performed. For individual i, where individual residual variance matrix V_i and individual vector of predicted \(f_i(\hat{\eta})\) are evaluated at the conditional mode or mean eta (designated as eta hat). The square root of the matrix V_i may be evaluated by using the square root of the eigenvalues, or by Cholesky decomposition when WRESCHOL option is used (see below). Similarly, the CIPRED is the individual predicted value \(f_i(\hat{\eta})\) at the conditional mode or mean eta, and CIRES=DV-\(f_i(\hat{\eta})\).

When INTERACTION is not set, then

\[ \text{CIWRES}_i = V(\eta=0)_i^{-1/2} (y_i - f_i(\hat{\eta})) \]

is evaluated, that is, the variance portion is evaluated using \(f_i(\eta=0)\). However CIWRESI (conditional individual weighted residual with interaction) is always evaluated as (except for FO, see below)

\[ \text{CIWRESI}_i = V(\hat{\eta})_i^{-1/2} (y_i - f_i(\hat{\eta})) \]

regardless of the INTERACTION setting.

For FO, the conditional individual weighted residual will not differ from the non-conditional weighted residual. That is, for FO, the CIWRES and CIPRED are evaluated using F(eta=0) for numerator and denominator terms, since this is what is done during estimation, and no EBE (eta-hat) is evaluated:

\[ \text{CIWRES}_i = V(\eta=0)_i^{-1/2} (y_i - f_i(\eta=0)) = \text{CIWRESI}_i \]

Even for FO with interaction, the predicted function (numerator) and residual variance (denominator) is still evaluated at eta=0, so CIWRESI=CIWRES. The interaction contribution is accounted for with additional first-order Taylor terms to make a linear projection of the contribution of eta-eps interaction. While it would be inappropriate to add these Taylor terms to CIWRESI, these Taylor terms are added to the population residual assessment WRESI, hence WRESI will differ from NWRESI with FO INTERACTION.

There are other individual residual values available, mostly as place holders in the system, but these have no additional statistical value. They are:

NIPRED=IPREDI=IPRD=NPRED
CIPREDI=CIPRED
EIPRED=EPRED

NIRES=IRESI=NRES=IRS
CIRESI=CIRES
EIRES=ERES

\[ \text{NIWRES}_i = V(\eta=0)_i^{-1/2} (y_i - f_i(\eta=0)) \]

IWRESI=NIWRES=IWRS

\[ \text{EIWRES}_i = \int_{-\infty}^{\infty} V(\eta)_i^{-1/2} (y_i - f_i(\eta)) p(\eta | 0, \Omega) d\eta \]

5.13 MDVRES=0 (NM73) (default)

Set MDVRES to 1 in the $ERROR or $PRED routine if you do not want to include a particular value for weighted residual assessment. This may be useful when, for example, this data point is assessed by a non-normal distribution likelihood such as the PHI() function for below detection limit values, in which F_FLAG is set. By default, if at least one data value of a given subject is fitted with a non-normal distribution likelihood, then population weighted residual diagnostics are not assessed for any of the data for that subject. By setting MDVRES=1 to these particular below detection values, the weighted residual algorithm can assess the remaining normally distributed values for that subject. For example,

$ERROR
SD = THETA(5)
IPRED = LOG(F)
DUM = (LOQ - IPRED) / SD
CUMD = PHI(DUM)
IF (TYPE .EQ. 1) THEN
          F_FLAG = 0
          Y = IPRED + SD * ERR(1)
ENDIF
IF (TYPE .EQ. 2) THEN
          F_FLAG = 1
          Y = CUMD
          MDVRES=1
ENDIF

MDVRES stands for missing data value (MDV) for residual (RES) assessment. Setting MDVRES to 1 is equivalent to temporarily declaring that data point as missing during the weighted residual assessments.

6 Handling Level of Quantitation (LOQ) and Level Above Quantitation (LAQ) Data, and their Incorporation in to NPDE and NPD Assessments (nm73, nm74)

To incorporate LOQ data into NPDE assessments (Nguyen et al. 2012), use the following method (as an example):

Here, TYPE and LOQ are user-defined in previous code, or data item (..\examples\loq\ad3tr4_loq0).

$ERROR
SD = THETA(5)
IPRED = LOG(F)
DUM = (LOQ - IPRED) / SD
CUMD = PHI(DUM)
IF (TYPE .EQ. 1.OR.NPDE_MODE.EQ.1) THEN
          F_FLAG = 0
          Y = IPRED + SD * ERR(1)
ENDIF
IF (TYPE .EQ. 2.AND.NPDE_MODE.EQ.0) THEN
          F_FLAG = 1
          Y = CUMD
          MDVRES=1
ENDIF
IF(TYPE.EQ.2) DV_LOQ=LOQ

By default, DV_LOQ is set to -1.0d+300 by the NONMEM routine that calls ERROR/PRED. If the user’s ERROR/PRED sets DV_LOQ to some other value and NPDE_MODE=1, then the NPDE is being evaluated during that time, and this censored value is to be treated as if it is a non-censored datum with value of LOQ (DV_LOQ=LOQ), in accordance with (Nguyen et al. 2012), utilizing a standard F_FLAG=0 definition for Y. Note that during estimation of the objective function (when NPDE_MODE=0), NPDE is not being evaluated, and censored values should be treated using F_FLAG=1, and Y must be defined as the integral of the normal density from -inf to LOQ.

New in nm743, you can specify an above quantifiable limit with the reserved variable DV_LAQ as well. For example, (..\examples\loq\ad3tr4a_loq0):

$ERROR
SD = THETA(5)
IPRED = LOG(F)
DUM = (LOQ - IPRED) / SD ; LOQ=lower level of quantifiable level
CUMD = PHI(DUM)+1.0E-10
DUMA = (LAQ - IPRED) / SD ; LAQ=Upper (above) quantifiable level
CUMDA = PHI(DUMA)-1.0E-10
IF(TYPE.EQ.2) DV_LOQ=LOQ ; Type=2 is concentration below LOQ
IF(TYPE.EQ.3) DV_LAQ=LAQ ; Type=3 is concentration above LAQ
IF (TYPE .EQ. 1.OR.NPDE_MODE==1) THEN
          F_FLAG = 0
          Y = IPRED + SD * ERR(1)
ENDIF
IF (TYPE .EQ. 2.AND.NPDE_MODE==0) THEN
          F_FLAG = 1
          Y = CUMD
          MDVRES=1
ENDIF
IF (TYPE .EQ. 3.AND.NPDE_MODE==0) THEN
          F_FLAG = 1
          Y = (1.0-CUMDA)
          MDVRES=1
ENDIF

New in nm74, for use with NPD, the user may supply the cumulative distribution function using the reserved variable CDF_L. For example, in a general likelihood modeled problem, essentially the previous example, but the Y values of all data are returned in their -2LL form (..\examples\loq\ad3tr4_loq6):

$ERROR
SD = THETA(5)
IPRED = LOG(F)
DUM2 = (DV - IPRED) / SD
DUM = (LOQ - IPRED) / SD
CUMD = PHI(DUM)+1.0E-30
CUMD2 = PHI(DUM2)+1.0E-30
IF(TYPE.EQ.1) THEN
Y=2.0*LOG(SD)+DUM2*DUM2
CDF_L=CUMD2
ENDIF
IF(TYPE.EQ.2) THEN
Y = -2.0*LOG(CUMD)
CDF_L=CUMD
DV_LOQ=LOQ
ENDIF
...
$EST METHOD=COND LAPLACE -2LL MAXEVAL=9999 NSIG=3 SIGL=9 SIGLO=9 PRINT=5 NOABORT MCETA=5
$TABLE ID TIME DV IPRED NPD NOAPPEND ONEHEADER ESAMPLE=1000
 FILE=ad3tr4_loq6.TAB NOPRINT

Note that only NPD can be evaluated without consideration of EWRES and EPRED constructs. NPDE, EWRES and EPRED cannot be evaluated for general non-normal likelihood data, or when normal data are modeled in a general log-likelihood manner, rather than assumed normal in the usual way. The CDF reserved variable associated with above quantitation level DV_LAQ is CDF_LA.

The DV_LOQ/CDF_L and DV_LAQ/CDF_LA reserved variables may be also used for evaluating NPD diagnostics for categorical data. In such cases, the DV_LAQ/CDF_LA variables must define the lower bound of a categorical level, while the DV_LOQ/CDF_L defines the upper boundary of a level, which is counter-intuitive, and its use is demonstrated as follows. In the following example, some data are normally distributed, and others are binomial (categorical). The NPDE will be evaluated only for those that are normal, while NPD are evaluated for both types of data. The CDF_L indicates to the diagnostics routine that this datum is to be treated as non-normal, with a cumulative distribution value of CDF_L, which it can use for evaluating the NPD. Because the probability is categorical, the lower bound CDF_LA needs also to be given, to map the probability of having data value DV (=0 or 1) be between CDF_LA and CDF_L for a random uniform variable (..\examples\loq\example10lcdf).

$ERROR
EXCL2=1.0-TYPE
EXCL=TYPE
EXCL3=0.0
IF(EVID/=0) EXCL=1.0
IF(EVID/=0) EXCL2=1.0
IF(EVID/=0) EXCL3=1.0
    EXPP=THETA(4)+F*THETA(5)
IPRED=F
; Use protected exponent PEXP, to avoid numerical overflow
A=PEXP(EXPP)
B=1.0+A
IF (TYPE.EQ.0.OR.NPDE_MODE==1) THEN
; PK Data
    F_FLAG=0
    Y=F+F*ERR(1) ; a prediction
ELSE
; Categorical data
    F_FLAG=1
    Y=DV*A/B+(1.0-DV)/B       ; a likelihood
    MDVRES=1
ENDIF
IF(TYPE==1) THEN
CDF_L=(1.0-DV)*1.0/B + DV
CDF_LA=DV*1.0/B
DV_LOQ=DV
DV_LAQ=DV-1.0
ENDIF

For DV=0 the lower bound of its value is DV_LAQ=-1 with a cumulative probability of CDF_LA=0 (that is, it is impossible for DV to have a value of -1), and upper bound DV_LOQ=0 has probability of CDF_L of 1/B (that is, the probability of it being 0 is 1/B, which is consistent with how the likelihood Y is defined). When DV=1, then range of its value is from lower bound DV_LAQ=0 with cumulative probability CDF_LA=1/B (which is consistent with the upper bound interpretation for DV=0), to upper bound DV_LOQ=1 with cumulative probability CDF_L= of 1 (consistent with the idea that the highest DV can be is 1).

When tabulated, we have the following:

DV	Lower Bound (DV_LAQ)	Upper Bound (DV_LOQ)	Lower Cumulative (CDF_LA)	Upper Cumulative (CDF_LOQ)
0	-1	0	0	1/B
1	0	1	1/B	1

The common thread of logic is to direct the NPD algorithm to create random positioned NPD values for any data that are below DV_LOQ, and or above DV_LAQ. For continuous data, values are accurately known as the observed value DV between DV_LOQ and DV_LAQ, but are not known below DV_LOQ, or above DV_LAQ, and so, the NPD should create a scrambled DV value whenever it is below DV_LOQ and/or above DV_LAQ. Categorical data have discrete values, so in order to turn them into continuous, normally distributed equivalents for NPD, scrambled values must be created, for values above DV_LAQ, and/or below DV_LOQ, that is, for all legitimate discrete values. Thus, when DV_LOQ<DV_LAQ, the DV_LOQ and DV_LAQ serve as the boundaries for the outer distribution region for scrambling, and when DV_LAQ<DV_LOQ, the DV_LAQ and DV_LOQ serve as boundaries for the inner distribution region for scrambling.

7 ESAMPLE=300

Number of random samples to be used to generate a Monte-Carlo based set of EPRED, ERES, ECWRES, NPDE (via EC), and EWRES. ESAMPLE should be specified only on the first $TABLE command. By default, ESAMPLE=300.

The ESAMPLE size should be at least the size of the largest number of individual observations (N) to be encountered among the subjects. It would be good to have ESAMPLE>=N+300, for example. This is because N equals the dimension size of the empirical random variance components (EC, EWRES) used for decorrelating, which will have a rank of min(ESAMPLE,N).

8 WRESCHOL (NM73)

Normally, population and individual weighted residuals are evaluated by square root of the eigenvalues of the population or individual residual variance. However, an alternative method is to Cholesky decompose the residual variance (suggested by France Mentre, personal communication), by entering the WRESCHOL option. This should be specified only on the first $TABLE command. The Cholesky form has the property of sequentially decorrelating each additional data point in the order of the data set.

9 SEED

Specify starting seed for Monte Carlo evaluations of EPRED, ERES, EWRES, ECWRES, and NPDE. The default seed is 11456. SEED should be specified only on the first $TABLE command.

10 CLOCKSEED=0 (default) (nm75)

As of nm75, the actual starting seed will be 10000*(seconds after midnight)+SEED (SEED may be set to 0 for this purpose), if CLOCKSEED=1. This allows a control stream to produce different stochastic results for automated replications, without the need to modify the seed value in the control stream file in each replication.

11 RANMETHOD=[n|S|m|P] (default n=3) (NM72)

By default, the random number generator used for Monte Carlo simulations of weighted residual items is ran3 of reference (Press et al. 1992). We feel this is the best random number generator for many purposes. However, you may choose alternative random number generators as follows:

• 0: ran0 of reference (Press et al. 1992), minimal standard generator.
• 1: ran1 of reference (Press et al. 1992), Bays and Durham.
• 2: ran2 of reference (Press et al. 1992).
• 3: ran3 of reference (Press et al. 1992), Knuth.
• 4: NONMEM’s traditional random number generator, used as default in $SIMULATION.

RANMETHOD should be specified only on the first $TABLE command. The RANMETHOD set in the $TABLE command does not propagate to $EST or $CHAIN.

As of NM73, the Sobol sequences with scrambling may be requested:

RANMETHOD=[n|S|m|P]

where n is the random number generator type, S is Sobol sequence, and m is the Sobol scrambler, and P may be specified to retain random number sequences for each subject, so that the random number sequence is retained regardless of single or parallel processing. See the description of RANMETHOD under I.38 Monte Carlo Importance Sampling EM.

Among the Sobol sequence methods, the S2 method appears to provide the least biased random samples, that is nearly uniform distribution, with good mixing in multi-dimensional spaces.

12 NOLABEL (NM73)

Do not print column labels. It may be combined with ONEHEADER to print only the title at the beginning of each table.

13 NOTITLE (NM73)

Do not print table titles. It may be combined with ONEHEADER to print only the column labels at the beginning of each table. NOLABEL NOTITLE is equivalent to NOHEADER.

14 ONEHEADERALL or ONEHEADERPERFILE (NM74)

Print header information only once per file. If a file already exists when the problem was started, and the FORWARD option was used, results are appended. In such cases, NONMEM assumes that no additional headers should be added to this file.

15 FORMAT=s1PE11.4 (default) (NM75)

This option defines the delimiter and number format for the present table, and subsequent tables, until a new FORMAT is specified. The first character defines the delimiter, which may be s for space, t for tab, comma (,) for a comma delimited file with aligned fields (so, padded with spaces), c for comma delimited file with no spaces, or q for comma delimited file with no spaces and double quotes around column names that have commas in their names (such as "OMEGA(2,1)").

The syntax for the number format is Fortran based, as follows.

For E field:

xPEw.d

indicates w total characters to be occupied by the number (including decimal point, sign, digits, E specifier, and 2 digit magnitude), d digits to the right of the decimal point, and x digits to the left of the decimal point.

Examples:

E12.5: -0.12345E+02
2PE13.6: -12.12345E+02

If you are outputting numbers that are less than 1.0E-99, such as 1.22345E-102, there will be one less significant digit displayed to make room for the extra digit in the exponent. To make room for a three digit exponent, you may set the format as follows:

xPEw.dEe

where e is the number of digits to be provided for the exponent. For example

1PE12.4E3: -2.3456E+002

For F field:

Fw.d

indicates w total characters to be occupied by the number (including decimal point, sign and digits), d digits to the right of the decimal point.

Examples:

F10.3: -0.012, 234567.123

For G field:

xPGw.d

For numbers >=0.1, will print an F field number if the value fits into w places showing d digits, otherwise will resort to xPEw.d format. For numbers <0.1, will always use xPEw.d format.

If the user-defined format is inappropriate for a particular number, then the default format will be used for that number.

An example $TABLE record could be:

$TABLE ID CMT EVID TIME NPRED NRES PREDI RESI WRESI CPRED CRES CWRES CPREDI
            CRESI CWRESI=ZABF EPRED ERES EWRES PRED RES WRES NPDE=PDERR ECWRES
            NOPRINT NOAPPEND FILE=myfile.tab ESAMPLE=1000 SEED=1233344

If you specify FORMAT=QCSV then this is equivalent to FORMAT=q1PG23.16. If you specify FORMAT=CSV is equivalent to FORMAT=c1PG23.16.

16 IDFORMAT=I (NM75)

This secifies the format (do not include the delimiter) for the ID column. By default the ID column has the same format as specified by FORMAT. However, sometimes you wish the ID to appear as an integer, in which case, you may set IDFORMAT as I.

Some examples:

• IDFORMAT=I: Integer value, left adjusted in the field.
• IDFORMAT=I6: Integer value, right adjusted in the first 6 characters of the field.
• IDFORMAT=F6.1: Floating value, with single digit to the right of the decimal.

If an improper format is given, it defaults to that of FORMAT.

17 LFORMAT, RFORMAT (NM72)

An alternative format description to FORMAT is RFORMAT and LFORMAT. RFORMAT (where R=real numbers) describes the full numeric record of a table, so that formats for specific columns may be specified. LFORMAT (where L=label) specifies the format of the full label record of a table. The formats must be enclosed in double quotes, and (), and have valid Fortran format specifiers. The RFORMAT and LFORMAT options can be repeated if the format specification is longer than 80 characters. Multiple RFORMAT and LFORMAT entries will be concatenated to form a single format record specification. For example,

LFORMAT="(4X,A4,4(',',4X,A8))"
RFORMAT="(F8.0,"
RFORMAT="4(',',1PE12.5))"

Will result in the following formats submitted to a Fortran write statement:

LFORMAT=(4X,A4,4(',',4X,A8))

for the table’s label record, and

RFORMAT=(F8.0,4(',',1PE12.5))

For the table’s numeric records. If RFORMAT and LFORMAT are given, then the FORMAT option will be ignored. By default, FORMAT, RFORMAT, LFORMAT specifications will be passed on to the next $TABLE record in a given problem unless new ones are given. To turn off an RFORMAT/LFORMAT specification in a subsequent table (and therefore use FORMAT instead), set

LFORMAT="NONE"
RFORMAT="NONE"

Here is an example of $TABLE statements designated in a control stream file:

$TABLE ID TIME PRED RES WRES CPRED CWRES EPRED ERES EWRES NOAPPEND ONEHEADER
 FILE=tabstuff.TAB NOPRINT,FORMAT=,1PE15.8

$TABLE ID CL V1 Q V2 FIRSTONLY NOAPPEND NOPRINT FILE=tabstuff.PAR
            LFORMAT="(4X,A4,4(',',4X,A8))"
            RFORMAT="(F8.0,"
            RFORMAT="4(',',1PE12.5))"

$TABLE ID ETA1 ETA2 ETA3 ETA4 FIRSTONLY NOAPPEND NOPRINT
 FILE=tabstuff.ETA,FORMAT=",F12.4"
            LFORMAT="NONE"
            RFORMAT="NONE"

There is no NMTRAN error checking on the RFORMAT and LFORMAT records, so the user must engage in trial and error to obtain a satisfactory table output. You should set MAXEVAL=0 or MAXEVAL=1 for the $EST step to do a quick check, so you don’t spend hours on estimation only to find the RFORMAT/LFORMAT were not appropriate, or use the model specification file system (MSFO/$MSFI) to preserve the estimation results and resume table output at the final estimates in a new control stream file.

A word of caution. The FORMAT descriptor 1P, which means move the decimal point to the left by 1, will be in effect for all remaining FORMAT components. For example, in

RFORMAT="(F8.0,37(',',1PE13.6),24(',',F7.2))"

the F field format that follows an E field format, in which 1P was used, will also have the decimal placed to the left, and a 1.00 would appear as a 10.00. To prevent this from occurring, revert to no decimal shift with 0P:

RFORMAT="(F8.0,37(',',1PE13.6),24(',',0PF7.2))"

18 PARAFILE (NM74)

As of NONMEM 7.4, computations of weighted residual diagnostics that are evaluated for table output are performed in a parallel computing setting, if parallel computing was requested, by the -parafile option on the command line, for example (see I.75 Parallel Computing (NM72)). If you wish to turn parallel computing off for the weighted residual computation, then set parafile to off on the first $TABLE record in the problem:

$TABLE PARAFILE=OFF

Remember that parallelization remains OFF until you turn it back ON with a PARAFILE option in a $EST, $COV or $TABLE record in subsequent problems. Also, use $TABLE RANMETHOD=P to assure that the sequence of generated random numbers remain consistent for evaluation of NPDE, etc., for repeated executions of the control stream file.

19 NOSUB=0 (Default) (NM74)

Subscripts of etas may be replaced with meaningful aliases, such as ETA(CL) in place of ETA(1). See I.8 Expansions on Abbreviated and Verbatim Code and Other Items (NM72,NM73,NM74,NM75) under $ABBR REPLACE feature for abbreviated code (NM73-NM75) for how the $ABBR REPLACE command can be used to do this. These aliases will also be used as labels for the column identifiers in tables produced by $TABLE. If you like the convenience of using meaningful labels in your abbreviated code, but do not wish to see them expressed in the table files for a given table, then set

$TABLE NOSUB=1

to turn this symbolic label substitution off. Scatter plots also are labeled with the aliases, which can be turned off or on for a given plot such as:

$SCAT NOSUB=1 ETA(V1) VS ETA(CL)
$SCAT NOSUB=0 ETA(V2) VS ETA(CL)

To set the default NOSUB for the entire problem, specify the following $DEFAULT (or $DEFAULTS) record:

$DEFAULT NOSUB=1

NOSUB also accepts a -1 to indicate revert to NONMEM default, which is to treat -1 as a 0. The NOSUB default setting is in effect throughout the control stream file, until another $DEFAULT record setting NOSUB is encountered.

20 FIXEDETAS=(number-list) (NM74)

It may be desired to treat certain etas, particularly $LEVEL etas that span groups of subjects, as if they were a fixed effect (theta) when evaluating populations characteristics during the $TABLE step, such as PRED, CWRES, NPDE, etc. In this way, the PRED evaluated will be, not of the total population, but of a given site level for that subject. For example, FIXEDETAS=(3-6,10-12) indicates etas 3 through 6, and 10 through 12 are to be placed at the Empirical Bayes position, while the other etas are set to 0, during evaluation of PRED, RES, WRES, CWRES, NPDE, EWRES, etc. Consider example ..\example\fixedetas.ctl:

$PK
...
IF(COMACT==1) THEN
PREDCL=CL
PREDV=V
ENDIF
...
$LEVEL
SID=(3[1],4[2])
CID=(5[3],6[4])

$EST METHOD=ITS INTERACTION PRINT=1 NSIG=3 NITER=10 SIGL=6 FNLETA=0 MCETA=3
         LEVCENTER=1

$COV MATRIX=R UNCONDITIONAL
$TABLE ID PREDCL PREDV DOSE RATE TIME CONC IPRED PRED IRES RES WRES CWRES
 EWRES NPDE NOAPPEND ONEHEADER FILE=FIXEDETAS.tab NOPRINT
     ESAMPLE=1000 SEED=1115678

$TABLE ID PREDCL PREDV DOSE RATE TIME CONC IPRED PRED IRES RES WRES CWRES
 EWRES NPDE NOAPPEND ONEHEADER FIXEDETAS=(3-6) FILE=FIXEDETAS2.tab NOPRINT
 ESAMPLE=1000 SEED 1115678

Notice that FNLETA=0 so that the super-individual level etas (3, 4, 5, and 6) evaluated during the estimation remain in memory (and not replaced by the FNLETA step), and so they are available for the table step. Notice also that the first table produces the standard population diagnostics, but the second table will produce those for which only Etas 1 and 2 (the inter-individual etas) are set to 0. Notice also that PREDCL and PREDV, because they are evaluated conditionally when COMACT=1, will be the population values (with their respective level interpretations for each table) of CL and V. If you apply VARCALC=3 (see Section 64.4 above), then you will obtain standard errors that incorporate errors associated with the super-individual level etas that are in PREDCL and PREDV.

Brendel, Karl, Emmanuelle Comets, Celine Laffont, Christian Laveille, and France Mentré. 2006. “Metrics for External Model Evaluation with an Application to the Population Pharmacokinetics of Gliclazide.” Pharmaceutical Research 23 (9): 2036–49. https://doi.org/10.1007/s11095-006-9067-5.

Comets, Emmanuelle, Karl Brendel, and France Mentré. 2008. “Computing Normalized Prediction Distribution Errors to Evaluate Nonlinear Mixed Effects Models: The npde Add-on Package for R.” Computer Methods and Programs in Biomedicine 90 (2): 154–66. https://doi.org/10.1016/j.cmpb.2007.12.002.

Hooker, Andrew C., Christine E. Staatz, and Mats O. Karlsson. 2007. “Conditional Weighted Residuals (CWRES): A Model Diagnostic for the FOCE Method.” Pharmaceutical Research 24: 2187–97. https://doi.org/10.1007/s11095-007-9361-x.

Nguyen, Thu Thuy, Emmanuelle Comets, and France Mentré. 2012. “Extension of NPDE for Evaluation of Nonlinear Mixed Effect Models in Presence of Data Below the Quantification Limit with Applications to HIV Dynamic Model.” Journal of Pharmacokinetics and Pharmacodynamics 39: 499–518. https://doi.org/10.1007/s10928-012-9264-2.

Press, William H., Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. 1992. Numerical Recipes: The Art of Scientific Computing. 2nd ed. Cambridge University Press.