59 NONMEM Utility Routines

NONMEM-internal utility routines.

1 ANEAL

MEANING: CONSTRAINT subroutine

CONTEXT: NONMEM utility routine

USAGE:

$SUBROUTINES ADVAN3 TRANS4 OTHER = aneal.f90

DISCUSSION:

Additional algorithmic constraints may be imposed upon the model parameters, by use of the subroutine CONSTRAINT. This feature is available only for the EM and Bayesian algorithms.

One use would be to slow the rate of reduction of the diagonal elements of the OMEGA values during the burn-in phase of the SAEM method. This is shown in example 9, where a user supplied annealing algorithm is used to replace the built-in one described earlier. By specifying OTHER=ANEAL.f90, where ANEAL.f90 was originally derived from a template of CONSTRAINT.f90 in the source directory, the user supplied CONSTRAINT subroutine can be incorporated into the model.

In example 9, whenever the argument iteration number (ITER_NO) of CONSTRAINT changes, a new OMEGA is evaluated that is larger than what was determined by the SAEM update. Typically, this expansion algorithm should be such that its impact decreases with each iteration. This example is given in example 9.

See INTRODUCTION TO NONMEM 7, $ANNEAL to facilitate EM search methods.

REFERENCES: Guide Introduction_7

2 CELS

MEANING: CELS subroutine

CONTEXT: NONMEM utility routine

This feature is not fully documented. The interested user may be able to obtain more information by studying the appropriate sections of NONMEM code and previous examples that may be available from advanced users.

DISCUSSION:

CELS is a NONMEM utility routine. It may be called by a user-supplied CCONTR routine. It computes the default contribution to the conditional objective function from a L2 record.

Certain NONMEM error messages refer to CELS, for example:

PROGRAM TERMINATED BY OBJ, ERROR IN CELS
WITH INDIVIDUAL 1 (IN INDIVIDUAL RECORD ORDERING)
INTRAINDIVIDUAL VARIANCE OF DATA FROM OBS RECORD 1 ESTIMATED TO BE 0

A common situation giving rise to this message is the use of a proportional intra-individual error model while some predicted values for actual observations are zero or close to zero. (For example, if the first dose is an infusion and there is a “base-line” observation at the start of the infusion, the predicted level will be zero.)

CELS must not be called when a user-supplied CRIT routine is used.

REFERENCES: None.

3 CHOL

MEANING: CHOL subroutine

CONTEXT: NONMEM utility routine

USAGE:

USE SIZES, ONLY:: ISIZE,DPSIZE
INTEGER(KIND=ISIZE) :: MODE,N,IER
REAL(KIND=DPSIZE) :: R
DIMENSION R(5)
CALL CHOL (MODE,N,R,IER)

DISCUSSION:

CHOL is a NONMEM utility routine that may be called by any user routine. It computes the Cholesky square root of a given positive definite matrix A, i.e., the lower triangular matrix B satisfying $BB' = A$.

OPTIONS:

MODE (input): An integer variable or integer constant.

0 - the required initialization call.
1 - a normal call.

N (input): An integer variable or integer constant between 1 and 5; the dimension of the matrix.
R (input/output): An array of dimension at least N. On input: the ith row of the matrix A, at the i+1 st. call to CHOL. Elements R(j) with j>i are ignored. On output: the ith row of B, upon return from the i+1 st. call to CHOL.
IER (input/output): Upon input when MODE is 0:

0 - CHOL should terminate NONMEM execution when the matrix is algorithmically not positive definite.
1 - CHOL should always return control to PRED.

On output:

0 - Normal return.
1 - matrix is algorithmically not positive definite. (This can only happen if IER was 1 at the call with MODE=0.)

If the dimension of A is NxN, then CHOL is called N+1 times. The first call initializes the routine. On the i+1 st. call, the ith row is passed to CHOL, and CHOL returns the ith row of B (i=1,…, N).

REFERENCES: Guide II, section D.4, D.4.2; Guide II, Figure 20

4 CONSTRAINT

MEANING: CONSTRAINT subroutine

CONTEXT: NONMEM utility routine

DISCUSSION:

CONSTRAINT is a NONMEM utility routine. Additional algorithmic constraints may be imposed upon model parameters by use of the subroutine CONSTRAINT. This feature is available only for the EM and Bayesian algorithms.

The default CONSTRAINT.f90 is identical to source\CONSTRAINT.f90.

Option CONSTRAIN of the $ESTIMATION record and the $ANNEAL record may be used to give information to the subroutine. It may be modified by the user.

See INTRODUCTION TO NONMEM 7, $ANNEAL to facilitate EM search methods.

See INTRODUCTION TO NONMEM 7, Imposing Thetas, Omegas, and Sigmas by Algebraic Relationships: Simulated Annealing Example.

REFERENCES: Guide Introduction_7

5 ELS

MEANING: ELS subroutine

CONTEXT: NONMEM utility routine

USAGE:

USE SIZES,        ONLY:: ISIZE,DPSIZE,NO

...

INTEGER(KIND=ISIZE) :: IER1,IER2

REAL(KIND=DPSIZE) :: CNT,D

DIMENSION :: D(NO,2)

...

CALL ELS (CNT,D,IER1,IER2)

DISCUSSION:

ELS is a NONMEM utility routine. It computes the Extended Least Squares contribution to the objective function for the data from a single individual (“L1”) record. It may be called by CONTR or SPTWO.

OPTIONS:

CNT (output)

Contribution to -2log likelihood for data from the L1 record as computed by extended least squares.

D(J,1) (output)

Value of RES for Jth observation record of the L1 record when called by SPTWO; otherwise 0.

D(J,2) (output)

Value of WRES for Jth observation record of the L1 record when called by SPTWO; otherwise 0.

IER1 (output)

=0 - Normal return.

>1 - Abnormal return.

IER2 (output): Undocumented error return code.

ELS obtains all its inputs from NONMEM read-only variables.

Certain NONMEM error messages refer to ELS, for example:

PROGRAM TERMINATED BY OBJ, ERROR IN ELS
 WITH INDIVIDUAL 1 (IN INDIVIDUAL RECORD ORDERING)
 VAR-COV OF DATA FROM INDIVIDUAL RECORD ESTIMATED TO BE SINGULAR

REFERENCES: Guide I, section C.1; Guide V, section 5.2.2

6 FILES

MEANING: FILES Subroutine

CONTEXT: NONMEM utility routine

USAGE:

USE SIZES, ONLY: ISIZE
INTEGER(KIND=ISIZE) :: IUNIT
CALL FILES (IUNIT)

DISCUSSION:

Whenever a FORTRAN OPEN statement or a CLOSE statement is executed in a user-supplied routine, NONMEM should be informed. This should be done as follows: Immediately after the OPEN or CLOSE statement is executed, a call should be issued to FILES.

Input argument:

IUNIT: An integer variable or integer constant; the number of the logical unit involved in the OPEN or CLOSE statement.

If reads or writes are issued to a given unit, but no OPEN statement for this unit is executed, then when ICALL is 0, a call to FILES should be issued (before any I/O statement involving this file is executed), where IUNIT is the number of the logical unit. However, reads and writes to a unit which has not been opened (connected) may not be allowed by some operating systems.

REFERENCES: Guide III, section V.10

7 FLU

MEANING: FLU subroutine

CONTEXT: NONMEM user-installation routine

USAGE:

