84  $COV Record

1 $COV: Unconditional Evaluation

The $COV step can be performed unconditionally even when an estimation terminates abnormally, by specifying:

$COV UNCONDITIONAL

2 $COV: Additional Options and Behavior

Example syntax:

$COV UNCONDITIONAL TOL=10 SIGL=10 SIGLO=11 NOFCOV ATOL=6 RESUME

If $COV is specified, then for IMP, IMPMAP, and ITS methods, standard error information will be supplied for every $EST statement.

Standard error information for the classical methods (METHOD=0, METHOD=1) will be given only if they are the last estimation method, and only if NOFCOV is not specified.

If UNCONDITIONAL is specified, then for the IMP and IMPMAP EM methods, if the R information matrix is not positive definite, the program will modify the matrix to be positive definite, will report that it has done so, and provide the standard errors. The user should use the standard error results with caution should a non-positive definite flag occur.

The ITS and SAEM methods can only evaluate the S matrix, and will do so even if MATRIX=R is requested. The banner information will show what type of variance was evaluated.

The BAYES method always supplies standard errors, correlation matrix, and covariance matrix, even when $COV step is not requested, as these results are a direct result of summarizing the accumulated NITER samples. Furthermore, the matrices are always positive definite, and therefore always successful.

To obtain the eigenvalues to the correlation matrix, even for the BAYES method, a $COV step must be issued with the PRINT=E feature.

2.1 TOL, SIGL, SIGLO (NM72)

The TOL (used by PREDPP when differential equations are integrated) and SIGL and SIGLO may be set specifically for the $COV step, distinct from those used during $EST. This special option for $COV is not so important for the new EM or BAYES methods, which are able to obtain suitable standard errors using SIGL, SIGLO, and TOL that are also used for estimation, but classical NONMEM methods in particular can require a different significant digits level of evaluation (usually more stringent) during the $COV step than during $EST. Keep in mind that when evaluating the R matrix, SIGL and TOL should be at least 4 times that of what one would normally set NSIG. If evaluating only the S matrix, then SIGL, SIGLO, TOL should be at least 3 times that of what one normally sets NSIG. For example, during $EST, NSIG=2, SIGL=6, TOL=6 may be sufficient, but during $COV, you may need SIGL=12 TOL=12 to avoid positive definiteness issues. The MATRIX, TOL, and SIGL have no relevance to the variance results for a BAYES method, which are derived from samples generated during the estimation step. If TOL is set in the $COV record, but SIGL and/or SIGLO are not, then the TOL is not changed. Also, if TOL is set for the $COV record, then this TOL is used for all compartments.

2.2 FAST (NM74)

This is equivalent to FAST for the $EST record (see I.50 The FAST Option for use with FOCE/ITS and Differential Equation ($DES) Models (NM74)), but for the $COV record. If $EST FAST is set, then $COV will be set to FAST, unless you specify SLOW or NOSLOW at the $COV record.

2.3 ATOL (NM72)

The absolute tolerance option pertains to using ADVAN13, and as of NM73, to ADVAN9 as well, where ATOL is the accuracy for derivatives evaluated near zero. The same ATOL value is set for all compartments. The ATOL by default is 12. Usually the problem runs quickly when using ADVAN13 with this setting. On occasion, however, you may want to reduce ATOL (usually equal to that of TOL), and improve speeds of up to 3 to 4 fold. ATOL may be set at the $EST or $COV command. Keep in mind that ATOL is changed for the $COV step only if SIGL and/or SIGLO are also specified at the $COV record.

2.4 KNUTHSUMOFF (NM74)

In NONMEM 7.4, the Knuth summing method is used to allow the most accurate summation of individual objective function values, even with large variations in values of the individual objective function. To turn this off, and allow a standard summation (not recommended except for comparison purposes from earlier versions), set KNUTHSUMOFF=1. If KNUTHSUMOFF was set in the $EST step, but not in the $COV step, the KNUTHSUMOFF value of the last $EST record will be used.

2.5 POSDEF=0 (default) (NM75)