SUBROUTINE FLU (I)
RETURN
END

DISCUSSION:

Routine FLU is a user-installation routine. It allows the user to cause NONMEM to execute statements needed with certain operating systems.

FLU is called by NONMEM whenever the buffer for an I/O unit should be emptied. Buffer emptying is controlled in an operating system specific way, and often a compiler that is supported under an operating system will support a FORTRAN language extension that allows a buffer to be emptied. Whenever possible, such a FORTRAN instruction should be coded into FLU. FLU has one argument, I, which is the I/O logical unit number. The statement that should be inserted prior to RETURN follows, for two different operating systems, is as follows:

(Sun UNIX):

CALL FLUSH (I)

(MS WINDOWS 95):

CALL COMMITQQ (I)

REFERENCES: None.

8 GAMLN

MEANING: GAMLN function

CONTEXT: NONMEM utility routine

USAGE:

USE SIZES, ONLY: DPSIZE
REAL(KIND=DPSIZE) :: GAMLN,X
REAL(KIND=DPSIZE) :: LN_OF_GAMMA

 . ..
LN_OF_GAMMA= GAMLN(X)

DISCUSSION:

Gives the natural logarithm of the gamma function. It is more accurate than the Stirling’s approximation.

May be used in abbreviated code. May be used in the computation of the objective function because it is continuous. The first derivative (GAMLND1) and second derivative (GAMLND2) are also evaluated, and are continuous for X>0.0.

Input argument:

X: (input argument)

EXAMPLE OF USAGE:

      FAC=EXP(GAMLN(X+1.))

If X takes only integer values, then

      FAC=X!=X*(X-1)*(X-2)... *1

The GAMMA function in effect “fills in” the values of the factorial (but with a shift of 1 in the argument) for non-integer values of x.

If the argument is a random variable, it should always be positive. Eta partial derivatives are computed.

REFERENCES: None.

9 GETETA

MEANING: GETETA subroutine

CONTEXT: NONMEM utility routine

USAGE:

DIMENSION ETA(LVR)
DOUBLE PRECISION ETA
CALL GETETA (ETA)

DISCUSSION:

The NONMEM utility routine GETETA is used to obtain subject-specific values of eta. It may be called only when the data are population data.

Output argument:

ETA: An array into which GETETA stores subject-specific eta values ETA(1), ETA(2), etc. If a conditional estimation method is not used, GETETA stores zeros in the ETA array.

GETETA must always be initialized at ICALL=1. This simply involves calling GETETA at ICALL=1. No values are stored in ETA at this call.

When the data are population and abbreviated code is used, values of etas are always obtained by a call to GETETA occurring in generated or Library PK or PRED subroutines. (In the case of PK, the array of etas is stored in a common and is available to the ERROR subroutine as well.) This call occurs even if a conditional estimation method is not implemented, in which case GETETA stores zeros in ETA. As a result, when the first-order method is used, apparent subject-specific values computed in NM-TRAN abbreviated code are the required typical values. Thus, an executable resulting from using an abbreviated code for PRED (for PK and ERROR if PREDPP is used) can be used with any estimation method.

EXAMPLE OF USAGE:

IF (ICALL == 1) THEN
 ...
 CALL GETETA (ETA)
 ...
 RETURN
ENDIF

...
IF (ICALL == 2) CALL GETETA (ETA)

It is also possible to obtain simulated subject-specific values of eta. (See Section 59.19).

REFERENCES: Guide IV, section V.C.5

REFERENCES: Guide VI, section III.E.2, IV.B.2

10 MULT

MEANING: MULT subroutine

CONTEXT: NONMEM utility routine

USAGE:

USE SIZES, ONLY: DPSIZE
REAL(KIND=DPSIZE) :: A,B
CALL MULT (A,B)

DISCUSSION:

MULT is a NONMEM utility routine that may be called by PRED. It allows reweighting of the objective function. That is, it allows terms A(i) and B(i) of the contribution to the (default) objective function of the ith individual to be set to values other than the defaults, which are 1. A(i) multiplies the logarithm of the determinant of the variance - covariance matrix for the ith individual and B(i) multiplies the sum of the squared weighted residuals for the ith individual.

Input argument:

A: A variable or constant.
B: A variable or constant.

To set A(i) and B(i) to values a and b for the ith individual, PRED should issue a single call to MULT during any one or more calls to it from NONMEM where ICALL is 2 and some data record from the ith individual is in DATREC.

At these calls, set A=a and B=b.

At least one such call to MULT must be issued for each individual. If more than one such call is issued for the same individual, then A and B should have the same values across all such calls.

MULT can be used whenever ELS is called. MULT can also be used whenever NCONTR is called and NCONTR calls CELS. NCONTR calls CELS whenever CCONTR is not called, or when CCONTR is called and CCONTR only calls CELS.

REFERENCES: Guide II, section D.4, D.4.3

REFERENCES: Guide II, Figure 20

11 NCONTR

MEANING: NCONTR subroutine

CONTEXT: NONMEM utility routine

USAGE:

USE SIZES, ONLY: ISIZE,DPSIZE
INTEGER(KIND=ISIZE) :: IER1,IER2,L2R
REAL(KIND=DPSIZE) :: CNT

 ...
CALL NCONTR (CNT,IER1,IER2,L2R)

DISCUSSION:

NCONTR is a NONMEM utility routine. It computes the normal-based contribution to the objective function from the data from a single individual (“L1”) record.

NCONTR may be called by a user-supplied CONTR. NCONTR can be used when eta is from a normal distribution with mean 0 and variance-covariance OMEGA. It can be used with population data. It should not be used with single-subject data.

Output argument:

CNT: Contribution to -2log likelihood for data from the L1 record.
IER1: IER1 return codes:

1: HESSIAN OF OBJ. FUNC. FOR COMPUTING SOME CONDITIONAL ESTIMATE IS SINGULAR, IF LAPLACIAN IS NOT USED, OR NONPOSITIVE DEFINITE, IF LAPLACIAN IS USED.
2: ETA-DERIVATIVE OF SOME CONDITIONAL -2LOG LIKELIHOOD VALUE IS TOO LARGE.
3: NUMERICAL HESSIAN OF OBJ. FUNC. FOR COMPUTING SOME CONDITIONAL ESTIMATE IS NON POSITIVE DEFINITE.
11: IF CCONTR IS NOT USED AND F = PREDICTION, SOME INTRAINDIVIDUAL VARIANCE-COVARIANCE IS SINGULAR. IF CCONTR IS NOT USED AND F = LIKELIHOOD, SOME CONDITIONAL LIKELIHOOD VALUE IS NEGLIGIBLE.
12: IF CCONTR IS NOT USED AND F = PREDICTION, WEIGHTED SUM OF “SQUARED” INDIVIDUAL RESIDUALS IS INFINITE.
10+: IF CCONTR IS USED, 10 IS ADDED TO A NONZERO VALUE OF IER1 FROM CCONTR.

IER2

0 if error-recovery is to be implemented when IER1 is nonzero.

1 if NONMEM is to stop when IER1 is nonzero.

L2R: Index of L2 record giving rise to error condition.

Other inputs: NCONTR obtains all its inputs from NONMEM read-only global variables.

REFERENCES: None.

12 NWPRI

MEANING: NWPRI subroutine

CONTEXT: NONMEM utility routine

USAGE:

USE SIZES, ONLY: ISIZE,DPSIZE
INTEGER (KIND=ISIZE) :: NTHETA,NETA,NEPS,NTHP,NETP,NEPP,NPEXP
INTEGER (KIND=ISIZE) :: `ITYP`,`NSAM`,ISS
REAL (KIND=DPSIZE), ONLY: CNT
REAL (KIND=DPSIZE), ONLY: `PLEV`
INTEGER (KIND=ISIZE) :: `ICALL`
CALL `NWPRI` (NTHETA,NETA,NEPS,NTHP,NETP,NEPP,NPEXP,`ITYP`,`PLEV`, &
                     `NSAM`,ISS,CNT)

DISCUSSION:

The user-written PRIOR subroutine allows a penalty function based on a prior to be specified and added to the -2log likelihood function (Gisleskog et al, JPP, 2002, p. 473-505). This serves as a constraint on THETA, OMEGA, and SIGMA estimates and thus as a way for stable estimates of these parameters to be obtained with insufficient data.

NWPRI may be called by PRIOR. (See Section 57.6). It computes a function based on a prior that has a multivariate normal form for THETA, and also, in the case of population data, an inverse Wishart form for OMEGA (independent from the normal for THETA). The parameters of these forms are called “hyperparameters”, and the modeler fixes these to values of his/her choice.

Using NWPRI, several penalties functions, each with a different set of values for the hyperparameters, may be used simultaneously. Each set of values may be thought to arise from a different prior experiment (study). In actuality, there is a way to combine the several penalties functions into one (See above reference), and this is what NWPRI computes.

When NWPRI is used during a Simulation Step, it produces a random value of THETA and a random value of OMEGA from the prior. (See Simulation example). (See nwpri example). If NWPRI is used at ICALL=2, it need not be used at ICALL=4, and vice-versa.

NWPRI should always be called at ICALL=0 or ICALL=1.

Use only NWPRI for the new $ESTIMATION methods of NONMEM 7.

Arguments

Input arguments:

NTHETA, NETA, NEPS

NTHETA = number of Thetas to be estimated.

NETA = number of Etas to be estimated.
NEPS = number of epsilons to be estimated.

These are the dimensions of the THETA, OMEGA, and SIGMA arrays of the parameters that enter into the model for the data.

Before NONMEM 7.3, NEPS was ignored. With odd-type population data or with non-odd-type single-subject data, where OMEGA takes the place of SIGMA, the input argument NETA must be set.

NTHP, NETP, NEPP

NTHP = number of thetas which have a prior.

NETP = number of Omegas with prior.
NEPP = Number of Sigmas with prior.

The prior will only affect the initial subvector of THETA of dimension NTHP and the initial submatrix of OMEGA of dimension NETP (i.e. the submatrix consisting of the intersection of the first NETP rows and the first NETP columns of OMEGA). The initial subvector and submatrix are called the prior-affected parts of THETA and OMEGA. If NTHP=0 (NETP=0), the prior does not affect THETA (OMEGA) at all. During the Simulation Step, simulated values for the affected parts of THETA and OMEGA are obtained according to the prior distribution, and “the simulated values” for the unaffected parts of THETA and OMEGA are simply taken to be the values given in the NONMEM control stream or input Model Specification record. Before NONMEM 7.3, NEPP was ignored.

NPEXP: The number of prior experiments.
ITYP: Relevant only if NWPRI is called during a Simulation Step. Values are:

0: The value of the THETA vector is obtained from simple random sampling.
1: Within the given problem, NWPRI is to be called a specified number (NSAM) of times to obtain this number of different values of the THETA vector. These values are obtained by generating a Latin sample of size NSAM from equiprobable partitions of an ellipsoid in THETA space (hyper-ellipsoidal sampling), followed by sampling a point “uniformly” from each partition. This scheme may be used, for example, when the problem has NSAM subproblems, in which case, NWPRI would be called NSAM times, once each time during the problem when ICALL=4, and at each of these calls, a different random value of THETA will be produced.
2: Just as with value 1, but the NSAM values are obtained by generating a Latin sample of size NSAM from equiprobable partitions of an ellipsoid in THETA space (hyper-ellipsoidal sampling), followed by taking the “center point” of each partition.

In all three cases, each new value of the OMEGA matrix is obtained from simple random sampling.

After each call to NWPRI, the simulated values for THETA, OMEGA and SIGMA may be found in global variables and thus they are communicated directly to NONMEM. (See Section 53.87).

PLEV: When NWPRI is being used at ICALL=0 or 1, but NWPRI will not be used at ICALL=4 (i.e. during the Simulation Step), PLEV can be set to 0. When it is being used at ICALL=2, PLEV can also be set to 0. When NWPRI is being used at ICALL=0 or 1 and it will be used at ICALL=4, or when it is being used at ICALL=4, then PLEV must be set to a fraction strictly less than 1, e.g. 0.999. PLEV is double precision with NONMEM 7, and is single precision with NONMEM VI.

A value of THETA will actually be obtained using a truncated multivariate normal distribution, i.e. from an ellipsoidal region R1 over which only a fraction of mass of the normal occurs. This fraction is given by PLEV. The distribution is further truncated to R2, the subregion of R1 lying within the rectangular boundaries defined on the $THETA record. Simple random sampling occurs in R2. Latin sampled partitions are partitions of R1. However, when ITYP=1, if a uniformly sampled point from a partition lies outside R2, it is replaced by a point obtained by simple random sampling from R2. When ITYP=2, if the center point of a partition lies outside R2, an abort occurs.

NSAM: Relevant only if NWPRI is called during a Simulation Step. Consider two cases.

Latin hyper-ellipsoid sampling is used with ITYP=1, and b) simple random sampling along with the adjustment for small sample correlation effect is used (See next input argument). In case a) NSAM should equal the exact total number of different values of THETA that must eventually be produced over the entire NONMEM problem. In case b) NSAM should be no less than this number.

ISS: Relevant only if NWPRI is called during a Simulation Step. A THETA value is obtained by transforming a value from the standard multivariate normal distribution - called here “the standard value”. The correlation matrix of the standard normal is the identity matrix. When NSAM is small, the estimated correlation matrix from the sampled standard values might not be quite close to the identity matrix - this is here called “the small sample correlation effect”.

1: An adjustment is made for the small sample correlation effect, by first transforming the NSAM standard values altogether into new values which are very nearly standard multivariate normal values and such that the sample correlation matrix of these new values is exactly the identity matrix.
0: No adjustment is made for the small sample correlation effect.

Output argument:

CNT: Relevant only if NWPRI is called at ICALL=2. CNT is the penalty.

Use of $PRIOR

A $PRIOR record may be used instead of a user-written PRIOR subroutine, in which case the input arguments listed above may be specified as options of $PRIOR, and the value of CNT may be displayed in the NONMEM output. (See Section 57.6).

Use of $THETA, $OMEGA, $SIGMA records

Specifying the multivariate normal form for THETA (NTHP>0):

After the usual set of $THETA records, add a second set of $THETA records giving the mean of the multivariate normal. The values on these records must be fixed. They should number NTHP altogether.

After the usual set of $OMEGA records, add a second set of $OMEGA records giving the variance-covariance matrix of the multivariate normal. (If there are no etas in the model, there will not be a usual set of $OMEGA records.) The values on these records must be fixed. They should number $(NTHP+1) \times NTHP / 2$ altogether, including implicit 0’s, which may occur because the variance-covariance matrix of the normal may be block-diagonal. If a (regular) theta is fixed, the corresponding values of the mean and variance-covariance matrix of the normal are ignored.