If POSDEF=1, then if the R matrix is not positive definite, then during its inversion, it will be forced positive definite, by method of shifting the eigenvalues by an amount slightly greater than the most negative eigenvalue. If POSDEF=2, then R matrix is made positive definite by making the negative valued eigenvalues approximately 1/100 of the least positive eigenvalue (courtesy of Larry Schaeffer, U. Guelph, Ontario, Ca, see ..\guides\PDforce.pdf). If POSDEF=3, then R matrix is made positive definite by using the absolute values of the eigenvalues. POSDEF is 0 (no correction) by default for classical covariance step, and 3 for EM methods. The preconditioning method, described below, is a more sophisticated, but time-consuming method of making the R matrix positive definite. The NOHABORT or NOABORT option during the $EST step must be set in order for POSDEF to function during an FO/FOCE/Laplace assessment of the $COV step.

2.6 CHOLROFF=0 (default) (NM74)

If CHOLROFF is set to 1, then one part of the R matrix evaluation (before inversion) will be evaluated in the manner of earlier versions of NONMEM. This is strictly for comparison with earlier versions for diagnostic purposes. If CHOLROFF=0, then a part of the R matrix evaluation will undergo cholesky decomposition, and if not positive definite, only then will it undergo evaluation according to earlier versions of NONMEM. If CHOLROFF=2, then if cholesky decomposition fails, it will be slightly adjusted to be positive definite. The CHOLROFF=0 is the safest option, in that the Cholesky decomposition (if positive definite) provides more precision in evaluating the R matrix. The CHOLROFF=2 setting is useful if you wish to increase the rate of success in obtaining an R matrix (followed by its inverse) that could be suitable for SIR sampling. If you use CHOLROFF=2, then positive definiteness has been corrected before the inverse occurs, so POSDEF will have no additional effect. If you use POSDEF, then CHOLROFF should be 0 or 1.

2.7 NOFCOV (NM72)

No $COV step for any classical estimation steps. This would be useful if you wanted EM estimation analyses with variance-covariance assessment performed, and a final FOCE analysis performed, but did not want the program to spend time on standard error assessments for FOCE, which can take a long time relative to the other methods.

2.8 RESUME (NM73)

If an MSFO=msffile specification was made in the $EST step, and analysis was interrupted during the $COV step for the FO/FOCE/Laplace method, then the $COV step may be resumed where it was interrupted by executing another control stream file that uses the $MSFI record specifying the MSFO file of the interrupted analysis, and the RESUME option is entered at the $COV record:

…
$MSFI=msffile
…
$COV RESUME

In addition to the main msf file and _ETAS file, files _RMAT, and _SMAT files will generated by the MSFO option of the previous control stream or problem, which are needed by the RESUME option of the new control stream or problem.

2.9 PARAFPRINT=1 (default, NM74)

The print iteration intervals to the parallelization log file can be controlled by this option during parallelization of the $COV step.

2.10 PARAFILE=OFF (NM75)

Specifically turn off parallelization during the $COV step.

2.11 THBND=1 (default) (NM74)

For covariance assessment for FOCE/FO/Laplace, if THBND=1, any theta boundaries specified in the $THETA record causes NONMEM to impose a non-linear transformation of the theta parameters so that the transformed parameters may vary from –infinity to infinity. It does this with logistic transformations. This is suitable during the estimation step, but may be desirable to have this off (THBND=0), and assess partial derivatives of the objective function with respect to the THETAs themselves, or some linear transformation of these THETAs. By default THBND=1, in keeping with the behavior of earlier NONMEM versions, which effectively has THBND=1. Usually boundaries that are fairly wide will not impact how the variance-covariance is assessed, such as when a lower bound of 0 is given, but if you have very narrow boundaries set, then this can impact the assessment of the variance-covariance of the estimates considerably, and you may wish to set THBND=0. If no lower or upper bounds are given to THETAs in $THETA record, this option has no impact.

The degree to which boundary settings can impact the variance of a parameter is linearly approximated by the relation:

\[I_c = I_u + \left[\frac{(H-L)^2}{(H-x)(x-L)}\right]\]

Where I refers to information (inverse variance), c=boundary constrained assessment (THBND=1), u=non-boundary constrained assessment (THBND=0), x is estimate of theta, H is high (upper) boundary, L is lower boundary. For the standard error to be decreased by no more than 10%, the second term should not be larger than 20% (as this is a 1/variance term), so no larger than 0.2. Assuming boundaries equidistant to the estimate, this comes to sqrt(4/0.2)=4.47, that is, the distance between estimate and nearest boundary should be no smaller than 4.47 times the unconstrained standard error.

IMP/SAEM/ITS do not evaluate the variance/covariance taking into account boundary settings, and act always effectively as THBND=0.

2.12 Importance Sampling of the Variance-Covariance of the Parameter Estimates (NM74)