Specifying the inverse Wishart form for OMEGA (NETAP>0; population data only):

After the second set of $OMEGA records (or if NTHP=0, after the first set of $OMEGA records), add another set of $OMEGA records giving the mode of the inverse Wishart. The values on these records must be fixed. They should number $(NETP+1) \times NETP / 2$ altogether, including implicit 0’s, which may occur because the mode of the inverse Wishart may be block-diagonal. If the prior-affected part of OMEGA is given as a block-diagonal matrix, then the mode must conform to this structure. The SAME attribute can be used.

With each diagonal block of (the prior-affected part of) OMEGA, there corresponds a number of “degrees of freedom” of the inverse Wishart. All blocks of a given block set are constrained to be equal (by using the SAME attribute), and therefore, to each of these blocks there corresponds the same number of degrees of freedom. After the second set of $THETA records (or if NTHP=0, after the first set of $THETA records), add another set of $THETA records, giving for each block set in turn the number of degrees of freedom for the blocks of the set. The values on these records must be fixed. There should be as many values as there are block sets with the prior-affected part of OMEGA.

The inverse Wishart for a given block of OMEGA may be explicitly given as “perfectly flat” by specifying the number of degrees of freedom to be the negative of the dimension of the block, minus 1. In this case the mode for the block may be simply taken to be the 0 matrix (or any positive definite matrix). If a block is fixed, the corresponding values of the mode and number of degrees of freedom are ignored.

With NONMEM 7.3 and higher, you can use more informative record names that obviate the need for specific ordering of the additional $THETA and $OMEGA records used by NWPRI, and provide an alternative source for values of arguments NTHETA, NETA, et al.

$THETAP for theta priors
$THETAPV for variance-covariance matrix for theta’s
$OMEGAP for OMEGA prior
$OMEGAPD for degrees of freedom (or dispersion factor) for omega prior
$SIGMAP for SIGMA prior
$SIGMAPD for degrees of freedom (or dispersion factor) for sigma prior

If the informative record names are used, the records may be in any order, and the options of $PRIOR need not be specified. Note that the name of the record describes the kind of information it gives, rather than the structure of the information. E.g., in the example below, $THETAPV is used instead of an $OMEGA record and $OMEGAPV is used instead of a $THETA record. (See Section 50.49, Section 50.29, Section 50.39)

Examples

Here are three examples of extra $THETA and $OMEGA records specifying prior information from an experiment. This information concerns all the regular elements of THETA and OMEGA. All examples are shown with a $PRIOR record although a PRIOR subroutine could be used instead.

$PRIOR NWPRI NTHETA = 3 NETA = 3 NTHP = 3 NETP = 3 NPEXP = 1
...
$`THETA` 3 FIX .08 FIX .04 FIX
$THETA 100 FIX
$`OMEGA` BLOCK (3) .494 .00207 .0000847 .000692 .0000471 .0000292 FIX
`$OMEGA` BLOCK (3) .7 .04 .05 .02 .06 .08 FIX

Perhaps it might be more perspicuous to organize the prior information thusly:

`$PRIOR` `NWPRI` NTHETA = 3 NETA = 3 NTHP = 3 NETP = 3 NPEXP = 1
...

;prior information for `THETA`
$THETA 3 FIX .08 FIX .04 FIX
$`OMEGA` BLOCK (3) .494 .00207 .0000847 .000692 .0000471 .0000292 FIX

;prior information for `OMEGA`
$OMEGA BLOCK (3) .7 .04 .05 .02 .06 .08 FIX
$`THETA` 100 FIX

With NONMEM 7.3 and higher, informative record names can be used:

`$PRIOR` `NWPRI`
...

;prior information for `THETA`
$THETAP 3 FIX .08 FIX .04 FIX
$THETAPV BLOCK (3) .494 .00207 .0000847 .000692 .0000471 .0000292 FIX

;prior information for `OMEGA`
$OMEGAP BLOCK (3) .7 .04 .05 .02 .06 .08 FIX
$OMEGAPD 100 FIX

Multiple Experiments

There may be prior information from a number of experiments, in which case the $THETA and $OMEGA records specifying this information for each experiment may be stacked, e.g.

;usual records
$THETA (0,4,10) (0,.09,.5) (.004,.01,.9)
$`OMEGA` BLOCK (3) .7 .04 .05 .02 .06 .08

;prior information from experiment 1
$THETA 3 FIX .08 FIX .04 FIX
$`THETA` 100 FIX
$OMEGA BLOCK (3) .494 .00207 .0000847 .000692 .0000471 .0000292 FIX
$`OMEGA` BLOCK (3) .7 .04 .05 .02 .06 .08 FIX

;prior information from experiment 2
$THETA 2 FIX .05 FIX .04 FIX
$`THETA` 50 FIX
$OMEGA BLOCK (3) .6 .003 .0001 .0004 .00001 .00003 FIX
$`OMEGA` BLOCK (3) .9 .02 .05 .01 .06 .09 FIX

Blocking on the variance-covariance matrix of the normal form for THETA need not be the same across experiments.

Limitation: There must be at least one THETA parameter and one OMEGA parameter in the model.

REFERENCES: None

13 PASS

MEANING: PASS subroutine

CONTEXT: NONMEM utility routine

USAGE:

USE SIZES, ONLY: ISIZE
INTEGER(KIND=ISIZE) :: MODE
CALL PASS (MODE)

DISCUSSION:

The NONMEM utility routine PASS can be used to read and/or to modify (“transgenerate”) data records. This can be done at both the beginning and ending of a problem.

When ICALL=0, 1, or 3, PASS may be called by PRED, or by a subroutine of PRED (such as PREDPP’s INFN), or PASS may be called from $PRED or $INFN abbreviated code. With abbreviated code, the variable MODE should be used as the argument to PASS and set as described below.

Input argument:

MODE (input)

MODE=0: Initialize PASS.

MODE=1: store next data record in DATREC (and in PASSRC of module NMPRD9); any modifications of this record performed by PRED are ignored by NONMEM.
MODE=2: store next data record in DATREC (and in PASSRC of module NMPRD9); any modifications of this record performed by PRED are placed into NONMEM’s internal data set.
MODE=3: terminate PASS (can be used only after MODE=1).

Output argument:

MODE (output): Upon return from a call to PASS with MODE=1 or 2, PASS sets MODE=0 when there are no more data records to be passed.

The first call to PASS with MODE=1 or 2 produces the first data record, the second call produces the second data record, etc. A value of NEWIND appropriate for the data record currently in DATREC and PASSRC is also stored by PASS in NWIND of module NMPRD_INT.

The NONMEM utility routine GETETA can be used to obtain subject-specific estimates of the eta’s. The values returned for the eta’s with each call to GETETA are appropriate for the individual whose data record is currently in DATREC. When PASS is used in abbreviated code, GETETA is used automatically. Any PRED-defined items that are displayed are available in module NMPRD4 with each data record.

The value of MIXEST is appropriate for the individual whose data record is currently in DATREC.

Since the call with ICALL=3 actually occurs before the Table and Scatterplot Steps, new data items can be tabled and scatterplotted.

The user should not use PASS to modify either the ID or MDV data items.

Within an initialization or finalization block, code to test MODE and call PASS is generated by NM-TRAN when abbreviated code such as the following is present:

IF (`ICALL` == 0) THEN
 DO WHILE (DATA)
 ... transgeneration statements ...
 ENDDO
 RETURN
ENDIF

(See Infn_pass_interpolation_example, Infn_example). (See Section 53.46).

REFERENCES: Guide II, section D.2.2; Guide VI, section VI.A, Figure 37

14 PHI

MEANING: PHI function

CONTEXT: NONMEM utility routine

USAGE:

USE SIZES, ONLY: DPSIZE
REAL(KIND=DPSIZE) :: PHI,X,CUMDIS
CUMDIS= PHI(X)

DISCUSSION:

Gives the value of the cumulative distribution function.

Input argument:

X: EXAMPLE OF USAGE:

PHI may be used in abbreviated code.

A = `THETA`(1) * EXP(ETA(1))
B = PHI(A)

The argument may be a random variable, but no eta partial derivatives are computed.

REFERENCES: None.

15 RANDMT, RANDMTU

MEANING: RANDMT and RANDMTU functions

CONTEXT: NONMEM utility routine

BACKGROUND:

The RANDOM subroutine may be used with all versions of NONMEM. The statement CALL RANDOM(K,R) obtains a random number R from random source K. It may be used in the Simulation Step (ICALL=4). The seed is described on the $SIMULATION record for the Kth source. It may also be called when data averages are being computed (ICALL=5). RANDOM can compute only from the NORMAL or UNIFORM distributions.

DISCUSSION:

Functions RANDMT and RANDMTU may be used with NONMEM 7.5 and later. They can be used to obtain random numbers during the Estimation or other steps, with seeds other than the ones used by NONMEM, or with other distributions.

Function RANDMT is simliar to the RANDOM subroutine, and allows the Student T distribution to be used in addition to normal and uniform distributions.

Function RANDMTU is independent of the RANDOM subroutine. It may be used with all NONMEM steps. The user has complete control and may specify starting seed values in the abbreviated code.

USAGE:

These two USE statements are required.

USE NM_INTERFACE,ONLY: RANDMT
USE NM_INTERFACE,ONLY: RANDMTU

They can be obtained with the statement include nonmem_reserved_general.

USAGE of RANDMT:

RVAL1=RANDMT(N)

Where

N: integer

0: get uniform random variable
N=1 or N>101: get normal random variable
N>1 and N<=100: get Student-t random variable with N-1 degrees of freedom

For more details about RANDMT, See INTRODUCTION TO NONMEM 7: RANDMT Function

USAGE of RANDMTU:

RVAL2=RANDMTU(N,SOURCE,STARTSEED,RANM)

Where

N: integer (same as with RANDMT)

0: get uniform random variable
N=1 or N>101: get normal random variable
N>1 and N<=100: get Student-t random variable with N-1 degrees of freedom

SOURCE: Integer between 0 and 20, allowing the user to maintain up to 20 distinct seed sources. IF SOURCE=0, it uses the default source.

STARTSEED: Integer, which if not 0, causes initialization with the designated starting seed. STARTSEED should be set to 0 after this initialization. If SOURCE=0, this argument is ignored. If STARTSEED is <0, then the starting seed will be 10000*(seconds after midnight)+ABS(STARTSEED). (This is similar to the CLOCKSEED option of the $SIMULATION record, which 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.) The actual startseed is stored in RANDMTU_STARTSEED(SOURCE), module NM_BAYES_INT, accessible by using nonmem_reserved_general (See example below).

RANM: Integer, indicating the type of random number generator algorithm to use, as follows:

0: ran0 of reference [5], minimal standard generator
1: ran1 of reference [5], Bays and Durham.
2: ran2 of reference [5].
3: ran3 of reference [5], Knuth.
4: NONMEM’s traditional random number generator, used as default in $SIMULATION

For more details bout RANDMTU, See INTRODUCTION TO NONMEM 7: RANDMTU Function

EXAMPLE:

In the following, 2 random sources are used, with different starting seeds

`$ERROR` ; may also be used in other records
...

include nonmem_reserved_general
...

; Make sure NMTRAN is aware of RVAL1 and RVAL2, so it is accessible
; for `$TABLE` output
RAVL1 = 0.0
RVAL2 = 0.0
RVAL3 = 0.0
IF(`NEWIND`==0) THEN
 ; Initialize only at beginning of the data set.
 ; 3=t-distribution with DF=2

 RVAL1 = RANDMTU(3,1,-1234,-1)
 write(*, *) RANDMTU_STARTSEED(1)
 RVAL2 = RANDMTU(3,2,5678,-1)
ENDIF
RVAL1 = RANDMTU(3,1,0,-1)
RVAL2 = RANDMTU(3,2,0,-1)
RVAL3 = RANDMT(3)

Note that calling RANDMT(N) (or equivalantly RANDMTU(N,0,...)) did not require an initialization call (and would be ignored if you did), as this is under the control of $EST, $SIML, or $TABLE.

REFERENCES: Guide Introduction_7

16 RANDOM

MEANING: RANDOM subroutine

CONTEXT: NONMEM utility routine

USAGE:

USE SIZES, ONLY: ISIZE,DPSIZE
INTEGER(KIND=ISIZE) :: KR
REAL(KIND=DPSIZE) :: R
CALL RANDOM (KR,R)

DISCUSSION:

The NONMEM utility routine RANDOM may be called (by PRED, PK or ERROR) during the Simulation Step (ICALL=4) and when data averages are being computed (ICALL=5) to obtain numbers from different random sources.

Input argument:

KR An integer variable or integer constant. The index of a random source.

Output argument:

R A number from the KRth source. R is double precision with NONMEM 7, and is single precision with NONMEM VI. Each time RANDOM is called, a new number is output.

Random sources are defined using the $SIMULATION record. Typically, the first source is reserved by NONMEM to generate realizations of the eta and epsilon variables and/or to randomly mix individuals into different subpopulations according to the mixing parameter. Thus, typically, RANDOM must be called with KR > 1. However, there are cases when NONMEM itself needs no random sources, and all defined random sources can be used by PRED (See references).

A random source that is specified on the $SIMULATION record as NONPARAMETRIC should not be used to obtain a value of R.

EXAMPLES OF USAGE:

WT (an item in the data record, and a local variable with the same name) is generated with a simulated value having mean 70 and standard deviation 7. The second random source is assumed to be pseudo-normal, which is the default distribution for each random source.

An abbreviated code may use RANDOM explicitly in either a simulation or a data average block. When RANDOM is called in abbreviated code, R is reserved and is used for the random number. E.g.,

IF (`ICALL` == 4 .AND. `NEWIND` /= 2) THEN
 CALL RANDOM(2,R)
ENDIF
IF (`ICALL` == 4) WT = 70 + 7 * R

REFERENCES: Guide IV, section III.B.13, IV.I REFERENCES: Guide V, section 12.4.8

17 SDE

MEANING: SDE subroutine

CONTEXT: NONMEM utility routine

Utility subroutine sde.f90 numerically evaluates SDE (stochastic differential equations).

USAGE:

$SUBROUTINE ADVAN6 TOL = 10 DP OTHER = sde.f90