Based on the Perl-Speaks-Nonmem Sampling-Importance-Resampling algorithm (SIR), after the variance-covariance matrix is assessed by the $COV step by the usual deterministic finite-difference evaluation of the Fisher information matrix for FO/FOCE/Laplace, it can be used as the variance-covariance to a proposal density to obtain Monte Carlo importance samples of population parameters, in the neighborhood of the minimum of the FOCE objective function. The samples are generated with normal distribution sampling in the unconstrained parameter domain (see Appendix K of NONMEM7_technical_guide.pdf describing unconstrained parameter domain), and then transformed back to the original domain. The Monte Carlo samples will be listed in the raw output file (.ext), along with their importance sampling weights. The weighted average parameter values are listed on the line labeled iteration -1000000000, and the weighted sample variance-covariances are listed on the line labeled -1000000001, in accordance with raw output file format. The other statistics labeled -100000000x, are also reported, and should be similarly interpreted as if they were results reported for a BAYES analysis. Iteration 0 contains the objective function minimum position, along with its importance sampling weight, which is 1. Re-sampling is not performed, however this can be done by a post-processing utility software (such as R), keeping in mind that each sample should be weighted according to its relative WEIGHT listing. For example, creating a uniform random variable r, when the sum of the normalized weights of the jth sample in the list first is equal to or exceeds r, that is, when:

\[\frac{\sum_{k=1}^{j} W_k}{(\text{sirsample}+1) W_{-1000000000}} \geq r\]

select that sample j (sirsample+1 is the number of total samples, including the 0th iteration one).

In this way, samples will be selected in proportion to their weight listings. Similarly, to determine a quantile position q for a particular parameter, sort the samples for that parameter, and adding their weights in order of the sorted samples, when the normalized summed weight first exceeds q:

\[\frac{\sum_{k=1}^{j} W_k}{(\text{sirsample}+1) W_{-1000000000}} \geq q\]

select that parameter as the qth quantile. Read about the utilities that perform frequency and quantile sorting (table_quant), or resampling (table_resample) in I.81 table_quant, and table_resample Utility Programs for Analyzing $COV Sampling-Importance-Resampling (SIR) data (NM74).

A simple example of the use if SIRSAMPLE is given in examples\sirsample.ctl.

The following options may be set for the $COV record to control the SIR sampling process. Note: the RESUME option only resumes an interrupted generation of finite difference of the covariance estimation. If your analysis is interrupted during a SIR sampling process, RESUME will start at the beginning of a new set of sample generation.

2.13 SIRSAMPLE=0 (default) (NM74)

By default SIRSAMPLE=0, so SIR process does not occur. Setting SIRSAMPLE to a value greater than 0 will produce SIRSAMPLE importance samples. These will be placed in the raw output file as a table with heading “Importance Sampling of Variance-Covariance (SIR): Goal Function=AVERAGE FITNESS OF PROPOSAL DENSITY”.

As of nm75, SIRSAMPLE may contain a list of integers, one for each of the SRNITER iterations (use double quotes around the list):

SIRSAMPLE="1000,200,500"

If there are fewer than SIRNITER values, the last value in the list will be used as the SIRSAMPLE value for the remaining iterations.

2.14 SIRNITER=1 (default) (NM74)

The number of times you want SIR sampling to be performed. While the proposal density for the first iteration comes from the main $COV step, subsequent iterations borrow the proposal density variance from the SIR sample variance of the previous iteration, centered around the SIR average of the parameters.

2.15 SIRCENTER=0 (default) (NM74)

Where the sampling (proposal) density is to be centered. On the first iteration, the mean of the sampling density is at the estimate. On subsequent iterations, the mean of the sampling density is at the estimate (SIRCENTER=0) or at the mean of the (transformed) samples of the previous iteration (SIRCENTER=1).

2.16 SIR_CAPCORR=1.0 (default) (NM75)

Limit the correlation of OMEGAs and SIGMAs generated from the proposal density to have |correlation| of SIR_CAPCORR or less. This is similar to PSN’s cap_correlation option.

2.17 SIRMINWT=0.001, SIRMAXWT=1000.0 (default) (NM75)

Limit the weight range that a sample may have relative to the proposal density, to avoid extreme values.

2.18 IACCEPT=1 (default) (NM74)