See sde9.ctl, sde9b.ctl and sde.ctl. They, along with the data files, can be found in the examples directory of the NONMEM 7 distribution medium. See INTRODUCTION TO NONMEM 7, Repeated Observation Records See INTRODUCTION TO NONMEM 7, Stochastic Differential Equation Plug-In

REFERENCES: Guide Introduction_7

18 SIMEPS

MEANING: SIMEPS subroutine

CONTEXT: NONMEM utility routine

USAGE:

USE SIZES, ONLY: DPSIZE,LVR
REAL(KIND=DPSIZE):: EPS(LVR)
CALL SIMEPS (EPS)

DISCUSSION:

The NONMEM utility routine SIMEPS can be called by PRED during the simulation step, to obtain simulated epsilon values. It may be called only when ICALL=4.

Output argument:

EPS An array into which SIMEPS stores simulated epsilon values EPS(1), EPS(2), …. The dimension of the array may be smaller than the maximum, e.g., it may equal the number of epsilons in the problem.

Simulated epsilon values arise from a multivariate normal pseudo-random distribution with mean 0 and variance-covariance as specified for SIGMA. With different calls to SIMEPS with different level-two records, new and different simulated epsilon values are obtained. (A level-two record has the same value of L2 and may also be called an L2 record.) (See L2). (When the L2 data item is not defined, all data records are level-two records, and so different data records are different level-two records.) By default, with different calls to SIMEPS with the same level-two record, the same simulated epsilon values are obtained - those obtained at the first call with the record. (There is an advanced feature whereby records are “repeated” (See Section 53.94), and when records of a level-two record are being repeated, with different calls to SIMEPS with the same level-two record, the values obtained are the last values stored in these variables when the record was previously passed to PRED.)

If, though, the NEW option is used with the first random source on the $SIMULATION record, then each time SIMEPS is called (with the same or different level-two record), new and different values are obtained. With any particular call to SIMEPS, the effect of the NEW option can, though, be overridden; (See Section 53.98)

So that simple simulation can be easily implemented with abbreviated code, values of epsilon are obtained by calls to SIMEPS occurring in the generated subroutine. When the data are population data and the Simulation Step is implemented, SIMEPS is called once with every call to PRED (or once at every call to ERROR if PREDPP is used). These calls are implemented so that even if, initially, the Simulation Step is not implemented, the NONMEM executable resulting from using an abbreviated code for PRED (or for ERROR if PREDPP is used) can be reused with a run implementing the Simulation Step.

Additional calls to SIMEPS may appear in simulation blocks of $PRED and $ERROR abbreviated code. There is an analogous routine SIMETA. (See Section 59.19).

REFERENCES: Guide IV, section III.B.13 REFERENCES: Guide V, section 12.4.8 REFERENCES: Guide VI, section IV.B.1, IV.B.2, Figure 10

19 SIMETA

MEANING: SIMETA subroutine

CONTEXT: NONMEM utility routine

USAGE:

USE SIZES, ONLY: DPSIZE,LVR
REAL(KIND=DPSIZE):: ETA(LVR)
CALL SIMETA (ETA)

DISCUSSION:

The NONMEM utility routine SIMETA can be called by PRED during the simulation step, to obtain simulated eta values. It may be called only when ICALL=4.

Output argument:

ETA An array into which SIMETA stores simulated eta values ETA(1), ETA(2), … The dimension of the array may be smaller than the maximum, e.g., it may equal the number of etas in the problem.

Simulated eta values arise from a multivariate normal pseudo-random distribution with mean 0 and variance-covariance as specified for OMEGA. With different calls to SIMETA with different individual records, new and different simulated eta values are obtained. (With population data, an individual record has the same value of ID and may also be called a level-one “L1” record.) (See ID). By default, with different calls to SIMETA with the same individual record, the same simulated eta values are obtained - those obtained at the first call with the record. (There is an advanced feature whereby records are “repeated” (See Section 53.94), and when records of an individual record are being repeated, with different calls to SIMETA with the same individual record, the values obtained are the last values stored in these variables when the record was previously passed to PRED.)

If, though, the NEW option is used with the first random source on the $SIMULATION record, then each time SIMETA is called (with the same or different individual record), new and different values are obtained. Thus, for example, when PRED is called with the first data record of an individual record, PRED can in turn call SIMETA multiple times until values are obtained such that, for example, ETA(2) is not larger than 5 in absolute value; that is, values can be obtained from a truncated distribution (See below). With any particular call to SIMETA, the effect of the NEW option can, though, be overridden; (See Section 53.98)

So that simple simulation can be easily implemented with abbreviated code, values of etas are obtained by calls to SIMETA occurring in the generated subroutine. When the data are population data, SIMETA is called once per individual record by PRED (or PK if PREDPP is used). (In the case of PK, the array of etas is stored in a common and is available to the ERROR subroutine as well.) When the data are single-subject data, SIMETA is called once at every call to PRED (once at every call to ERROR if PREDPP is used).

These calls are implemented so that even if, initially, the Simulation Step is not implemented, the NONMEM executable resulting from using an abbreviated code for PRED (for PK or ERROR if PREDPP is used) can be reused with a run implementing the Simulation Step.

EXAMPLES OF USAGE:

In this example, the value of ETA(2) used by PRED will be less than 5 in absolute value. For this code to have the desired effect, the option NEW must be used in the $SIMULATION record.

IF (`ICALL` == 4 .AND. `NEWIND` /= 2) THEN
 DO WHILE (ABS(ETA(2)) > 5)
 CALL SIMETA (ETA)
 ENDDO

ENDIF

Suppose there are two etas to be selected in this manner. Each one needs its own CALL SIMETA loop, because each CALL SIMETA replaces all the etas.

IF (`ICALL` == 4 .AND. `NEWIND` /= 2) THEN
 DO WHILE (ABS(ETA(2)) > 5)
 CALL SIMETA (ETA)
 ENDDO
 ETA2 = ETA(2)
 DO WHILE (ABS(ETA(1)) > 0.52)
 CALL SIMETA (ETA)
 ENDDO
 ETA1 = ETA(1)

ENDIF

Another way this can be implemented is as follows:

IF (`ICALL` == 4 .AND. `NEWIND` /= 2) THEN
 DO WHILE (ABS(ETA(2)) > 5 .OR. ABS(ETA(1)) > 0.52)
 CALL SIMETA (ETA)
 ENDDO
 ETA2 = ETA(2)
 ETA1 = ETA(1)

ENDIF

There is an analogous routine SIMEPS. (See Section 59.18).

REFERENCES: Guide IV, section III.B.13, IV.I REFERENCES: Guide V, section 12.4.8 REFERENCES: Guide VI, section III.E.2, IV.B.2

20 SUPP

MEANING: SUPP subroutine

CONTEXT: NONMEM utility routine

USAGE:

USE SIZES, ONLY:: ISIZE
INTEGER(KIND=ISIZE):: IE,IC

 ...
CALL SUPP (IE,IC)

DISCUSSION:

The NONMEM utility routine SUPP is used to suppress portions of the NONMEM output report. SUPP may be called only when ICALL is 0, 1, or 3.

Input argument:

IE An integer variable or integer constant. If the value is 1, the output from the Estimation Step is not included in the output report.

IC An integer variable or integer constant. If the value is 1, the output from the Covariance Step is not included in the output report.

SUPP is typically used during simulation with numerous sub-problems. A user-supplied INFN routine, or similar code in a user-supplied PRED routine, presumably calculates and prints summary statistics based on the final parameter estimates, etc. Printed values from each sub-problem are of little interest and may be suppressed to reduce the size of the report file.

When IE=1, pages with these headings are suppressed:

MINIMUM VALUE OF OBJECTIVE FUNCTION
FINAL PARAMETER ESTIMATE

When IC=1, pages with these headings are suppressed:

STANDARD ERROR OF ESTIMATE
COVARIANCE MATRIX OF ESTIMATE
CORRELATION MATRIX OF ESTIMATE
INVERSE COVARIANCE MATRIX OF ESTIMATE

The option PRINT=0 of the $ESTIMATION record is used to suppress the Intermediate printout from the Estimation Step (“MONITORING OF SEARCH:…”)

EXAMPLES OF USAGE:

CALL SUPP(1,1)

This suppress both Estimation and Covariance Step output.

REFERENCES: None.

21 TNPRI

MEANING: TNPRI subroutine

CONTEXT: NONMEM utility routine

USAGE:

USE SIZES, ONLY: ISIZE,DPSIZE
INTEGER (KIND=ISIZE) :: `IFND`,MODE,`ITYP`,`NSAM`,ISS,IVAR
REAL (KIND=DPSIZE) :: `PLEV`
REAL (KIND=DPSIZE) :: CNT
REAL (KIND=ISIZE) :: `ICALL`
CALL `TNPRI` (`IFND`,MODE,`ITYP`,`PLEV`,`NSAM`,ISS,IVAR,CNT)

DISCUSSION:

The model parameters THETA, OMEGA and SIGMA parameters may be constrained in various ways. Those parameters whose values are not fixed are transformed in a one-to-one manner to parameters whose values are not at all constrained a priori (the “unconstrained parameters”, or UCP). The user-written PRIOR subroutine allows a penalty function based on a frequency prior to be specified and added to the $-2\log$ likelihood function (Gisleskog et al, JPP, 2002, p. 473-505). This function serves as a constraint on the estimates of THETA, OMEGA and SIGMA and thus as a way for stable estimates to be obtained with insufficient data. TNPRI may be called by PRIOR. (See Section 57.6). It computes a penalty function based on a frequency prior that has a multivariate normal form for all the UCP. The one-to-one transformation between THETA, OMEGA and SIGMA and the UCP induces a frequency prior for THETA, OMEGA and SIGMA itself. This latter form is called the “transformed normal” form.

Both forms have the same parameters, which are called “hyperparameters”. The values of the hyperparameters are produced with a prior NONMEM problem (called the “prior problem”) used to analyze a prior data set. Thus the frequency prior has an empirical nature. The prior problem should have implemented the Covariance Step (the Estimation Step could have been implemented in yet an earlier problem). Whenever a problem implements the Covariance Step and a model specification file has been output, this file automatically contains information - referred to as the “prior information”, and including values for the hyperparameters - which can be used when TNPRI is used in a subsequent problem.

The abilty to use TNPRI suggests that with any serious analysis undertaken in the future, though it may not be necessary to save the data that were analyzed, a model specification file should be produced and safely stored along with a control stream showing the particulars that gave rise to this file. There is a special way to input the prior information from the model specification file for its use with a subsequent application of TNPRI. (Alternatively, a NONMEM control stream may consist of multiple problems, where one problem (A) uses TNPRI, and the prior problem appears as an earlier problem (B) in the same control stream. The prior information from problem B can be available to problem A without the need to use a model specification file.)

When TNPRI is used during a Simulation Step, it produces a random value of the vector of all model parameters (whose values are not fixed in the parameter records) from the frequency prior. (See Simulation example). (See tnpri example). If TNPRI is used during a given problem at ICALL=2, it should not be used during the same problem at ICALL=4, and vice-versa.

TNPRI may be used at ICALL=0 or ICALL=1, in which case only the values of its input arguments are checked. If it is not used at ICALL=0 or ICALL=1, then the first time it is used at ICALL=2 or ICALL=4, checking will occur.

Do not use TNPRI with the new $ESTIMATION methods of NONMEM 7.

OPTIONS:

Input argument:

IFND: 0: Indicates that a prior problem has been specified either earlier in the current control stream or in a previously run control stream, and that this problem has output a model specification file. The prior information will have automatically been stored in the model specification file. The prior problem should have implemented the Covariance Step (which should have computed either the default variance-covariance matrix or two times the inverse of the R matrix). The current control stream should consist of at least two problems: the problem that uses TNPRI and a preceeding problem that serves only to input the prior information from the model specification file, such as the following:

$PROB READ THE MODEL SPECIFICATION FILE
$DATA data
$INPUT ID TIME AMT DV TYPE SS II
$`MSFI` msf1 ONLYREAD

This simple problem (B) should be the last problem in the current control stream that appears before the problem using TNPRI (problem (A)) and inputs a model specification file. The problem specification should not include any task records, but it may include e.g. a $SUBROUTINES record and/or abbreviated code, if problem B is the first problem in the control stream, and these elements will be needed for a subsequent problem. If problem A inputs a model specification file msf2, any prior information in msf2 will be ignored. If problem A outputs a model specification file msf3 and implements the Covariance Step, the prior information stored in msf3 is based not only on the data being analyzed, but also on the prior information in msf1.

1: Indicates that along with the current problem (with which TNPRI is being used), a prior problem has been specified earlier in the same control stream. The prior problem is taken to be the last problem in the control stream that appears before the current one and that implements the Covariance Step (which should have computed either the default variance-covariance matrix or two times the inverse of the R matrix).

MODE: Discussion: In any NONMEM problem specification, all of the last contiguously listed parameters on the $THETA record, whose values on this record are fixed, and assuming the value of the very last parameter itself is fixed, are called the “terminal fixed THETA parameters”. There may be no such parameters. E.g. if the $THETA record is

`$THETA` 3.1 FIX 6.3 .01 400 FIX 7 FIX 90 FIX

then THETA parameters 4-6 are the terminal fixed THETA parameters. If the record is

`$THETA` 3.1 FIX 6.3 .01 400 FIX 7 FIX 90

then there are no terminal fixed THETA parameters, because the value 90 of the very last parameter is not fixed. Similarly, all the last contiguous fixed parameters in OMEGA (SIGMA), assuming the values in the very last block set are fixed, are called the “terminal fixed OMEGA (SIGMA) parameters”. (Recall that in a diagonal $OMEGA record, all the values form separate block sets.)

There are some rules concerning the parameter records of a NONMEM control stream that relate the parameters of the prior problem to those of the current problem:

I. With the exception of the terminal fixed THETA parameters in the prior problem, which are altogether ignored in the current problem, the THETA parameters of the prior problem must be the first THETA parameters listed in the $THETA record of the current problem, whether or not their values are fixed in the prior problem.
1. The order of these first THETA parameters on the $THETA record must be the same as that in the $THETA record of the prior problem. The initial estimates of these parameters need not be the same, but any upper or lower bounds must be the same.
1. Among the first THETA parameters, values that are fixed on the $THETA record of the prior problem must also be fixed on the $THETA record of the current problem, although the values need not be the same, but not vice-versa; see discussion below concerning parameter interpretation.
1. THETA parameters peculiar to the current problem may be listed at the end of the $THETA record of the current problem, after the parameters that are shared between the prior and current problems.

E.g. if the $THETA record with the prior problem is

`$THETA` 3.1 FIX 6.3 (0,.01) 400 FIX 7 FIX 90 FIX