The acceptance ratio acts similarly to importance sampling in EM analysis. For objective function profiles that are particularly heavy tailed (this could happen when there are relatively few subjects, and or data), then IACCEPT=0.4 may be more suitable, so the proposal density variance is made to be greater than the original variance-covariance, and more sampling is done farther away from the center (objective function minimum).

2.19 IACCEPTL=0 (default) (NM74)

The IACCEPTL option performs the same as listed for IACCEPTL=0 (default) (NM74) in section I.38 Monte Carlo Importance Sampling EM. It adds a normal density sampling scaled according to the IACCEPTL, alternating with normal density sampling scaled according to IACCEPT.

2.20 SIRDF=n (NM74)

The proposal density is to be a t distribution with n degrees of freedom. Default is 0, a normal density. You may wish to utilize a t-distribution with SIRDF degrees of freedom, to provide heavy tail sampling. This has a similar purpose as setting IACCEPT less than 1.

2.21 SIRSEED=11456 (default) (NM75)

Starting seed for SIR analysis.

2.22 SIRCLOCKSEED=0 (default) (nm75)

As of nm75, the actual starting seed will be 10000*(seconds after midnight)+SIRSEED (SIRSEED may be set to 0 for this purpose), if SIRCLOCKSEED=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.

2.23 RANMETHOD=blank (default) (NM74)

The RANMETHOD settings for the SIR analysis perform the same actions as described for RANMETHOD in the Importance Sampling EM section. See RANMETHOD=[n|M|S|m|P] (default n=3) in I.38 Monte Carlo Importance Sampling EM. If RANMETHOD with a ‘P’ attribute is set, then parallelization occurs by the usual parsing of subjects among the nodes, and retaining the seed pattern for each subject, so that the sample pattern is re-produced regardless of how the job is split up. As of NM75, if RANMETHOD does not have the ‘P’ attribute, then parallelization occurs by parsing the SIRSAMPLE objective function evaluations evenly among the nodes. In such a case, the sampling pattern cannot be made to be reproducible. There are efficiency advantages to splitting the job up among the samples rather than among subjects within each objective function (OFV) evaluation, if the OFV is rapidly calculated, or there are fewer subjects for parsing than there are nodes available.

2.24 SIRPARAFPRINT=1 (default, NM75)

The print iteration intervals to the parallelization log file can be controlled by this option during parallelization of the SIR analysis.

2.25 SIRPARAFILE=OFF (NM75)

Specifically turn off parallelization during the SIR analysis.

2.26 SIRPRINT=0 (default) (NM74)

Set the console print iterations interval. This does not impact the iterations listed in the raw output file.

2.27 FILE=blank (default) (NM74)

By default, the raw output file is whatever was listed in the $EST step, or root.ext, where root is the root name of the control stream file. You can re-direct SIR sample listings to an alternative file with this option.

2.28 FORMAT=s1PE12.5 (default) (NM75)

By default, the raw output file format is whatever was listed in the $EST step, or s1PE12.5. You can change its format with the above option. 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 (such as "OMEGA(2,1)").

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

2.29 SIRTHBND=THBND (default) (NM74)

As with the deterministic covariance assessment step, when SIRTHBND=1, the transformed parameters are sampled, so that no sample is below the $THETA lower bound specification, and no higher than the $THETA upper bound specification. To allow a boundariless search in the original theta domain, SIRTHBND should be 0. By default, SIRTHBND is the value of THBND, which in turn is 1 by default.

2.30 Using $RCOV with SIR sampling (NM75)

You can use $RCOV (or $RCOVI) (see [I.69 $RCOV and \(RCOVI Record For Inputting Variance-Covariance information from another problem (NM75)](init-records.qmd)) to introduce a covariance of estimates matrix from a `*.cov` file of any previous analysis, or generated by any other means. The following shows how one could generate the covariance matrix by the usual ``\)COV`step from one control stream problem, sayrcov.ctl`:

…
$EST METHOD=1 INTERACTION NSIG=3 PRINT=10 NOABORT MAXEVAL=9999
         FORMAT=,1PE23.15
$COV MATRIX=R PRINT=E UNCONDITIONAL

Which generates rcov.cov, rcov.ext, and rcov.phi. These files in turn can be used in a subsequent analysis (Suppressing estimation with FNLETA=2):

…
$CHAIN FILE=rcov.ext NSAMPLE=0 ISAMPLE=-1000000000 FORMAT=,1PE23.15
$RCOV FILE=rcov.cov TBLN=1 DELIM=,
$ETAS FILE=rcov.phi TBLN=1 DELIM=,
$EST METHOD=1 INTERACTION NSIG=3 PRINT=10 NOABORT MAXEVAL=9999 FNLETA=2
         FORMAT=,1PE23.15