then the $THETA record with the current problem may be

`$THETA` 10 FIX 7.0 (0,.03) 80.2 100.8

in which case the values 80.2 and 100.8 are the values of the THETA parameters that are peculiar to the current problem. The 2nd and 3rd parameters are used in both prior and current problems, but have different values in the two problems. The first parameter may be used in the current problem, in which case its value is fixed to 10.

The same rules apply to the $OMEGA and $SIGMA records of the current problem. A model specification file from a previous instance of the current problem may be used, as usual, providing that if with the previous instance, parameter records were used, the above rules were followed.

Parameter Interpretation: The prior information concerns all parameters appearing in the prior problem whose values are not fixed in that problem. However, a parameter with a nonfixed value in the prior problem may have a fixed value in the current problem, and this gives rise to some ambiguity. Such a parameter may have an interpretation in the prior problem that remains unchanged in the current problem. Then it needs to be identified as a “shared parameter”, in which case, because the value of the parameter is fixed with the current model, the frequency prior needs to be adjusted in order to approximate the $-2\log$ likelihood function for the prior data when the parameter value is regarded as being fixed in the prior model. (or when TNPRI is being called in the Simulated Step, in order to use a conditional prior distribution, given the parameter value is the fixed value). A shared parameter may not actually be used in the current problem, but the salient point is that with the current problem, the interpretation of this parameter remains unchanged. Alternatively, the parameter’s interpretation may change in the current problem. Then it needs to be identified as a “prior-specific” parameter, in which case if it is actually used in the current model, its value will be the fixed value, but for the purpose of applying the prior information, TNPRI will regard this parameter as one specific to the prior model. The MODE argument concerns this distinction.

0: All parameters with a nonfixed value in the prior problem but a fixed value in the current problem are identified as being prior-specific parameters. However, the final estimates of these parameters are their fixed values.

1: All parameters with a nonfixed value in the prior problem but a fixed value in the current problem are identified as being prior-specific parameters. If the Estimation Step is implemented, the final estimate for such a parameter is the same one that would result from not fixing the value of the parameter in the current problem and letting it be estimated, assuming that the parameter is not at all used in the current model. This value will be similar to the final estimate of the parameter in the prior problem. (A difference from this prior estimate may result, reflecting the fact that if there is a sequence of models, each a submodel of the one succeeding it, then the data used to fit the last of the models may bear on a parameter peculiar to the first model, if only slightly.) If indeed the parameter’s value were not fixed, but the parameter’s estimate obtained along with the estimates of the non-prior-specific parameters, then the latter estimates and the minimum value of the objective function would be unchanged, but the search for the parameter’s estimate would require additional computer time. With the value MODE=1, in fact, the parameter’s estimate is obtainable posthoc, after the search, and this additional computer time is saved. The output from the Covariance Step includes all the usual type of information about the estimators of the prior-specific parameters along with that about the estimators of the non-prior-specific parameters. However, note that an output model specification file will not contain any information about the estimates of prior-specific parameters.

2: All parameters with a nonfixed value in the prior problem but a fixed value in the current problem are identified as being shared parameters (i.e. the frequency prior is adjusted). The final estimates of these parameters are their fixed values.

ITYP: Relevant only if TNPRI is called during a Simulation Step. Values are:

0: The value of the UCP (unconstrained parameter) vector is obtained from simple random sampling.

1: Within the given problem, TNPRI is to be called a specified number (NSAM) of times to obtain this number of different values of the UCP vector. These values are obtained by generating a Latin sample of size NSAM from equiprobable partitions of an ellipsoid in UCP space (hyper-ellipsoidal sampling), followed by sampling a point “uniformly” from each partition. This scheme may be used, for example, when the problem has NSAM subproblems, in which case, TNPRI would be called NSAM times, once each time during the problem when ICALL=4, and at each of these calls, a different random value of the UCP vector will be produced.

2: Just as with value 1, but the NSAM values are obtained by generating a Latin sample of size NSAM from equiprobable partitions of an ellipsoid in UCP space (hyper-ellipsoidal sampling), followed by taking the “center point” of each partition.

In all three cases, a value of the UCP vector is transformed into THETA-OMEGA-SIGMA space.

After each call to NWPRI, the simulated values for THETA, OMEGA and SIGMA may be found in global variables and thus they are communicated directly to NONMEM. (See Section 53.87).

PLEV: When TNPRI is being used at ICALL=0 or 1, but TNPRI will not be used at ICALL=4 (i.e. during the Simulation Step), PLEV should be set to 0. When it is being used at ICALL=2, PLEV should also be set to 0. When TNPRI is being used at ICALL=0 or 1 and will also be used at ICALL=4, or when it is being used at ICALL=4, then PLEV must be set to a fraction strictly less than 1, e.g. 0.999. PLEV is double precision with NONMEM 7, and is single precision with NONMEM VI.

A UCP value will actually be obtained using a truncated multivariate normal distribution, i.e. from an ellipsoidal region $R_1$ over which only a fraction of mass of the normal occurs. This fraction is given by PLEV. Simple random sampling occurs in $R_1$. Latin sampled partitions are partitions of $R_1$.

NSAM: Relevant only if TNPRI is called during a Simulation Step. Consider two cases.

Latin hyper-ellipsoid sampling is used with ITYP=1, or b) simple random sampling along with the adjustment for small sample correlation effect is used (See next input argument). In case a) NSAM should equal the exact total number of different values of the parameter vector that must eventually be produced over the entire NONMEM problem. In case b) NSAM should be no less than this number.

ISS: Relevant only if TNPRI is called during a Simulation Step. A value of the UCP vector is obtained by first sampling a value from the standard multivariate normal distribution - called here “the standard value” and then transforming this value to one from the appropriate multivariate normal. The correlation matrix of the standard normal is the identity matrix. When NSAM is small, the estimated correlation matrix from the sampled standard values might not be quite close to the identity matrix - this is here called “the small sample correlation effect”.

1: An adjustment is made for the small sample correlation effect, by first transforming the NSAM standard values altogether into new values which are very nearly standard multivariate normal values and such that the sample correlation matrix of these new values is exactly the identity matrix.

0: No adjustment is made for the small sample correlation effect.

IVAR: Relevant only if TNPRI is called during a Simulation Step. When TNPRI simulates a parameter value, this value is that of an estimate of the parameter from possible data under the prior model, and this simulation is based on asymptotic statistical theory. The partial derivatives comprising the R and S matrices referred to below are taken with respect to the UCP (See Section 50.8). The possible values of IVAR are:

0: The variance-covariance matrix of the multivariate normal on the UCP is taken to be two times the inverse R matrix obtained from the prior problem. This could be appropriate if the estimated asymptotic variance-covariance matrix from the prior problem is also based only on the R matrix.

1: The variance-covariance matrix of the multivariate normal on the UCP is taken to be Rinv*S*Rinv matrix from the prior problem.

2: The variance-covariance matrix of the multivariate normal on the UCP is taken to be S matrix from the prior problem. (NM75)

Briefly:

IVAR=0: Uses Rinv of a former problem.
IVAR=1: Uses Rinv*S*Rinv of a former problem.
IVAR=2: Uses S of a former problem.

Output argument:

CNT: Relevant only if TNPRI is called at ICALL=2. CNT is the penalty.

REFERENCES: Guide II, section D.2.5