$COV MATRIX=R UNCONDITIONAL SIRSAMPLE=1000 NITER=3

Thus the $RCOV record serves the same purpose as the option -covmap_input in Perl-Speaks-NONMEM.

2.31 Preconditioning the R Matrix to Improve Precision and Success Rate of $COV Step (NM74)

Preconditioning of the R matrix, as described by Aoki et al (Aoki et al. 2016), and implemented in Perl Speaks NONMEM, can improve the precision, as well as increase the chance of success of the $COV step.

2.31.1 PRECOND=0 (default) (NM74)

By default, PRECOND is 0, and no preconditioning of the R matrix is performed. When PRECOND=n, then up to n preconditioning cycles are performed. This is used in combination with the PFCOND setting.

2.31.2 PRECONDS=TOS (default) (NM74)

By default, if preconditioning is performed, it is done on THETAs (T), OMEGAs (O), and SIGMAs (S). Specify PRECONDS=T to do only THETAs, PRECONDS=TO to do only THETAs and OMEGAs, etc.

2.31.3 PFCOND=0 (default) (NM74)

PFCOND means “forced” preconditioning. Preconditioning occurs exactly PFCOND times, without testing if the R matrix is positive definite or not on each preconditioning cycle. On the remaining PRECOND-PFCOND cycles, the R matrix is tested for positive definiteness, and upon success, will terminate the preconditioning cycles.

2.31.4 PRETYPE=0 (default) (NM74)

By default (PRETYPE=0), the R matrix corrector is V*square_root(eigenvalue), as described in Aoki. If you set PRETYPE=1, then corrector is V*square_root(eigenvalue)*Vtranspose. If you set PRETYPE=2, then corrector is the correlation version of PRETYPE=1.

2.31.5 FPOSDEF=0 (default) (NM74)

If FPOSDEF>0, then if the R matrix is not positive definite, it will be forced positive definite. If PRECOND>0, this will occur after the PRECONDth try. The FPOSDEF operates on preconditioned R matrix, while POSDEF described earlier operates on any non-preconditioned R matrix, that is, when MATRIX when PRECOND=0.

The positive definite algorithm used is as follows. If FPOSDEF=1, then if the R matrix is not positive definite, it will be forced positive definite, by method of shifting the eigenvalues by an amount slightly greater than the most negative eigenvalue. If FPOSDEF=2, then R matrix is made positive definite by making the negative valued eigenvalues approximately 1/100 of the least positive eigenvalue. If FPOSDEF=3, then R matrix is made positive definite by using the absolute values of the eigenvalues. Setting FPOSDEF=3 appears to be the best method here.

3 A Note on Covariance Diagnostics

There are several conditions that can occur in assessing the variance-covariance matrix of the estimates, which are best defined according to eigenvalues that it detects in them.

  1. Positive definite means there are only positive eigenvalues. NONMEM outputs proper variance-variance matrices.
  2. Non-positive definite means there is at least one eigenvalue that is less than or equal to zero.
  3. Positive-semidefinite means there are no negative eigenvalues, but at least one zero valued eigenvalue (singular).
  4. Non-positive-semidefinite means there is at least one negative eigenvalue.
  5. Non-positive-semidefinite and singular means there is at least one negative eigenvalue, and at least one zero valued eigenvalue. Non-inverted matrices may be outputted by NONMEM.
  6. Non-positive-semidefinite and non-singular means there is at least one negative eigenvalue, and no zero valued eigenvalue. Alternative diagnostic matrices may be outputted by NONMEM.
  7. Negative-definite means there are only negative eigenvalues.
  8. Non-negative-definite means there is at least one eigenvalue that is greater than or equal to zero.

NONMEM tests for conditions 1), 5), and 6), and outputs appropriate result matrices, or diagnostic matrices, as it is able.

Alternative expressions would be unsuitable to describe the condition of the matrices. For example, non-positive-definite (2) does not mean the same as positive-semi-definite (3). Similarly, non-positive-definite (2) is not exactly the same as non-positive-semidefinite (4). The set of non-negative-definite matrices (8) includes matrices that are positive-definite (1), positive-semi-definite (3), and a subset of non-positive-semidefinite (4) not including those with all negative eigenvalues.