72  EM and Monte Carlo Methods

1 Introduction to EM and Monte Carlo Methods

Expectation-maximization methods use a two step process to obtain parameters at the maximum of the likelihood. In the expectation step, the THETAs, OMEGAs, and SIGMAs are fixed, while for each individual, expected values (conditional means) of the eta’s and their variances are evaluated. If necessary, expected values of gradients of the likelihood with respect to the THETAs and SIGMAs are also evaluated, integrated over all possible values of the etas. From these constructs, the THETAs and SIGMAs are updated during the maximization step using these conditional means of the etas and/or the gradients. The OMEGAs are updated as the sample variance of the individual conditional means of the etas, plus the average conditional variances of the etas. The maximization step is therefore typically a single iteration process, requiring very little computation time. The more accurately these constructs are evaluated during the expectation step, the more accurately the total likelihood will be maximized.

2 Iterative Two Stage (ITS) Method

Iterative two-stage evaluates the conditional mode (not the mean) and first order (expected) or second order (Laplace) approximation of the conditional variance of parameters of individuals by maximizing the posterior density. This integration step is the same as is used in FOCE or Laplace. Population parameters are updated from subjects’ conditional mode parameters and their approximate variances by single iteration maximization steps that are very stable (usually converging in 50-100 iterations). Because of approximations used, population parameters almost, but not quite, converge towards the linearized objective function of FOCE. Iterative two stage method is about as fast as FOCE with simple one or two compartment models, and when set up with MU referencing (described below) can be several fold faster than FOCE with more complex problems, such as 3 compartment models, and differential equation problems.

The iterative two stage method is specified by

$EST METHOD=ITS INTERACTION NITER=50

where NITER (default 50) sets maximum number of iterations. For all new methods, it is essential to set INTERACTION if the residual error is heteroscedastic.

3 Monte Carlo Importance Sampling EM

Importance sampling evaluates the conditional (posterior) mean and variance of parameters of individuals (etas) by Monte Carlo sampling (integration, expectation step). It uses the posterior density which incorporates the likelihood of parameters relative to population means (THETAs) and variances (etas) with the individual’s observed data. By default, for the first iteration, the mode and first order approximation of the variance are estimated (called mode a posteriori, or MAP estimation) as is done in ITS or FOCE, and are used as the parameters to a normal distribution proposal (sampling) density. From this proposal density Monte Carlo samples are generated, then weighted according to the posterior density as a correction, since the posterior density itself is generally not truly normally distributed, and conditional means and their conditional variances are evaluated. For subsequent iterations, the normal density near the mean of the posterior (obtained from the previous iteration) is used as a proposal density. Population parameters (THETAs, SIGMAs, and OMEGAs) are then updated from subjects’ conditional mean parameters, gradients, and their variances by single iteration maximization steps that are very stable, and improve the objective function. The population parameters converge towards the minimum of the objective function, which is an accurate marginal density based likelihood (exact likelihood). A series of options defined at the $EST command are available to the user to control the performance of the importance sampling, such as the number of Monte Carlo samples per individual (ISAMPLE), and scaling of the proposal density relative to the posterior density (IACCEPT). Termination criteria (CITER, CALPHA, CTYPE, and CINTERVAL) may also be set, which are explained in detail in a later section. Typically, 300 Monte Carlo samples are needed, and 50-200 iterations are required for a randomly stationary objective function, that is, when the objective function does not vary in a directional manner beyond the Monte Carlo fluctuations.

The Importance sampling method is specified by

$EST METHOD=IMP INTERACTION

Followed by one or more of the following options:

3.1 NITER/NSAMPLE=50

Sets maximum number of iterations (default 50). Typically, 50-100 iterations are need to for a problem to have a randomly stationary objective function.

3.2 ISAMPLE=300

Sets number of random samples per subject used for expectation step (default 300). Usually 300 is sufficient, but may require 1000-3000 for very sparse data, and when desiring objective function evaluation with low Monte Carlo noise.

3.3 ISAMPEND=n, STDOBJ=d (NM73)

For importance sampling and direct sampling only, if ISAMPEND is specified as an integer value greater than ISAMPLE, and STDOBJ is set to a real value greater than 0, then NONMEM will vary the number of Monte Carlo samples under each subject between ISAMPLE and ISAMPEND, until the stochastic standard deviation of the objective function falls below STDOBJ.

3.4 IACCEPT=0.4

Expand proposal (sampling) density variance relative to conditional density so that on average conditional density/proposal density=IACCEPT (default 0.4). For very sparse data or highly non-linear posterior densities (such as with categorical data), you may want to decrease to 0.1 to 0.3.

3.5 IACCEPT=0.0 (NM7.3)

For importance sampling only, you may set IACCEPT=0.0, and NONMEM will determine the most appropriate IACCEPT level for each subject, and if necessary, will use a t-distribution (by altering the DF for each subject) as well. If IACCEPT=0, the individual IACCEPT values and DF values will be listed in root.imp, where root is the name of the control stream file.

3.6 ISCALE_MIN=0.1 (defaults for IMP, NM72), ISCALE_MAX=10.0 (NM72)

In importance sampling, the scale factor used to vary the size of the variance of the proposal density in order to meet the IACCEPT condition, is in NM72 by default bounded by ISCALE_MIN of 0.1, and ISCALE_MAX=10.0. On very rare occasions, the importance sampling objective function varies widely, and the scale factor boundary may need to be reduced (perhaps ISCALE_MIN=0.3, ISAMPLE_MAX=3). After the importance sampling estimation, remember to revert these parameters to default operation on the next $EST step:

ISCALE_MIN=-100 ISCALE_MAX=-100.

Note: the values to ISCALE_MIN and ISCALE_MAX for the IMP method in NONMEM 7.1 and earlier were 0.01, 100, respectively, and were not changeable by the user.

3.7 SEED=11456 (default)

The seed for the random number source used in Monte Carlo integration is initialized (default seed is 11456).

3.8 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.

3.9 MAPITER=1 (default) (NM72)

By default, MAP estimation is performed only on the first iteration, to obtain initial conditional values (modes and approximate variances) to be used for the sampling density. Subsequently, the Monte Carlo assessed conditional means and variances from the previous iteration are used as parameters to the sampling density. However, the user can select the pattern by which MAP estimations are intermittently done, and their conditional statistics used for the sampling density. MAPITER=n means the first n iterations are to use MAP estimation to assess parameters for the sampling density. After these n iterations, the Monte Carlo conditional means and variances of the pervious iteration are used for the sampling density parameters of the present iteration. If MAPITER=0, then the first iteration will rely on conditional means and variances that are in memory. These may have come from an MSF file, or from a previous estimation step.

3.10 MAPINTER=0 (default) (NM72)

Every nth iteration, the MAP estimation should be used to provide parameters to the sampling density. Thus, if MAPITER=20 and MAPINTER=5, then for the first 20 iterations, MAP estimation is used, and thereafter, every 5th iteration the MAP estimation is used. If MAPINTER=-1 (NM73), then mapinter will be turned on only if the objective function increases consistently over several iterations.

Setting an option to -100 will force NONMEM to select the default value for that parameter.

3.11 MAPCOV=1 (default) (NM74)

For iterations for which the MAP estimation is performed, by default (MAPCOV=1), the MAP estimated mode is used as the center (mean) for the sampling density, and the first order (or second order if Laplace option is used) approximate conditional variance is used as the variance of the sampling density. If MAPCOV=0, then only the mode is used for the sampling density’s center, and the Monte Carlo assessed variance of the previous iteration is used as the sampling density’s variance. If MAPCOV=2, then the Monte Carlo assessed conditional mean of the previous iteration is used, but the MAP first order (or second order if Laplace option) assessed variance is used for the sampling density. This option has been added for experimental purposes, and has no value for the user. It should be left at its default value of 1.

3.12 EONLY=1

Evaluate the objective function by performing only the expectation step, without advancing the population parameters (default is 0, population parameters are updated). When this method is used, NITER should equal 5 to 10, to allow proposal density to improve with each iteration, since mean and variance of parameters of normal or t distribution proposal density are obtained from the previous iteration. Also it is good to get several objective function values to assess the Monte Carlo noise in it.

As of NM74, if EONLY=2, then not only are the population parameters not updated with each iteration, neither are the individual conditional means/modes and conditional variances until the last iteration, if the variances of the population parameter estimates are estimated ($COV is requested). If EONLY=3, then the conditional modes and approximate variances from MAP estimation will be saved on the first iteration, and used for the sampling density of all subsequent iterations. This improves efficiency when selecting MAPITER=1, MAPINTER=0, so that the MAP estimation did not need to be repeatedly performed.

3.13 DF=4

The proposal density is to be t distribution with 4 degrees of freedom. Default DF=0 is normal density. The t distribution has larger tails, and is useful for situations where the posterior density has a highly non-normal distribution. For very sparse data or highly non-linear posterior densities (such as with categorical data), you may want to set DF to somewhere between 2 and 10.

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

Where

n=0-4
m=0-3

By default, the random number generator used for all Monte Carlo EM and Bayesian methods use the Knuth method, 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 (n) as follows (n=0-4):

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 in $SIMULATION

For special purposes, a Sobol (Press et al. 1992) sequence method with or without scrambling (Hong and Hickernell 2003) may be called upon, and only for the purpose of creating quasi-random samples of eta vectors. To select the Sobol method without scrambling, add an S to RANMETHOD. For example,

RANMETHOD=2S

Selects random number generator ran2 for general purposes, and Sobol sequence for the eta vector generation. The number m is reserved for the type of scrambling desired (m=0-3):

0: no scrambling (so S0 is the same as S)

1: Owen type scrambling

2: Faure-Tezuka type scrambling

3: Owen plus Faure-Tezuka type scrambling.

Other examples:

RANMETHOD=S1

Indicates Sobol sequence with Owen scrambling for eta vector generation. Since there is no integer in the first position of RANMETHOD indicated, the general random number generator remains unchanged from the RANMETHOD specification previously specified, or ran method 3, if none was specified earlier.

RANMETHOD=1S2

Indicates ran1 type random number generator for general purposes, Sobol sequence with Faure-Tezuka scrambling for eta vector generation.

The Sobol sequence method of quasi-random number generation can reduce the Monte Carlo noise in the objective function evaluation during importance sampling under some circumstances. When the sampling density fits the posterior density well, such as with rich, continuous data, the Sobol sequence method does not reduce the Monte Carlo noise by much. If you are fitting categorical data, or sparse data, and perhaps you are using the t distribution (DF>0) for the importance sampling density, then Sobol sequence generation may be helpful in reducing Monte Carlo noise. The RANMETHOD specification propagates to subsequent $EST records in a given problem, but does not propagate to $CHAIN or $TABLE records.

In NM72, only DIRECT and IMP/IMPMAP methods could utilize the Sobol quasi-random method. As of NM73, Sobol may be used for BAYES and SAEM methods as well. From experience, The S0 and S1 methods produce considerable bias for SAEM and BAYES, whereas S2 and S3 perform better.

As of NM73, if you add a P descriptor to RANMETHOD, such as

RANMETHOD=P
RANMETHOD=3P
RANMETHOD=3S2P

then each subject will receive its own random number sequence, that will stay with that subject regardless of whether the job is run as a single process or parallel process. This assures that stochastically similar answers will be obtained for Monte Carlo estimation methods, regardless of the number of processes or different kinds of parallelization setups used to solve the problem. There is additional memory cost in using this option because the seed and seed status (additional internal variables of the random number generator that establish the random sequence) must be stored for each subject, and for SOBOL/QR sampling there may even be a reduction in speed because the random sampling algorithm has to be re-set for each subject. To reiterate, a single job run without the P descriptor will not be stochastically similar to a single job run with the P descriptor (although they will be statistically similar), or to any parallel job run. But, a single job run using the P descriptor will be stochastically similar to any parallel job run also using the P descriptor. If maintaining stochastic similarity regardless of how the job is run (single or any parallel profile) is important to you, then always set the P descriptor (so, RANMETHOD=P, at least).

3.14.1 Note on the t-Distribution Sampling Density (DF>0), and its Use With Sobol Method (RANMETHOD=S)

When using the t-distribution sampling density (DF>0), by default the algorithm creates a random vector from n independent univariate t-distributed samples. This is called the U algorithm, and the most efficient use of the U type t-distribution is when DF=1,2,4,5,8, or 10. These algorithms were designed to work well with the Sobol method’s ability to reduce Monte Carlo noise.

As of NM74, another way of producing the vector is from a multi-variate t-distribution algorithm (suggested by Robert Leary), which can be selected by placing an M in the RANMETHOD descriptor, placed after the random number method number, for example:

RANMETHOD=3MS2

(the default setting is U for composite univariate). The multivariate t-distribution algorithm (M) produces samples that have radially symmetric densities (that is, the density is a function of the sum of squares of vector elements) may provide a better fitting sampling density for some kinds of models, and hence, more efficient sampling. However, the individual random elements in the vector are not statistically independent, and when used with the Sobol method, does not result in as much reduced Monte Carlo noise as when the composite univariate t-distribution vector (U) is used. An alternative method is by a mixture of two normal densities, using IACCEPTL.

3.15 IACCEPTL=0 (default) (NM74)

If IACCEPTL is set to greater than 0 then NONMEM uses this value as a scale to a second multi-variate normal density, to cover long tails in the posterior density (hence L for long tails), in combination with the ACCEPT value to cover the posterior density near the mode. For one half of ISAMPLE samples, IACCEPT is used to scale a multivariate-normal proposal density, and for the other half of ISAMPLE samples, IACCEPTL is used to scale another multivariate normal proposal density. Thus, a mixture of two normal densities, with two different variance scales, are used as the proposal density. This serves as a pseudo t-distribution, but assuring radial symmetry as well as statistically independent samples, that may be useful when using with the Sobol method. This method has been suggested by Robert Leary, recommending IACCEPT=1.0 and IACCEPTL=0.01.

If IACCEPT is set to 0, and IACCEPTL is set to a value greater than 0, then a search for the best IACCEPTL for each subject is made, starting the testing at the IACCEPTL value given by the user, while IACCEPT is fixed to 1. The root.imp file will contain the final values selected for each subject, listing the two IACCEPT scale values, the first one being 1 for near the mode, and the second for the long tails.

3.16 GRDQ=0 (default) (NM74)

The gradient quick option, called GRDQ, allows THETAs that must be gradient assessed (such as those that are not mu-referenced) and SIGMAS to be more quickly evaluated by not evaluating the gradients for every one of the ISAMPLE random samples, but chooses a subset of the most important samples. This reduces the computational cost, since gradients of the objective function with respect to the THETAs require more objective function calls than is usually required when evaluating mu-referenced THETAs. If GRDQ>=1, then this is interpreted as the number of important samples to be used for theta gradient assessment per subject. If GRDQ<1, then GRDQ is interpreted as the fraction of ISAMPLEs to be used (GRDQ*ISAMPLE samples are used for theta gradient assessment). When GRDQ<0.0, then the number of samples used is ABS(GRDQ)*ISAMPLE/(Number of subjects with observations). When GRDQ=0 (default), then all ISAMPLE samples are used to evaluate the theta gradients. Some experience suggests that if GRDQ is too low (<30 samples), the quality of standard error assessments may deteriorate, so some trial and error may be needed to determine to what extent the GRDQ can be reduced. Suggestion of the GRDQ algorithm courtesy of Robert Leary.

During estimation, the iteration output for importance sampling contains the objective function, average effective number of samples (eff.), actual samples (smpl., which is usually ISAMPLE), and average fitness (fit.) of proposal (sampling) density to the posterior density. The average effective number of samples per total samples is essentially the average ratio (or weights) of posterior to proposal density at all sample positions, and NONMEM tries to adjust the scaling of the variance to the proposal density so that on average, eff./smpl. approximates the IACCEPT value. The average fitness is the average absolute difference between the posterior density and proposal density among the samples, and indicates to what extent that the proposal density matches the posterior density. Typically if the IACCEPT is not too low, then the closer the shape of the posterior density is to a normal density, the higher the fitness, with perfect fitness being 1 if IACCEPT is set to 1.

4 Monte Carlo Importance Sampling EM Assisted by Mode a Posteriori (MAP) estimation

Sometimes for highly dimensioned PK/PD problems with very rich data the importance sampling method does not advance the objective function well or even diverges. For this the IMPMAP method may be used. At each iteration, conditional modes and conditional first order variances are evaluated as in the ITS or FOCE method, not just on the first iteration as is done with IMP method. These are then used as parameters to the multivariate normal proposal density for the Monte Carlo importance sampling step. This method is implemented by:

$EST METHOD=IMPMAP INTERACTION

This is equivalent to

$EST METHOD=IMP INTERACTION MAPITER=1 MAPINTER=1

5 Stochastic Approximation Expectation Maximization (SAEM) Method

As in importance sampling, random samples are generated from normal distribution proposal densities. However, instead of always centered at the mean or mode of the posterior density, the proposal density is centered at the previous sample position. New samples are accepted with a certain probability. The variance of the proposal density is adjusted to maintain a certain average acceptance rate (IACCEPT). This method requires more elaborate sampling strategy, but is useful for highly non-normally distributed posterior densities, such as in the case of very sparse data (few data points per subject), or when there is categorical data.

In the first phase, called the burn-in or stochastic mode, SAEM evaluates an unbiased but highly stochastic approximation of individual parameters (semi integration, usually 2 samples per individual). Population parameters are updated from individual parameters by single iteration maximization steps that are very stable, and improves the objective function (usually in 300-5000 iterations). In the second mode, called the accumulation mode, individual parameter samples from previous iterations are averaged together, converging towards the true conditional individual parameter means and variances. The algorithm leads to population parameters converging towards the maximum of the exact likelihood.

The SAEM method is specified by

$EST METHOD=SAEM INTERACTION

Followed by one or more of the following options:

5.1 NBURN=2000

Maximum number of iterations in which to perform the stochastic phase of the SAEM method (default 1000). During this time, the advance of the parameters may be monitored by observing the results in file specified by the FILE parameter (described later in the Format of Output Files section), and the advance of the objective function (SAEMOBJ) at the console may be monitored. When all parameters or the SAEMOBJ do not appear to drift in a specific direction, but appear to bounce around in a stationary region, then it has sufficiently “burned” in. A termination test is available (described later), that will give a statistical assessment of the stationarity of objective function and parameters.

The objective function SAEMOBJ that is displayed during SAEM analysis is not valid for assessing minimization or for hypothesis testing. It is highly stochastic, and does not represent a marginal likelihood that is integrated over all possible eta, but rather, is the likelihood for a given set of etas.

5.2 NSAMPLE/NITER=1000

Sets maximum number of iterations in which to perform the non-stochastic/accumulation phase (default 1000).

5.3 ISAMPLE options

ISAMPLE=2     (defaults listed)
ISAMPLE_M1=2
ISAMPLE_M1B=2 (NM74)
ISAMPLE_M1A=0 (NM72)
ISAMPLE_M2=2
ISAMPLE_M3=2
IACCEPT=0.4

These are options for the MCMC Bayesian Metropolis-Hastings algorithm for individual parameters (ETAS) used by the SAEM and BAYES methods. For each ISAMPLE, SAEM performs ISAMPLE_M1 mode (or kernel) 1 iterations using the population means and variances as proposal density, followed by ISAMPLE_M1B mode 1B iterations using the individual conditional mean and individual conditional variance collected from previous iterations as proposal density, followed by ISAMPLE_M1A mode 1A iterations, testing model parameters from other subjects as possible values (by default this is not used, ISAMPLE_M1A=0), followed by ISAMPLE_M2 mode 2 iterations, using the present parameter vector position as mean, and a scaled variance of OMEGA as variance (Lavielle 2007). Next, ISAMPLE_M3 mode 3 iterations are performed, in which samples are generated for each parameter separately. The scaling is adjusted so that samples are accepted IACCEPT fraction of the time. The final sample for a given chain is then kept. The average of the isample parameter vectors and their variances are used in updating the population means and variances. Usually, these options need not be changed.

The ISAMPLE_M1A method of sampling has limited use to assist certain subjects to find good parameter values by borrowing from their neighbors, in case the neighbors had obtained good values while the present subject has difficulty finding good samples. This mode should generally not be used, and can be inaccurate if not all subjects share the same μ and Ω, such as in covariate modeling. Alternatively, use mode 1A sampling at the beginning of an SAEM analysis for a few burn in iterations, then continue with a complete SAEM analysis with mode 1A sampling turned off, with more burn in and accumulated sampling iterations, for example:

$EST METHOD=SAEM INTERACTION NBURN=500 NITER=0 ISAMPLE_M1A=2
$EST METHOD=SAEM INTERACTION NBURN=500 NITER=1000 ISAMPLE_M1A=0

As of NM75, If the option MCETA is set to greater than ISAMPLE_M1, then for the first iteration of SAEM, ISAMPLE_M1, ISAMPLE_M2, and ISAMPLE_M3 will be set to MCETA, to facilitate a robust initial search for reasonable parameter values.

5.4 ISAMPEND=n (NM73)

For SAEM, if ISAMPEND is specified as an upper integer value (usually 10), then NONMEM will perform a ISAMPLE preprocess to determine the best ISAMPLE value. For the ISAMPLE preprocessing the used entered ISAMPLE value must be at least 2. It will perform 200 iterations during the ISAMPLE preprocess, and the last 50 iterations will be used to obtain average conditional variance/OMEGA (eta shrinkage) for each subject. The largest etashrinkage fraction*10 is the ISAMPLE for that subject. Thus,

ISAMPLE=2 ISAMPEND=10

Will assess a best ISAMPLE for each subject. The ISAMPLE will not be higher than 10 or lower than 2.

5.5 ISCALE_MIN=1.0E-06 (defaults for SAEM, BAYES, NM72), ISCALE_MAX=1.0E+06 (NM72)

In MCMC sampling, the scale factor used to vary the size of the variance of the proposal density in order to meet the IACCEPT condition, is by default bounded by ISCALE_MIN of 1.0E-06, and ISCALE_MAX=1.0E+06. This should left alone for MCMC sampling, but on occasion there may be a reason to reduce the boundaries (perhaps to ISCALE_MIN=0.001, ISAMPLE_MAX=1000). After the SAEM estimation method, remember to revert these parameters back to default operation on the next $EST step:

ISCALE_MIN=-100 ISCALE_MAX=-100

The default operation is that NONMEM sets (ISCALE_MIN, ISCALE_MAX) to (0.1, 10) for importance sampling (as described earlier), and to (1.0E-06, 1.0E+06) for MCMC sampling.

5.6 NOCOV=[0,1] (NM73)

If covariance estimation is not desired for a particular estimation step, set NOCOV=1. It may be turned on again for the next estimation step with NOCOV=0. If NOCOV=1 is set for an FOCE/Laplace/FO method, this is equivalent to $COV NOFCOV setting. For ITS and IMP, covariance estimation can take some time for large problems, and you may wish to obtain only the objective function, such as in the case of $EST METHOD=IMP EONLY=1 after an SAEM estimation. NOCOV has no effect on BAYES analysis, as no extra time is required in assessing covariance for BAYES.

By default, standard error information for the classical methods (FO/FOCE/Laplace) will be given only if they are the last estimation method, even if NOCOV=0 for an intermediate estimation step. If NOCOV=1 for the FOCE/LAPLACE/FO method, and it is the last estimation step, then standard error assessment for it will be turned off.

5.7 DERCONT=[0,1] (NM73)

By the default value of the derivative continuity (DERCONT) is 0. When it equals 1, the partial derivative of the objective function with respect to THETAs will perform an additional test to determine if a backward difference assessment is more accurate than a forward difference assessment. The forward difference assessment can differ greatly from the backward difference assessment in cases of extreme discontinuity when varying certain THETAs by even just a small amount in the model results in a large change in objective function, (such as a viral model in which a very small change in the potency of an anti-viral agent results in widely varying time of return of viral load). This results in standard errors being poorly assessed for THETAs that do not have inter-subject variances associated with them. Setting DERCONT to 1 slows the analysis, but can provide more accurate assessments of SE in such models. The DERCONT works only for the Monte Carlo EM algorithms such as IMP and SAEM.

5.8 CONSTRAIN=1 (NM72)

A built-in simulated annealing algorithm has been put in place for NONMEM 7.2.0. Simulated annealing slows the rate of reduction of the elements of the OMEGA values during the burn-in phase of the SAEM method, allowing for a more global search of parameters. The subroutine CONSTRAINT performs this algorithm when the option CONSTRAIN is set to 1 or 5, where 1 is the default setting. This is by the constraint algorithm starting the OMEGAs at 1.5 times the initial values, and then controlling the rate at which the OMEGAs shrink during each iteration. CONTRAIN=2 or 6 performs simulated annealing on sigma parameters, CONSTRAIN=3 or 7 performs simulated annealing on both OMEGA and SIGMA parameters. CONSTRAIN=0 or 4, performs no simulated annealing on non-zero valued OMEGAS.

The user may modify the subroutine CONSTRAINT that performs the simulated annealing algorithm. The source code to the CONSTRAINT subroutine is available from the ..\source directory as constraint.f90, and the user may copy this to their run directory, and as convenient, to rename it. Then, specify OTHER=name_of_source.f90 in the $SUBROUTINE record, as shown in example 9.

As of NM73, when CONSTRAIN>=4, simulated annealing is also performed on diagonal elements of OMEGAS that are fixed to 0 to facilitate estimation of any associated THETAs. See I.57 $ANNEAL to facilitate EM search methods for this additional annealing technique. The subroutine CONSTRAINT may also be used to provide any kind of constraint pattern on any parameters.

As of NM72, the SAEM setting produces first order approximation standard errors, that is, MATRIX=S type, but not proper objective function for hypothesis testing.

5.9 PHITYPE=0 (NM74)

By default, after an estimation is performed, the phi(), conditional means of the individual parameters, and their variances, are reported in the root.phi file, where root is the root name of the control stream file. If you wish to have conditional mean etas reported, set PHITYPE=1. See root.phi in section I.65 $EST: Additional Output Files Produced for more information.

5.10 MAPITERS=0 (default) (NM75)

By default, no MAP estimation is performed with SAEM. To get good individual parameter values near the mode of the posterior density for the first iteration of SAEM, you can set MAPITERS=1. Alternatively, you can insert the record:

$EST METHOD=ITS NITER=0

Followed by

$EST METHOD=SAEM …

5.11 Monolix / NONMEM SAEM parameter mapping

The mapping of parameters between Monolix and NONMEM SAEM is as follows:

Monolix	NONMEM SAEM
Number of Chains	ISAMPLE
K0	CONSTRAINT subroutine may be user modified to provide any constraining pattern on any population parameters
K1	NBURN
K2	NITER
Auto K1	CTYPE=1,2,3
Population Parameter settings menu:
rho	IACCEPT
m1	ISAMPLE_M1
m2	ISAMPLE_M1A
m3	ISAMPLE_M2
m4	ISAMPLE_M3
No simulated annealing	CONSTRAIN=0
Simulated Annealing	CONSTRAIN=1,2,3 (User may also define algorithm)
SEED	SEED

5.12 Obtaining the Objective Function for Hypothesis Testing After an SAEM Analysis

After the analysis, suitable objective functions for hypothesis testing and second order standard errors can be obtained by importance sampling at the final population parameter values. Thus, one could issue this sequence of commands:

$EST METHOD=SAEM INTERACTION NBURN=2000 NITER=1000
$EST METHOD=IMP EONLY=1 ISAMPLE=1000 NITER=5

Here, after SAEM is performed, importance sampling, with MAP estimation done on its first iteration, is performed, but without updating the main population parameters. Sometimes the MAP estimation is problematic, and/or, the user wishes to use the SAEM’s last conditional mean and variances as the parameters to the importance sampler’s sampling density for the first iteration, so one may try:

$EST METHOD=SAEM INTERACTION NBURN=2000 NITER=1000
$EST METHOD=IMP EONLY=1 ISAMPLE=1000 NITER=5 MAPITER=0

For very large dimensioned problems (many OMEGAs), the IMP evaluated objective function can have a lot of stochastic variability (more than plus or minus 10 units), or continually increase with each iteration even though the population parameters are kept fixed. One way to reduce this volatility is to use IMPMAP instead of IMP, if the MAP estimation is not an issue:

$EST METHOD=IMPMAP EONLY=1 ISAMPLE=1000 NITER=5 MAPITER=0

Another way is to increase the ISAMPLE to 3000:

$EST METHOD=IMP EONLY=1 ISAMPLE=3000 NITER=5 MAPITER=0

and sometimes, using the combination of IMPMAP with ISAMPLE=3000 is needed. Using IMPMAP or increasing ISAMPLE do increase computation time, and it is a choice of which is more efficient.

As of NM74, another choice is to set EONLY=2, and the Monte Carlo variability of the objective function can be significantly reduced if, like the population parameters, the conditional means and variances from the SAEM estimation are also not updated after each IMP iteration:

$EST METHOD=SAEM INTERACTION NBURN=2000 NITER=1000
$EST METHOD=IMP EONLY=2 ISAMPLE=1000 NITER=5 MAPITER=0

Or, have the first iteration evaluate conditional modes and conditional MAP variances, and then use them for the subsequent iterations, by setting EONLY=3 and MAPITER=1:

$EST METHOD=SAEM INTERACTION NBURN=2000 NITER=1000
$EST METHOD=IMP EONLY=3 ISAMPLE=1000 NITER=5 MAPITER=1

Another set of commands for SAEM is the following, which begins with a short iterative two stage run to provide good initial eta estimates for each subject, followed by the SAEM analysis, which uses these initial eta estimates as a starting point for its Markov Chain Monte Carlo scan of each subject’s conditional (posterior) density, followed by objective function evaluation:

$EST METHOD=ITS INTERACTION NITER=5
$EST METHOD=SAEM NBURN=1000 ISAMPLE=2 NITER=1000
$EST METHOD=IMP EONLY=2 ISAMPLE=1000 NITER=5 MAPITER=0

Values of NBURN, NITER, and ISAMPLE may be changed as needed.

If you want conditional mean values (values listed in root.phi) evaluated by MCMC sampling used in the SAEM method, but at a constant set of the final fixed parameters, then you could invoke EONLY=1 with the SAEM method as well:

$EST METHOD=ITS INTERACTION NITER=5
$EST METHOD=SAEM NBURN=1000 ISAMPLE=2 NITER=1000
$EST METHOD=SAEM EONLY=1 NBURN=200 ISAMPLE=2 NITER=1000
$EST METHOD=IMP EONLY=2 ISAMPLE=1000 NITER=5 MAPITER=0

5.13 ETASAMPLES=0 (default) (nm74)

As of NM74, one can obtain random samples of individual etas, while keeping population parameters fixed, and uses these for covariate and model diagnostics (see (Lavielle and Ribba 2016)). To do this, perform an SAEM analysis (or BAYES) after the primary estimation, keeping the population parameters fixed (ETYPE=1), only performing a BURN period so samples do not get accumulatively averaged, and set ISAMPLE to 10 or higher to collect suffient samples per subject. For example, note the $EST line in bold, after the usual SAEM analysis:

$SIZES ISAMPLEMAX=250
…
$EST METHOD=SAEM NBURN=500 NITER=1000 MASSRESET=1 ISAMPLE=2
         PRINT=20 NOPRIOR=1 RANMETHOD=P
$EST METHOD=SAEM NBURN=0 NITER=0 MASSRESET=0 ETASAMPLES=1 ISAMPLE=200 EONLY=1
$EST METHOD=IMP NITER=5 MASSRESET=1 ETASAMPLES=0 ISAMPLE=1000
         EONLY=1 PRINT=1 MAPITER=0

ETASAMPLES=1 causes individual ISAMPLE random eta samples per subject, to be written to root.ets, where root is the root name of the control stream file. MASSRESET is set to 1 to initialize the internal burn-in coefficients, and collect that information during the population parameter estimation. Then, have the next SAEM step accumulate eta samples (ETASAMPLES=1), and set MASSRESET=0 so the internal burn-in coefficients are note reset, but use those obtained during the previous SAEM step. After this, you could request the IMP step (IMP needs to use its own burn-in or shaping coefficients, so re-initialize with MASSRESET=1).

Suppose your result was from a previous analysis, called precevample.ctl. You could do the following for the present problem:

; Obtain final result of previous example from its ext file
$EST METHOD=CHAIN NSAMPLE=0 TBLN=1 ISAMPLE=-1000000000 FILE=prevexample.ext
; Have a burn-in process with just 2 chains (ISAMPLE=2), to create appropriate
; burn-in coefficients
$EST METHOD=SAEM NBURN=0 NITER=0 MASSRESET=1 ISAMPLE=2 EONLY=1 RANMETHOD=P
; Use burned-in coefficients (MASSRESET=0) and collect many samples, and store
in .ets file.
$EST METHOD=SAEM NBURN=0 NITER=0 MASSRESET=0 ETASAMPLES=1 ISAMPLE=200 EONLY=1

The ETASAMPLES option of collection of individual parameter samples is comparable to the BIONLY option for BAYES method (BIONLY=0 (default) (NM75) in section I.41 Full Markov Chain Monte Carlo (MCMC) Bayesian Analysis Method).

5.14 SAEM_EXTRA_REQUEST and SAEM_EXTRA Reserved Variables (nm760)

Occasionally one wishes to obtain a distribution of any individual parameter defined in the control stream code, after an SAEM analysis is completed, and not just the phis/etas. As of nm760, a file can be generated during an SAEM ETASAMPLES=1 individual parameter collection process, that collects the posterior distribution of desired parameters, such as individual weighted residuals, as shown here (shapiro_iwr2.ctl):

$SIZES ISAMPLEMAX=1100
;Model Desc: Two compartment Model, Using ADVAN3, TRANS4
;Project Name: example1
;Project ID: gm00-001

$PROB RUN# shapiro_iwr2 (from samp5l)
$ABBR DECLARE INTEGER FIRST_WRITE
$INPUT C SET ID JID TIME DV=CONC AMT=DOSE RATE EVID MDV CMT
            CLX V1X QX V2X SDIX SDSX
$DATA shapiro_iwr2.csv IGNORE=C
$SUBROUTINES ADVAN3 TRANS4

$PK
include nonmem_reserved_general
SAEM_EXTRA_REQUEST=1 ; Request NONMEM to trigger SAEM-EXTRA events
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
CL=DEXP(MU_1+ETA(1))
V1=DEXP(MU_2+ETA(2))
Q=DEXP(MU_3+ETA(3))
V2=DEXP(MU_4+ETA(4))
S1=V1

$ERROR
include nonmem_reserved_general
SAEM_EXTRA_REQUEST=1
IPRE=F
IWRES=(DV-IPRE)/IPRE/SQRT(SIGMA(1,1))
Y = F + F*EPS(1)
; When SAEM_extra=1, then this particular set of individual
; parameters were "accepted" So you may record them if you wish
IF(SAEM_EXTRA==1 .AND. EVID==0 .AND. FIRST_WRITE==0) THEN
IF(PNM_RUN_MODE==PNM_SINGLE) THEN ; best if a pnm file with -awnf is used.
" OPEN(UNIT=52,FILE='shapiro_iwr2.iwr')
ELSE ;If parallel, need to have different file name for each node
" OPEN(UNIT=52,FILE='shapiro_iwr2_'//trim(tfi(pnm_node_number))//'.iwr')
ENDIF
IF(PNM_NODE_NUMBER==1) THEN ; Only one file should have header, for easy
merging later
" WRITE(52,96) 'SAMPLE  ','ID  ','TIME                    ','IWRES        '
" 96 FORMAT(A12,1X,A14,2(1X,A12))
ENDIF
FIRST_WRITE=1
ENDIF

IF(SAEM_EXTRA==1 .AND. EVID==0 ) THEN
" WRITE(52,97) SAEM_SAMPLE,ID,TIME,IWRES
" 97 FORMAT(I12,1X,F14.0,2(1X,1PG12.5))
ENDIF

; Initial values of THETA
$THETA
(2.0) ;[LN(CL)]
(2.0) ;[LN(V1)]
(2.0) ;[LN(Q)]
(2.0) ;[LN(V2)]
;INITIAL values of OMEGA

$OMEGA BLOCK(4)
0.15 ;[P]
0.01 ;[F]
0.15 ;[P]
0.01 ;[F]
0.01 ;[F]
0.15 ;[P]
0.01 ;[F]
0.01 ;[F]
0.01 ;[F]
0.15 ;[P]
;Initial value of SIGMA

$SIGMA
(0.6 ) ;[P]

$EST METHOD=SAEM NBURN=500 NITER=300 MASSRESET=1 PRINT=20 CTYPE=3
$EST METHOD=SAEM NBURN=0 NITER=0 MASSRESET=0 ETASAMPLES=1 ISAMPLE=1000
EONLY=1

You will notice that SAEM_EXTRA and SAEM_EXTRA REQUEST parameters are available (as of nm760). NONMEM will set SAEM_EXTRA to 1 only when ETASAMPLES=1 is set, so you will only get samples outputted during the special period when you want to collect individual parameters with population parameters fixed, so it is not necessary to test for ETASAMPLES=1 in the control stream code.

6 Full Markov Chain Monte Carlo (MCMC) Bayesian Analysis Method

The goal of the MCMC Bayesian analysis [12,13] is not to obtain the most likely THETAs, SIGMAs, and OMEGAs, but to obtain a large sample set of probable population parameters, usually 10000-30000. The samples are not statistically independent, but when analysis is properly performed, they are uncorrelated overall. Various summary statistics of the population parameters may then be obtained, such as means, standard deviations, and even confidence (or credible) ranges. The mean population parameter estimates and their variances are evaluated with considerable stability. Maximum likelihood parameters are not obtained, but with problems of sufficient data, these sample mean parameters are similar to maximum likelihood values, and the standard deviations of the samples are similar to standard errors obtained with maximum likelihood methods. A maximum likelihood objective function is also not obtained, but, a distribution of joint probability densities is obtained, from which 95% confidence bounds (assuming a type I error of 0.05 is desired) can be constructed and tested for overlap with those of alternative models.

As with the SAEM method, there are two phases to the BAYES analysis. The first phase is the burn-in mode, during which population parameters and likelihood may change in a very directional manner with each iteration, and which should not be used for obtaining statistical summaries. The second phase is the stationary distribution phase, during which the likelihood and parameters tend to vary randomly with each iteration, without changing on average. It is these samples that are used to obtain summary statistics.

The Bayesian method is specified by

$EST METHOD=BAYES INTERACTION

Followed by one or more of the following parameter options:

6.1 NBURN=4000

Maximum number of iterations in which to perform the burn-in phase of the MCMC Bayesian method (default 4000). During this time, the advance of the parameters may be monitored by observing the results in file specified by the FILE parameter, and/or the objective function displayed at the console. The objective function progress is also written in OFV.TXT, and the report file. Full sets of population parameters and likelihood functions are also written in the file specified with the FILE= option. When all parameters and objective function do not appear to drift in a specific direction, but appear to bounce around in a stationary region, then it has sufficiently “burned” in. A termination test may be implemented to perform a statistical assessment of stationarity for the objective function and parameters. As mentioned earlier, the objective function (MCMCOBJ) that is displayed during BAYES analysis is not valid for assessing minimization or for hypothesis testing in the usual manner. It does not represent a likelihood that is integrated over all possible eta (marginal density), but the likelihood at a given set of etas.

6.2 NSAMPLE/NITER=10000

Sets number of iterations in which to perform the stationary distribution for the BAYES analysis (default 10000).

6.3 ISAMPLE options

ISAMPLE_M1=2 (defaults listed)
ISAMPLE_M1A=0 (NM72)
ISAMPLE_M1B=2 (NM74)
ISAMPLE_M2=2
ISAMPLE_M3=2
IACCEPT=0.4

These are options for the MCMC Bayesian Metropolis-Hastings algorithm for individual parameters (ETAS) used by the SAEM and BAYES methods. For Bayesian analysis, the MCMC algorithm performs ISAMPLE_M1 mode 1 iterations using the population means and variances as proposal density, followed by ISAMPLE_M1B mode 1B iterations using the individual conditional mean and individual conditional variance collected from previous iterations as proposal density, followed by ISAMPLE_M1A mode 1A iterations, testing model parameters from other subjects as possible values (by default this is not used, ISAMPLE_M1A=0), followed by ISAMPLE_M2 mode 2 iterations, using the present parameter vector position as mean, and a scaled variance of OMEGA as variance (Lavielle 2007). Next, ISAMPLE_M3 mode 3 iterations are performed, in which samples are generated for each parameter separately. The scaling is adjusted so that samples are accepted IACCEPT fraction of the time. The final sample for a given chain is then kept. Usually, these options need not be changed. There is only one chain of samples produced for a given NONMEM run (ISAMPLE is not used for MCMC, only for SAEM). If you would like additional chains, then create separate control stream files with different starting seed numbers.

As of NM75, if the option MCETA is set to greater than ISAMPLE_M1, then for the first iteration of SAEM, ISAMPLE_M1, ISAMPLE_M2, and ISAMPLE_M3 will be set to MCETA, to facilitate a robust initial search for reasonable parameter values.

6.4 ISCALE_MIN=1.0E-06 (defaults for SAEM, BAYES, NM72), ISCALE_MAX=1.0E+06 (NM72)

In MCMC sampling, the scale factor used to vary the size of the variance of the proposal density in order to meet the IACCEPT condition, is by default bounded by ISCALE_MIN of 1.0E-06, and ISCALE_MAX=1.0E+06. This should left alone for MCMC sampling, but on occasion there may be a reason to reduce the boundaries (perhaps to ISCALE_MIN=0.001, ISAMPLE_MAX=1000). After the SAEM estimation method, remember to revert these parameters back to default operation on the next $EST step:

ISCALE_MIN=-100 ISCALE_MAX=-100

The default operation is that NONMEM sets (ISCALE_MIN, ISCALE_MAX) to (0.01, 100) for importance sampling (as described earlier), and to (1.0E-06, 1.0E+06) for MCMC sampling.

6.5 IKAPPA=1 (default)

The individual parameters are averaged using a weight N⁻ᴵᴷᴬᴾᴾᴬ for the Nth iteration (so a simple average with each iteration’s value equally weighted would be IKAPPA=1), in obtaining the mean and variance-covariance for the ISAMPLE_M1B mode. A value of 0.75 can sometimes provide an improved decorrelation efficiency when performing standard Bayesian analysis.

6.6 PSAMPLE options

PSAMPLE_M1=1 (defaults listed)
PSAMPLE_M2=-1
PSAMPLE_M3=1
PACCEPT=0.5

These are the options for the MCMC Metropolis-Hastings algorithm. These options only have meaning for population parameters (theta/sigma) that are not Gibbs sampled. Normally NONMEM determines whether THETA and SIGMA parameters are Gibbs sampled or not, based on the model setup (see MU_ Referencing section below). For each iteration, a vector of THETAs/SIGMAs are generated using a multivariate normal proposal density that has mean/variances based on the previous samples, done PSAMPLE_M1 times. Next, a vector of parameters are generated using a multivariate normal proposal density with mean at the present parameter position, and variance scaled to have samples accepted with PACCEPT frequency. This is done PSAMPLE_M2 times (if PSAMPLE_M2<0, then program performs this as many times as there are M-H parameters). Finally, each parameter is individually sampled PSAMPLE_M3 times. The final accepted parameter vector is kept. Usually these options do not need to be changed from their default values, listed above.

6.7 PSCALE_MIN=0.01 (NM73), PSCALE_MAX=1000 (NM73)

In MCMC sampling, the scale factor used to vary the size of the variance of the proposal density population parameters (theta/sigma) that are not Gibbs sampled, in order to meet the PACCEPT condition, is by default bounded by PSCALE_MIN of 0.01, and PSCALE_MAX=1000. This should left alone for MCMC sampling, but on occasion there may be a reason to expand the boundaries (perhaps to PSCALE_MIN=1.0e-06, PSCALS_MAX=1.0E+06).

6.8 OSAMPLE options

OSAMPLE_M1=-1 (defaults listed)
OSAMPLE_M2=-1
OSAMPLE_M3=-1
OACCEPT=0.5

These are the options for the MCMC Metropolis-Hastings algorithm for OMEGA sampling. If OSAMPLE_M1<0 (default), then the OMEGAs are Gibbs sampled using the appropriate inverse Wishart proposal density, and the other options (OSAMPLE_M2 and OACCEPT) are not relevant. Otherwise, for each iteration, a matrix of OMEGAs are generated using an inverse Wishart proposal density that has variance based on the previous samples, done OSAMPLE_M1 times. Next, a matrix of OMEGAS are generated using a Wishart proposal density at the present OMEGA values postion, and degrees of freedom (dispersion factor for variances) scaled to have samples accepted with OACCEPT frequency. This is done OSAMPLE_M2 times (if OSAMPLE_M2<0, then program performs this as many times as there are non-fixed omega elements). Then, individual cholesky elements of OMEGA are varied, each OSAMPLE_M3 times (if OSAMPLE_M3<0, then program performs this as many times as there are non-fixed omega elements). The final OMEGA matrix is kept. Usually these options do not need to be changed from their default values, listed above.

6.9 NOPRIOR=[0,1]

If prior information was specified using the $PRIOR statement (available since NM 6, release 2.0, and described in the Help manual ..\guides\VIII.pdf: use only NWPRI option for the new $EST methods), then normally the analysis is set up for three stage hierarchical analysis. By default NOPRIOR=0, and this prior information will be used. However, if NOPRIOR=1, then for the particular estimation, the prior information is not included in the analysis. This is useful if you want to not use prior information during a maximization (METHOD=IMP, CONDITIONAL, IMPMAP, SAEM, or ITS), but then use it for the Bayesian analysis (METHOD=BAYES).

As of NM73, when NOPRIOR=1 is set, the estimation will not use TNPRI prior information (TNPRI should only be used with FO/FOCE/Laplace estimations). In previous versions of NONMEM, NOPRIOR=1 did not act on TNPRI priors.

6.10 MAPITERS=0 (default) (NM75)

As with SAEM, by default, no MAP estimation is performed with BAYES. To get good individual parameter values near the mode of the posterior density for the first iteration of BAYES, you can set MAPITERS=1. Alternatively, you can insert the record:

$EST METHOD=ITS NITER=0

Followed by

$EST METHOD=BAYES …

6.11 THIN=1 (default, NM74)

As of nm74, the Bayesian records retained in the raw output file may be adjusted by every THINth iteration. So, if THIN=10, then every 10th iteration is recorded in the raw output file. The PRINT option controls only the iterations printed to the console and NONMEM report file.

6.12 BAYES_PHI_STORE=0 (default) (NM75)

As of nm75, samples of phi/eta are collected at each BAYES iteration, and summarized to provide mean phi and phc() in the root.phi table, as described above. By default, the individual phi values from each iteration are not stored. However, if you set

$EST … BAYES_PHI_STORE=1

then phi and eta values from each BAYES iteration will be stored in root.iph. For non-mixture problems, only records of SUBP=0 are recorded, as there are no sub-population divisions. For mixture problems, the SUBP=0 records contain the composite phis and etas (the average of these across all non-negative iterations are in the root.phi table), and the SUBP>0 records contain the phis and etas appropriate to each sub-population SUBP (the average of these across all non-negative iterations are in the root.phm table). The root.iph file can become quite large, so it should be used only on the final analysis. This root.iph file can be read by the utility program NEFFI (see I.102 NEFF and NEFFI Utility Programs (NM74)).

6.13 BIONLY=0 (default) (NM75)

BIONLY stands for Bayesian individual parameters only, and when set to 1, will create new samples of individual parameters only, but will keep the population parameters fixed. This is equivalent to EONLY=1 used for SAEM and IMP estimation. Of course, for this to be useful, either write statements capturing the individual parameters samples in the control stream must be entered (see example 8 near end of this manual), or the BAYES_PHI_STORE option must be set to 1, to collect individual samples in the .iph file.

Here is an example showing its use:

$EST METHOD=BAYES INTERACTION NBURN=200
         NITER=1000 PRINT=100 NOPRIOR=1 CINTERVAL=100 BIONLY=1
         BAYES_PHI_STORE=1

The BIONLY option of BAYES, with collection of individual samples in the .iph file, is comparable to the ETASAMPLES option of SAEM, with collection of individual samples in the .ets file (ETASAMPLES=0 (default) (nm74) in section I.40 Stochastic Approximation Expectation Maximization (SAEM) Method).

7 No U-Turn Sampling (NUTS) Markov Chain Monte Carlo (MCMC) Bayesian Analysis Method (NM74)

The No U-Turn sampling algorithm was developed by Hoffman, Gelman, and others of the STAN Development Team ([19,20]). The algorithm developed in NONMEM is a limited form of the No-U-Turn MCMC Bayesian environment that is available in STAN, focused on analyzing population PK/PD models, with normally distributed THETA priors (or t-distributed theta priors, see TTDF), and Wishart/Gamma distributed OMEGA and SIGMA priors, or using the LKJ correlation prior (see OLKJDF and SLKJDF options) for OMEGAS and SIGMAS. The algorithm in NONMEM is based on references [19,20], and helpful suggestions by Bob Carpenter, Andrew Gelman, Matt Hoffman, Michael Betancourt, and Sebastian Weber.

Typical Bayesian algorithms search for individual parameters (phis or etas) and population parameters (THETAs, OMEGAs, and SIGMAs) in separate stages. While this provides for rapid generation of an MCMC sample, the samples can be heavily correlated, especially if there are high correlations in the OMEGA matrix. The NUTS sampler uses a directed search using partial derivatives and scaling techniques using posterior density knowledge from previous samples to reduce the correlation of the parameters from one iteration to the next. While each iteration takes longer to generate with NUTS, the samples may be 10-100 times decorrelated relative to a standard MCMC sampling. Thus, 10000 samples with NUTS may be worth 100000 samples with traditional MCMC algorithms. See the NEFF (Number of EFFective samples) utility for analyzing the quality of an MCMC run, I.102 NEFF and NEFFI Utility Programs (NM74).

Because the NUTS algorithm relies on derivatives, it is best if analytical derivatives are created for each of the estimated parameters. For OMEGAs and SIGMAs, these are done automatically, but for THETAs, analytical derivatives are created only if they are MU referenced (see MU Reference section below). So MU reference all THETAs that are to be estimated. It is okay to set their OMEGA to 0, the analytical derivative will still be utilized.

The easiest way to use the NUTS algorithm in NONMEM is to use the AUTO option (see section I.45 Some General Options and Notes Regarding EM and Monte Carlo Methods for more details on the AUTO feature, AUTO=0 (default) (NM73)). For example, stanrb42.ctl uses AUTO=1:

$EST METHOD=NUTS AUTO=1 PRINT=20

The example stanrb10.ctl shows another setup for using NUTS:

$EST METHOD=ITS NITER=0 file=stanrb10_its.ext
$EST METHOD=NUTS NBURN=1000 NITER=2000 PRINT=20 OLKJDF=3.0 file=stanrb10.ext

The OLKJDF option specifies the degrees of freedom to the LKJ decorrelation density for the OMEGA prior. The OLKJDF should be set to a value greater than 0 to use the LKJ decorrelation prior for OMEGAs, as recommended by the STAN group. Also see comment below about using LKJ decorrelation versus inverse Wishart prior. Experience has shown that a low non- or weakly-informative OLKJDF should be at least 2, and no greater than the number of the Omega diagonals. If OLKJDF is set to 1, then there may be considerable pull towards high correlations which reduces efficiency of sampling. A single iteration of ITS is helpful to center the initial etas at their modes, as a good facilitator for initiating the NUTS run that follows it (NITER=0, so population parameters are not advanced).

An alternative method is to use the traditional BAYES method to rapidly generate samples for an initial mass matrix, which can then be passed on to the NUTS algorithm (example stanrb9):

$EST METHOD=BAYES NBURN=2000 NITER=2000 PRINT=50 MASSRESET=1
file=stanrb9_bayes.ext KAPPA=0.75
$EST METHOD=NUTS NBURN=500 NITER=2000 PRINT=20 file=stanrb9.ext
         OLKJDF=3.0 MASSRESET=0 KAPPA=1.0 MADAPT=250

Notice that the MASSRESET is set to 1 to initialize the mass matrix accumulator at the BAYES step, and then set MASSRESET=0 at the NUTS step, so that the mass matrix is not re-initialized, but rather, carried over, from the BAYES step. Also, you may wish to set KAPPA=0.75 during BAYES so that the accumulated mass matrix favors values collected during the latter portion of the BAYES analysis (for the NUTS step itself, KAPPA should be set back to 1). This technique can sometimes make the burn-in (warm-up) period for NUTS execute faster, and/or improve the de-correlation. MADAPT is set to ½ of NBURN, during which the mass matrix is continuously updated every NUTS_BASE iterations (which in this case by default is 0.025*NBURN=25).

Other stanrb* examples in the examples directory show the various ways in which the problem may be analyzed.

Several versions of a differential equation problem, example6, are also in the examples directory, example6hmt*, to be compared with example6classic2 and example6classico3. Comparing example6hmto26.ctl using pure NUTS algorithm with example6hmto19 which uses a pre-warmup from a previous BAYES estimation, you can see that the pre-warmup can reduce time of computation.

The following is a list of options the user may play with.

7.1 METHOD=NUTS

This sets the Bayesian analysis for the NUTS method.

7.2 NUTS_MASS=B (default)

By default, NUTS_MASS=B, the NUTS method is to use a block diagonal mass matrix for scaling its search. The THETAs/SIGMAs/OMEGAs and their correlations will be scaled with one block matrix, and parameter sets of each individual will have their own block matrix correlations. Correlations between THETAs and individual parameters will not be accounted for. Most efficient.

F: Full mass matrix. A full correlation matrix between THETAs, OMEGAs, SIGMAs, and individual parameters among all individuals will be accounted for. Computationally very expensive. When using METHOD=BAYES as a preparation for NUTS, as shown in example stanrb_177.ctl, make sure you set the NUTS_MASS (if you will be using something different from the default B value) at the METHOD=BAYES record, as well as MASSRESET=1, as this is required to set the appropriate memory allocation, and store the posterior variance-covariance (mass matrix) information that the NUTS algorithm will then use.

D: Diagonal mass matrix. No correlations between parameters will be considered.

BD: Block mass matrix covering THETAs, SIGMAs, OMEGAs, and diagonal mass matrix on individual parameters.

DB: Diagonal mass matrix covering THETAs, SIGMAs, OMEGAs, and block mass matrix on individual parameters.

BBD: THETAs and SIGMAs are blocked together, OMEGAs are in their own block, and individual parameters are diagonal.

BBB: THETAs and SIGMAs are blocked together, OMEGAs are in their own block, and individual parameters are blocked within each subject.

The mass matrix is generated by accumulating previous samples of parameters and obtaining their variance-covariance, so that the NUTS algorithm performs an efficient search in the domain of the empirical posterior density. It is best to acquire this mass matrix by first performing a couple of thousand iterations using standard Gibbs Bayesian analysis, followed by a NUTS process.

7.3 MASSRESET=-1 (default)

By default mass matrix information accumulation is turned off. However, MASSRESET=1 should be set when performing a Gibbs/MH Baysian analysis to initialize accumulation of mass matrix information, followed by the NUTS algorithm, with MASSRESET=0 set, so the mass matrix accumulator information is not reset to 0, but adds to the information acquired during the previous Gibbs/MH process. IF you have a BAYES record setup previous to the NUTS record, make sure to set the NUTS_MASS value at the BAYES record first, so it accumulates the correct type of mass matrix information, and allocates the appropriate memory for its storage. See the description of MADAPT=-1 (default), below.

7.4 NUTS_MAXDEPTH=10 (default)

This sets the maximum number of total branchings to try in the NUTS algorithm in the search for the next decorrelated sample. If many messages are received of reaching the maximum buildtree level, increase NUTS_MAXDEPTH.

7.5 MADAPT=-1 (default)

If MADAPT/=-1 (also called PMADAPT), the mass matrix is accumulative average updated throughout the NUTS analysis during the burn-in period, the update occurring every NUTS_BASE iterations, for the first MADAPT iterations for the parameters, then changes no further after that. For example, if NBURN=300, NUTS_BASE=20, and MADAPT=200, then during the first 20 burn-in iterations, an identity mass matrix is used because there is no past history of samples collected yet, then during the next 20 burn-in iterations the mass matrix derived from the first twenty burn-in iterations is used, then for following 20 burn-in iterations, the mass matrix from the first 40 burn-in iterations are used, etc. If KAPPA/=1, the samples of each are not evenly weighted, but are weighted according to the description of KAPPA below. If MASSRESET=0 because the previous $EST was a METHOD=BAYES … MASSRESET=1 estimation from which an initial mass matrix was derived, then the first 20 burn-in iterations use this BAYES pre-prepared mass matrix, the next 20 burn-in iterations use a mass matrix derived from an average of the pre-prepared BAYES mass matrix (weighted as 2number of population parameters to be estimated) plus the first 20 burn-in iterations (weighted as 20), the next 20 burn-in iterations use a mass matrix derived from an average of the pre-prepared BAYES mass matrix (weighted as 2number of population parameters to be estimated) plus the first 40 burn-in iterations (weighted as 40), etc. This weighting is modified according to the KAPPA value used.

In the following code portion (from examples\nuts_example6h.ctl), a brief iterative 2 stage is performed to center the parameters quickly in the stationary zone, followed by a BAYES estimation to quickly create a starting mass matrix, followed by the NUTS estimation, that uses the BAYES derived mass matrix as a starting point:

$EST METHOD=ITS … NITER=15 …
$EST METHOD=bayes … NBURN=200 NITER=0 … MASSRESET=1
$EST METHOD=NUTS … NBURN=100 NUTS_BASE=20 … NITER=200 MASSRESET=0
         PMADAPT=100

In the above setting, the mass matrix is updated every 20 (NUTS_BASE) iterations during the entire 100 NBURN iterations (since PMADAPT=NBURN).

In another example, (from examples\nuts_example6g.ctl), the mass matrix from the BAYES estimation is used directly for the stationary sampling phase of the NUTS estimation, as there is no burn-in period for the NUTS estimation (NBURN=0):

$EST METHOD=ITS … NITER=15 …
$EST METHOD=bayes … NBURN=200 NITER=0 … MASSRESET=1
$EST METHOD=NUTS … NBURN=0 … NITER=200 MASSRESET=0 PMADAPT=100

If MADAPT=-1, then the STAN method of warmup and mass matrix accumulation is used (according to the STAN manual (Stan Development Team 2015), Optimization Algorithms section). When using MADAPT=-1, the tuning options, NUTS_INIT, NUTS_BASE, and NUTS_TERM are useful. Use of the AUTO option allows for easy setup of the STAN method of matrix accumulation without having to deal with every one of these options. An example of this is examples\nuts_example6f.ctl:

$EST METHOD=NUTS AUTO=1 NITER=200 NUTS_INIT=25 …

The AUTO=1 option (see section I.45 Some General Options and Notes Regarding EM and Monte Carlo Methods for more details on the AUTO feature, AUTO=0 (default) (NM73)) automatically sets the various options to appropriate default values (including setting MADAPT=-1) for STAN type mass matrix development, unless over-written by the user, as is done above for NUTS_INIT (default is 75), to reduce computation time. Using the STAN burn-in method for creating the mass matrix takes longer than using the BAYES pre-analysis BAYES as done in nuts_example6g.

7.6 KAPPA=1 (default)

The parameters are averaged using a weight N⁻ᴷᴬᴾᴾᴬ (also called PKAPPA) for the Nth iteration (so a simple average with each iteration’s value equally weighted would be KAPPA=1), in obtaining the mass matrix. A value of 0.75 gives the best results when preparing a mass matrix during the BAYES step, in anticipation of the NUTS step.

7.7 NUTS_GAMMA=0.05 (default)

Gamma factor in the NUTS algorithm. Should not be changed, and NUTS experts (the Stan developers) recommend 0.05.

7.8 NUTS_DELTA=0.8 (default)

This is essentially the sample acceptance rate for the NUTS sampling process, equivalent to PACCEPT in standard MH sampling. NUTS experts recommend 0.8.

7.9 TTDF=0 (default)

TTDF stands for Theta t-density degrees of freedom. When 0, the usual normal density prior is used as a prior density for THETAs. When TTDF>0, then a t-distributed prior is used. TTDF may be set >0 when using METHOD=BAYES as well, but THETAs will then be M-H sampled using the PSAMPLE_M1, PSAMPLE_M2, and PSAMPLE_M3 settings. TTDF may be a real number. See record $TTDF below to specify degrees of freedom for each theta.

7.10 TPU=0 (default, NM75)

This is a $EST option that specifies, if >0, to use the THETA_PRIORU routine in ..\source\THETA_PRIORU.f90 as a template to provide the desired probability density function. To do this, you can make a copy of THETA_PRIORU.f90, call the file USERPRIORT.f90, for example, modify its THETA_PRIORU subroutine, then specify this file as an OTHER routine:

$SUBR … OTHER=USERPRIORT.f90

7.11 OLKJDF=0 (default)

OLKJDF stands for Omega LKJ density degrees of freedom. When 0, the usual inverse Wishart prior is used for OMEGAs. When OLKJDF>0, then the LKJ density is used as the prior, with OLKJDF degrees of freedom for all omega blocks. In addition, only diagonal elements of the OMEGA prior are used, assuming a density dependent on the OVARF value. OLKJDF may be set >0 when using METHOD=BAYES as well, but OMEGAs will then be M-H sampled using the OSAMPLE_M1, OSAMPLE_M2, and OSAMPLE_M3 settings. See record $OLKJDF below to specify LKJ correlation degrees of freedom for each omega block.

7.12 OVARF=1 (default)

OVARF is the weight factor to STD prior to the log sqrt OMEGA diagonal elements, the normal density of the log square root of OMEGA centered about log square root of Omega prior, and scaled with OVARF (see below). That is,

\[\log(\sqrt{\Omega_i}) \sim \text{Normal}(\log(\sqrt{\Omega_{\text{prior},i}}), 1/\text{OVARF}).\]

If OVARF<0, then a half-t-distribution of degrees of ABS(OVARF) is used as the prior to the sqrt of OMEGA diagonal elements. Use OVARF=-1 for the half-Cauchy distribution.

7.13 SLKJDF=0 (default)

SLKJDF stands for Sigma LKJ density degrees of freedom. When 0, the usual inverse Wishart prior is used for SIGMAs. When SLKJDF>0, then the LKJ density is used as the prior, with SLKJDF degrees of freedom. In addition, only diagonal elements of the Sigma prior are used. SLKJDF may be set >0 when using METHOD=BAYES as well, but SIGMAs (in cholesky format) will then be M-H sampled using the PSAMPLE_M1, PSAMPLE_M2, and PSAMPLE_M3 settings (choleskys of sigma elements are treated as extensions of the THETA parameters in M-H sampling methods). See record $SLKJDF below to specify LKJ correlation degrees of freedom for each sigma block.

7.14 SVARF=1 (default)

SVARF is the weight factor to STD prior to the log sqrt Sigma diagonal elements, the normal density of the log square root of Sigma centered about log square root of Sigma prior, and scaled with SVARF (see below). That is,

\[\log(\sqrt{\Sigma_i}) \sim \text{Normal}(\log(\sqrt{\Sigma_{\text{prior},i}}), 1/\text{SVARF}).\]

If SVARF<0, then a half-t-distribution of degrees of ABS(SVARF) is used as the prior to the sqrt of SIGMA diagonal elements. Use SVARF=-1 for the half-Cauchy distribution.

7.15 $OLKJDF (0 default, NM75)

The $OLKJDF is a separate record that allows the user to specify LKJ decorrelation degrees of freedom for each OMEGA block. For example:

$OLKJDF 4.5 3.5 -2.0

Where 4.5 degrees of freedom are specified for the first omega block, 3.5 for the second, and -2.0 specifies 2 degrees of freedom for the third omega block, but a user-defined definition of the standard deviations of the diagonals for the third omega block. For any blocks without degrees of freedom defined, it is 0. Use the OMEGA_STD_PRIORU.f90 file in the ..\source directory as a template to modify the OMEGA_STD_PRIORU to provide the desired probability density function. Set $SUBR OTHER=OMEGA_STD_PRIORU.f90 to use the user modified template (you may rename the file for organizational purposes). $EST OLKJDF over-rides $OLKJDF.

7.16 $OVARF 0 (default, NM75)

The OVARF is a separate record that allows the user to specify the weighting (inverse variance) to the standard deviations LKJ decorrelation degrees of freedom for each OMEGA block. For example:

$OVARF 2.0 5.0

Where 2.0 is specified for the first omega block, 5.0 for the second. If the corresponding $OLKJDF value is negative then this is argument STDSSP in user-defined OMEGA_STD_PRIORU.f90. $EST OVARF over-rides $OVARF.

7.17 $SLKJDF 0 (default)

The SLKJDF is a separate record that allows the user to specify LKJ decorrelation degrees of freedom for each SIGMA block. For example:

$SLKJDF 4.5 3.5 -2.0

Where 4.5 degrees of freedom are specified for the first sigma block, 3.5 for the second, and -2.0 specifies 2 degrees of freedom for the third sigma block, but a user-defined definition of the standard deviations of the diagonals for the third sigma block. For any blocks without degrees of freedom defined, it is 0. Use the SIGMA_STD_PRIORU.f90 file in the ..\source directory as a template to modify the SIGMA_STD_PRIORU to provide the desired probability density function. Set $SUBR OTHER=SIGMA_STD_PRIORU.f90 to use the user modified template. $EST SLKJDF over-rides $SLKJDF.

7.18 $SVARF 0 (default, NM75)

The SVARF is a separate record that allows the user to specify the weighting (inverse variance) to the standard deviations LKJ decorrelation degrees of freedom for each SIGMA block. For example:

$SVARF 2.0 5.0

Where 2.0 is specified for the first sigma block, 5.0 for the second. If the corresponding $SLKJDF value is negative then this is argument STDSSP in the user-defined SIGMA_STD_PRIORU.f90. $EST SVARF over-rides $SVARF.

7.19 $TTDF 0 (default, NM75)

The ttdf is a separate record that allows the user to specify the t-distribution degrees of freedom for each theta. So,

$TTDF (2.0)x2 (3.0)x2

sets t-distribution degrees of freedom for THETAs 1 and 2 to value of 2, and DF for THETAs 3 and 4 are set to 3.0. By default $TTDF value is 0, for any THETAs whose degrees of freedom is not specified. These will be used in estimation (integer or real) and simulation (truncated integer only) unless $EST TTDF or the $SIM TTDF is set, respectively, which set the t degrees of freedom to a single TTDF value for all THETAs.

7.20 NUTS_TRANSFORM=0 (default)

When NUTS_TRANSFORM=0, model parameters are transformed using the mass matrix, for population parameters. If NUTS_TRANSFORM=1, the momentum parameters are transformed using the mass matrix. It is best to set NUTS_TRANSFORM to 1 if NUTS_TEST is set to 1, and NUTS_TRANSFORM should be set to 0 when NUTS_TEST is set to 0.

7.21 NUTS_EPARAM=0 (default)

When NUTS_EPARAM=0, parameters are entered into the NUTS algorithm parameterized as THETAs and phis, Cholesky. When NUTS_EPARAM=1, parameters are entered into the NUTS algorithm parameterized as THETAs and etas. When NUTS_EPARAM=2, parameters are entered into the NUTS algorithm parameterized as THETAs and Choleksy of Omega inverse*eta.

7.22 NUTS_OPARAM=1 (default)

When NUTS_OPARAM=1, Omega elements are parameterized in a correlation cholesky format that constrains correlations to be between -1 and 1. When NUTS_OPARAM=0, then the full Omega elements are parameterized in cholesky format.

7.23 NUTS_SPARAM=1 (default)

When NUTS_SPARAM=1, Sigma elements are parameterized in a correlation cholesky format that constrains correlations to be between -1 and 1. When NUTS_SPARAM=0, then the full Sigma elements are parameterized in cholesky format.

7.24 NUTS_REG=0.0 (default)

By default, the mass matrix is made slightly diagonal dominant by adding a fraction of the diagonal element. If NUTS_REG>0.0, then the mass matrix is made slightly diagonal dominant by adding the value of NUTS_REG. When OLKJDF>0, then NUTS_REG=1.0 may provide a more efficient sampling process.

7.25 NUTS_STEPITER=1 (default)

An initial step size is calculated for the first NUTS_STEPITER iterations.

7.26 NUTS_STEPINTER=0 (default)

An initial step size is calculated every NUTS_STEPINTER iterations.

7.27 NUTS_TEST=0 (default)

The acceptance of a sample is tested using an algorithm as originally in algorithm 6 of Hoffman and Gelman (Hoffman and Gelman 2014) (default, NUTSTEST=0), or as performed in STAN (NUTSTEST=1)

7.28 NUTS_INIT=0.075 (default)

When using the STAN algorithm (MADAPT=-1) for mass matrix and step size development during the burn-in (warmup) stage, when NUTS_INIT<1 serves as the fraction of NBURN iterations for Stage I of the warmup period ((Stan Development Team 2015)). When NUTS_INITS>1, then the explicit number of iterations is interpreted. Similarly, when MADAPT>0, this period is also used to accumulate NUTS_INIT*NBURN iterations before using the mass matrix.

7.29 NUTS_BASE=0.025 (default)

When using the STAN algorithm (MADAPT=-1) for mass matrix and step size development during the burnin-in (warmup) stage, NUTS_BASE (if NUTS_BASE>=1) or NUTS_BASE*NBURN (if NUTSBASE<1) serves as the number of iterations for the first Stage II of the warmup period ((Stan Development Team 2015)), and doubles in iteration number with each subsequent segment of Stage II. The total number of stage II iterations is bounded by NBURN-NUTSINIT-NUTS_TERM. If NUTS_INIT=75, NUTS_BASE=25, NUTS_TERM=150, and NBURN=1000, then you have 5 segments of stage II, so that

NUTS_INIT+NUTS_TERM+NUTS_BASE+NUTS_BASE*2+NUTS_BASE*4=NBURN=
75 + 150 +25 +25*2+25*4+25*8+25*16=1000

When MADAPT>0, this period is also used to update the mass matrix every NUTSBASE iterations until MADAPT total iterations have been performed.

If NUTS_BASE<=-1.0, then NUTS_BASE will be set to the largest block section of the mass matrix plus 10. This assures that a large enough base set of samples are collected before the mass matrix is used.

If NUTS_BASE<-1, then in addition, the number of stage II iterations is ABS(NUTS_BASE). The actual NBURN will be based on the above equation, but not to exceed the user specified NBURN, which serves as the max NBURN. With NUTS_BASE<-1.0, set NBURN to a large number (4000 or so). The AUTO feature set NUTS_BASE to -3.

7.30 NUTS_TERM=0.05 (default)

When using the STAN algorithm (MADAPT=-1) for mass matrix and step size development during the burn-in (warmup) stage, NUTS_TERM serves as the number of iterations for Stage III of the warmup period ((Stan Development Team 2015)), to make final adjustments in step size.

7.31 MUFIRSTREC, OBJQUICK

For simple problems, the NO U-Turn process can have excessive overhead and run slowly. To have it run faster, you can do the following:

Set MUFIRSTREC=1 in $PRED or $PK. MUFIRSTREC=1 selects the covariate of the first record of the subject, rather than averaging among its records when using that covariate in a MU reference.

OBJQUICK=0. Standard NONMEM processing of the model occurs.

OBJQUICK=1. Certain tests and initializations are skipped.

OBJQUICK=2. A simplified modeling process occurs, but which cannot be used when $LEVEL or $MIX is used in the model. Also, parallelization is not performed.

Usage:

$PRED
include nonmem_reserved_general
MUFIRSTREC=1
OBJQUICK=1

The OBJQUICK and MUFIRSTREC can also speed up the other analysis methods, such as ordinary BAYES, FAST FOCE, ITS, and the EM methods.

7.32 Note on Combinations of Option Settings NUTS_MASS, NUTS_EPARAM, NUTS_OPARAM, NUTS_SPARAM

Certain combinations of option settings work well, others do not.

Default: NUTS_EPARAM=0 NUTS_MASS=B

These settings are the most efficient for many of the problems tested so far. They offer the greatest speed efficiency and sampling (Neff/Nsample) efficiency. On occasion, one or two THETAs will have low efficiencies relative to the rest. The AUTO=1 option allows an easy setup of this configuration (see section I.45 Some General Options and Notes Regarding EM and Monte Carlo Methods for more details on the AUTO feature, AUTO=0 (default) (NM73)). Example stanrb42.ctl uses the AUTO=1 feature.

NUTS_EPARAM=1 NUTS_MASS=D

When the problem is submitted to the NUTS algorithms with etas rather than phis (NUTS_REPARAM=1), the NUTS_MASS=B does not yield good efficiencies on THETAs. Therefore, NUTS_MASS=D needs to be used. However, this reduces speed efficiencies by about 5 fold, but evens out the theta efficiencies so the lowest Neff/Nsample efficiency is about 3x higher than the lowest sampling efficiency of the NUTS_EPARAM=0 setting. The AUTO=3 option allows an easy setup of this configuration.

NUTS_EPARAM=2 NUTS_MASS=BD

With NUTS_EPEARM=2, this is called the “Matt trick” in the Stan community, and offers very high sampling efficiencies of 3 to 5 fold than that of the default settings. However, the speed is about 4-8 times slower from the default settings when performing NUTS in NONMEM, so there may or may not be greater overall efficiency, in terms of number of independent samples per unit time. Furthermore, this method uses a conditional likelihood equation that differs from the standard that the population analysis community is used to: It drops the NIND*LOG(DET(OMEGA)) term (where NIND=number of subjects), resulting in a conditional likelihood that has very different distribution properties and does not fully represent the contribution of the individual etas to the likelihood. An example of using NUTS_EPARAM=2 is ..\examples\stanrb19. The AUTO=2 option allows an easy setup of this configuration, and example ..\examples\stanbrb39 uses this option.

The Stan community also supports the notion of using LKJ correlation priors for supporting OMEGAs, rather than inverse Wishart priors. This is reasonable when using LKJ decorrelation as an uninformative prior, and there is no previous knowledge of the scale of the variances. For example, one could use Identity matrix for the Omega priors (diagonal OMEGA values=1, and very low degrees of freedom, and the LKJ decorrelation prior will not introduce much bias, whereas the inverse Wishart prior would introduce considerable bias. However, the uninformative inverse Wishart prior is reasonable to use if the diagonal Omega values are in a reasonable range of where the data are. In PK/PD modeling, we have the benefit of obtaining reasonable variances from first performing a maximum likelihood analysis, (using FOCE, ITS, IMP, ITS, or SAEM), which often do not require any priors, and then supplying these results as priors for the NUTS analysis, as long as the degrees of freedom is set to <=D, the dimension of the block matrix. In such cases, the Inverse Wishart as an uninformative prior (DF<=D), but with the variances obtained from an earlier maximum likelihood analysis on the same data, is equivalent to LKJ correlation prior in terms of quality and lack of bias.

When using Omega information from a previous study to supply as an informative prior to a present study, the Inverse Wishart format of the prior information conjugates well with the information in the cross products of etas provided by the present study, and the informative prior information from the previous study offers a natural statistical support, as if the data of the previous study were added to the present study. Furthermore, the inverse Wishart prior supplies the information of the entire block, off-diagonals and diagonals, whereas the LKJ correlation prior method (OLKJDF>0) only uses the diagonal elements, and some general notion of correlation in the OLKJDF value. Such a natural interpretation for the inverse Wishart is evident in the mathematical structure of the total likelihood, when using NUTS_EPARAM=0 or NUTS_EPARAM=1, and Omega priors with inverse wishart distribution.

7.33 Prepare a Single Burn-In for Multiple Stationary Chains

Because the burn-in of NUTS is computationally expensive, you can use the MSF file to perform a one-time burn-in in one control stream file, and then have subsequent control stream files use the burn-in information (including the mass density information) and evaluate their stationary phases in parallel (at different starting seeds). For example, start with the following burn-in control stream, where NITER=0:

$PROB RUN# Example 1 (from samp5l)
$INPUT C SET ID JID TIME DV=CONC AMT=DOSE RATE EVID MDV CMT CLX
            V1X QX V2X SDIX SDSX
$DATA example1.csv IGNORE=C
$SUBROUTINES ADVAN3 TRANS4
$PK
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
CL=DEXP(MU_1+ETA(1))
V1=DEXP(MU_2+ETA(2))
Q=DEXP(MU_3+ETA(3))
V2=DEXP(MU_4+ETA(4))
S1=V1
$ERROR
Y = F + F*EPS(1)
$THETA (2.0)X4
$OMEGA BLOCK(4) VALUES(0.15,0.01)
$SIGMA (0.6 ) ;[P]
$PRIOR NWPRI
$THETAP (2.0 FIX)X4
$THETAPV BLOCK(4) FIX VALUES(10000,0.0)
$OMEGAP BLOCK(4) FIX VALUES(0.2,0.0)
$OMEGAPD (4 FIX)
$SIGMAP 0.06 FIX
$SIGMAPD (1 FIX)
$EST METHOD=ITS NITER=0 PRINT=10 file=nuts_resume.ext
$EST METHOD=NUTS INTERACTION AUTO=1 NITER=0 PRINT=10 file=nuts_resume.txt
         MSFO=nuts_resume.msf
$COV MATRIX=R PRINT=E UNCONDITIONAL

You can use the MSF file to resume an analysis, in this case with NITER non-zero (the stationary phase), and specify a seed:

$PROB RUN# Example 1 (from samp5l)
$INPUT C SET ID JID TIME DV=CONC AMT=DOSE RATE EVID MDV CMT CLX
            V1X QX V2X SDIX SDSX
$DATA example1.csv IGNORE=C
$SUBROUTINES ADVAN3 TRANS4
$PK
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
CL=DEXP(MU_1+ETA(1))
V1=DEXP(MU_2+ETA(2))
Q=DEXP(MU_3+ETA(3))
V2=DEXP(MU_4+ETA(4))
S1=V1
$ERROR
Y = F + F*EPS(1)
$PRIOR NWPRI NTHETA=4, NETA=4, NTHP=4, NETP=4, NEPS=1, NEPP=1
$MSFI nuts_resume.msf
$EST METHOD=NUTS AUTO=1 NITER=1000 PRINT=10 file=nuts_resume2.txt
         SEED=2234556
$COV MATRIX=R PRINT=E UNCONDITIONAL

You can run additional control streams as the one above, in parallel, each with a different starting seed. These will serve as separate chains.

8 A Note on Setting up Prior Information

Prior information is important for MCMC Bayesian analysis, but not necessary for maximization methods. Of greatest importance are priors to the OMEGAs. As a general rule, if your data set consists of fewer subjects than 100 times the dimension of the Omega matrix to be estimated, then you should have at least uninformative OMEGA prior information. Priors to THETAS are assumed multivariate normal, and priors to OMEGAS and SIGMAS are assumed inverse Wishart distributed. Alternatively, a residual variance in the form of its square root, may be modeled via THETA (a sigma-like Theta parameters is set up in example 2). For a thorough reference to the options in the $PRIOR record, see ..\guides\VIII.pdf. The following describes the setup for most Bayesian analysis purposes.

To set up the $PRIOR NWPRI statement, keep in mind the following:

• NTHETA=number of THETAs to be estimated
• NETA=number of Etas (OMEGAs) to be estimated (and is to be described by an NETAxNETA OMEGA matrix)
• NEPS=number of epsilons (SIGMAs) to be estimated (and is to be described by an NEPSxNEPS SIGMA matrix)
• NTHP=number of THETAs which have a prior
• NETP=number of OMEGAs with prior
• NEPP=number of SIGMAs with prior (NM73). Before NM73, the NEPP option was ignored, as supplying priors for SIGMAs was not activated.

For example:

$PRIOR NWPRI NTHETA=4, NETA=4, NEPS=1 NTHP=4, NETP=4, NEPP=1

Then the $THETA records list the parameters, in order, the following:

• NTHETA of initial THETAs
• NTHP of Priors to THETAS
• Degrees of freedom to each OMEGA block Prior
• Degrees of freedom to each SIGMA block Prior

The $OMEGA records list the variances, in order, the following:

• NETAxNETA of initial OMEGAS
• NTHPxNTHP of variances of Priors to THETAS
• NETPxNETP of priors to OMEGAS, matching the block pattern of the initial OMEGAS

The $SIGMA records list the variances, in order, the following:

• NEPSxNEPS of initial SIGMAS
• NEPPxNEPP of priors to SIGMAS, matching the block pattern of the initial SIGMAS (NM73).

So we may have the following example control stream file portion:

$THETA 2.0 2.0 4.0 4.0 ; Initial THETAs
$OMEGA BLOCK(4) ; Initial Parameters for OMEGA
0.4
0.01 0.4
0.01 0.01 0.4
0.01 0.01 0.01 0.4
$SIGMA 0.1

$PRIOR NWPRI NTHETA=4, NETA=4, NEPS=1, NTHP=4, NETP=4, NEPP=1

; Prior information of THETAS (NTHP=4 of them)
$THETA (2.0 FIX) (2.0 FIX) (2.0 FIX) (2.0 FIX)

; Variance to prior information of THETAS (NTHPxNTHP=4x4 of them).
; Because variances are very large, this means that the prior
; information to the THETAS is highly uninformative. Note that the
; order of $THETA values among the THETA records, and the order
; of $OMEGA values among the OMEGA records, is very important,
; But $THETAs and $OMEGAs can be interspersed.
$OMEGA BLOCK(4)
10000 FIX
0.00 10000
0.00 0.00 10000
0.00 0.00 0.0 10000

; Prior to OMEGA (NETPxNETP=4x4 if them)
$OMEGA BLOCK(4)
0.2 FIX
0.0 0.2
0.0 0.0 0.2
0.0 0.0 0.0 0.2

; Set degrees of freedom of OMEGA Prior (one value per OMEGA block)
; Uninformative Omega prior is designated by having a DF that is equal to
; the dimension size of the Omega block.
$THETA (4 FIX)

; Prior to SIGMA (NEPPxNEPP=1x1 if them)
$SIGMA 0.05 FIX

; Set degrees of freedom of SIGMA Prior (one value per SIGMA block)
; Uninformative SIGMA prior is designated by having a DF that is equal to
; the dimension size of the Sigma block.
$THETA (1 FIX)

By default, the number of prior experiments is 1. However, perhaps you have more than one previous study, and you wish to average their contribution, forming a composite average set of prior parameters to influence the present analysis. In this case, add NEXP=n to the $NWPRI record above, where n is the number of experiments. Then, add the prior information of each additional study with additional $THETA, $OMEGA, and $SIGMA statements. The order is then:

$THETA records list the parameters, in order, the following:

• NTHETA of initial THETAs
• Exp 1:
  • NTHP of Priors to THETAS
  • Degrees of freedom to each OMEGA block Prior
  • Degrees of freedom to each SIGMA block Prior
• Exp 2:
  • NTHP of Priors to THETAS
  • Degrees of freedom to each OMEGA block Prior
  • Degrees of freedom to each SIGMA block Prior
• …

The $OMEGA records list the variances, in order, the following:

• NETAxNETA of initial OMEGAS
• Exp 1:
  • NTHPxNTHP of variances of Priors to THETAS
  • NETPxNETP of priors to OMEGAS, matching the block pattern of the initial OMEGAS
• Exp 2:
  • NTHPxNTHP of variances of Priors to THETAS
  • NETPxNETP of priors to OMEGAS, matching the block pattern of the initial OMEGAS
• …

The $SIGMA records list the variances, in order, the following:

• NEPSxNEPS of initial SIGMAS
• Exp 1:
  • NEPPxNEPP of priors to SIGMAS, matching the block pattern of the initial SIGMAS
• Exp 2:
  • NEPPxNEPP of priors to SIGMAS, matching the block pattern of the initial SIGMAS

Additional examples of setting up prior information for various problems are shown in the example problems listed at the end of this document.

As of NM73, you can use more informative names as follows:

• $THETAP for theta priors
• $THETAPV for variance to theta priors
• $OMEGAP for omega priors
• $OMEGAPD for degrees of freedom (or dispersion factor) for omega priors
• $SIGMAP for SIGMA priors
• $SIGMAPD for degrees of freedom (or dispersion factor) for SIGMA priors

This allows you to intersperse these records at will in the control stream files, but it also gives NM-TRAN an alternative source for values to NTHETA, NETA, NTHT, NETP, NEPS, and NEPP that is typically given in the $PRIOR NWPRIOR record. However, if these values are also listed in $PRIOR NWPRI, then these values are chosen over what is surmised from the informatively labeled theta/omega/sigma records. Thus, the above control stream file could be structured as follows, with the various records in any order, and a shortened $PRIOR record (in the following example uninformative priors are used):

$PRIOR NWPRI

; Prior information of THETAS (NTHP=4 of them)
$THETAP (2.0 FIX) (2.0 FIX) (2.0 FIX) (2.0 FIX)

$THETA 2.0 2.0 4.0 4.0 ; Initial THETAs
$OMEGA BLOCK(4) ; Inital Parameters for OMEGA
0.4
0.01 0.4
0.01 0.01 0.4
0.01 0.01 0.01 0.4

; Set degrees of freedom of SIGMA Prior (one value per SIGMA block)
$SIGMAPD (1 FIX)

;intial parameters to sigma
$SIGMA 0.1

; Set degrees of freedom of OMEGA Prior (one value per OMEGA block)
$OMEGAPD (4 FIX)

; Prior to OMEGA (NETPxNETP=4x4 if them)
$OMEGAP BLOCK(4)
0.2 FIX
0.0 0.2
0.0 0.0 0.2
0.0 0.0 0.0 0.2

; Variance to prior information of THETAS (NTHPxNTHP=4x4 of them).
$THETAPV BLOCK(4)
10000 FIX
0.00 10000
0.00 0.00 10000
0.00 0.00 0.0 10000

; Prior to SIGMA (NEPPxNEPP=1x1 if them)
$SIGMAP 0.05 FIX

Informative prior information may come from a previous study. Typically, they are used as follows:

The theta priors for the present analysis are obtained from the estimates of THETAs from the previous study. For example, in the report file of the previous study:

FINAL PARAMETER ESTIMATE

THETA - VECTOR OF FIXED EFFECTS PARAMETERS

TH 1        TH 2          TH 3  TH 4

1.64E+00 1.57E+00 7.58E-01 2.35E+00

would be placed in the present study control stream file as:

$THETAP (1.64 FIXED) (1.57 FIXED) (0.758 FIXED)
                  (2.35 FIXED)

The variance-covariance to theta priors of the present analysis are obtained from the variance-covariance submatrix pertaining to the theta estimates from the previous study. For example, the information in the report file of the previous study:

COVARIANCE MATRIX OF ESTIMATE

            TH 1          TH 2   TH 3           TH 4
                                           2.98E-03
   TH 1

+        2.33E-03

   TH 2

+        4.76E-04 2.86E-03

   TH 3

+        7.87E-04 1.27E-04 5.35E-03

   TH 4

+        7.80E-05 2.36E-04 1.76E-03

would be placed in the control stream file of the present study as:

$THETAPV BLOCK(4)                5.35E-03  2.98E-03
              2.33E-03 FIXED     1.76E-03
              4.76E-04 2.86E-03
              7.87E-04 1.27E-04
              7.80E-05 2.36E-04

The omega priors of the present analysis are obtained from the estimates of OMEGAs from the previous study. For example, from the report file of the previous study:

OMEGA - COV MATRIX FOR RANDOM EFFECTS - ETAS

            ETA1          ETA2   ETA3      ETA4

   ETA1

+        1.75E-01

   ETA2

+        8.33E-03 1.51E-01

   ETA3

+        2.98E-02 1.74E-02 2.41E-01

   ETA4

+        -8.05E-03 1.84E-02 5.14E-02 1.62E-01

you transpose as follows to the control stream of the present study:

$OMEGAP BLOCK(4)          2.41E-01  1.62E-01
       1.75E-01 FIXED     5.14E-02
       8.33E-03 1.51E-01
       2.98E-02 1.74E-02
     -8.05E-03 1.84E-02

Similarly for Sigma priors, the results of the previous study:

SIGMA - COV MATRIX FOR RANDOM EFFECTS - EPSILONS ***

            EPS1

   EPS1

+        5.28E-02

Is transposed to the present study control stream as:

$SIGMAP (5.28E-02 FIXED)

The degrees of freedom to the omega priors of the present analysis are at most the total number of subjects in the previous study. Dr. Mats Karlsson has proposed the following formula for selecting degrees of freedom (justification is described in appendix G of the NONMEM7 Technical Guide), which is appropriate when the previous analysis was an nlme (maximum likelihood) form:

\[\text{DF} = 2 \times \left[\frac{\text{Omega estimate of previous analysis}}{\text{SE of omega of previous analysis}}\right]^2\]

Or

\[\text{DF} = 2 \times \left[\frac{\text{Omega estimate of previous analysis}}{\text{SE of omega of previous analysis}}\right]^2 + 1\]

to adjust for degrees of freedom loss in the estimate of Omega of the previous study.

For an OMEGA block, use the smallest DF calculated among the OMEGA diagonal estimates in that block.

A similar formula would apply for SIGMA priors, with the proviso that the DF be no larger than the total number of data points that apply for that sigma in the previous study (for example, if there are two SIGMAs, one for PK data, and another for PD data, then the sigma for PK data gets no more than total number of PK data points in the previous study).

If the previous analysis was that of a Bayesian analysis, the final mean OMEGAs (or SIGMAs) are expanded relative to the modal values from a maximum likelihood analysis, and the ratio of mean OMEGAs to variance OMEGAs from a Bayesian analysis is smaller than the ratio of model OMEGAs to its Fisher information assessed variance that one obtains from a maximum likelihood/nlme analysis. As a result, the recommended transformation when incorporating results from a previous Bayesian analysis are as follows (where p is the dimension of the omega block, Q is the variance matrix from the previous analysis):

OMEGAPD:

\[\text{DF} = 2 \times \left[\frac{\text{Omega estimate of previous analysis}}{\text{SE of omega of previous analysis}}\right]^2 + p + 3\]

OMEGAP:

\[R = \left(\frac{\text{DF} - p - 1}{\text{DF}}\right) \times Q\]

where R is the matrix that you would place as OMEGAP. The justification of the above is based on the mean and variance formulas for an inverse Wishart (see appendix G of the NONMEM7 Technical Guide).

The SCALE option available for $OMEGA(P) and $SIGMA(P) records allows the user to conveniently input a re-scaling of the input matrix without having to hand calculate every element (see SCALE option for $OMEGA and $SIGMA (NM760) in section I.8 Expansions on Abbreviated and Verbatim Code and Other Items (NM72,NM73,NM74,NM75).

For convenient transfer of the entire information from a previous analysis to prior information for a subsequent analysis, see section I.100 priorget: Transfer Results of an Analysis to NMTRAN Prior Information (NM75).

As yet another alternative, the mean of the inverse OMEGAS or SIGMAS from MCMC Bayesian samples and their variances are such that they are similar to those obtained from an nlme analysis. Therefore, the inverse of the mean inverse OMEGA can be used as the input $OMEGAP without re-scaling, and the effective DF for $OMEGAPD would be similar to that of an nlme analysis:

\[ \text{DF} = 2 \times \left[ \frac{ \text{average Omega inverse of previous analysis} }{ \text{SE of omega inverse of previous analysis} } \right]^2 \]

The priorget utility will calculate the inverse of the average omega inverse from the non-negative iterations in the .ext file of an MCMC Bayesian analysis, and provide the appropriate prior information from these.

As of nm74, the degrees of freedom to the inverse wishart algorithms used for OMEGAS and SIGMAS may be any real number greater than 0. Thus, the inverse wishart matrix distribution can substitute for inverse gamma matrix distribution as follows. The parameter beta is the rate parameter (with inverse units of the variable), and alpha is the shape parameter, to a gamma distribution. This gamma distribution to the inverse residual variance can be expressed with an equivalent Wishart distribution to the inverse residual variance. Set 2*alpha for the $SIGMAPD, and beta/alpha for $SIGMA. The gamma distribution to the inverse residual variance is equivalent to the inverse gamma distribution of the residual variance.

9 Monte Carlo Direct Sampling (NM72)

On rare occasions, direct Monte Carlo sampling may desired. This method is the purest method for performing expectation maximization, in that it creates completely independent samples (unlike MCMC), and there is no chance of causing bias if the sampling density is not similar enough to the conditional density (unlike IMP). However, it is very inefficient, requiring ISAMPLE values of 10000 to 300000 to properly estimate the problem. The method can be implemented by issuing a command such as

$EST METHOD=DIRECT INTERACTION ISAMPLE=10000 NITER=50

On occasion it can have some use in jump starting an importance sampling method, especially if the first iteration of importance sampling fails because it relies on MAP estimation, and the problem is too unstable for it. Thus, one could perform the following, where just a few iterations of direct sampling begin the estimation process:

$EST METHOD=DIRECT INTERACTION ISAMPLE=10000 NITER=3
$EST METHOD=IMP INTERACTION ISAMPLE=1000 NITER=50 MAPITER=0

Notice that since MAPITER=0, the first iteration of IMP method relies on starting parameters for its sampling density that came from the DIRECT sampling method.

10 Some General Options and Notes Regarding EM and Monte Carlo Methods

10.1 AUTO=0 (default) (NM73)

If option AUTO=1 is selected, then several options will be set by NONMEM that will allow best settings to be determined. The user may still override those options set by AUTO, by specifying them on the same $EST record. For example,

$EST METHOD=ITS AUTO=1 PRINT=10
$EST METHOD=SAEM AUTO=1 PRINT=50
$EST METHOD=IMP PRINT=1 EONLY=1 NITER=5 ISAMPLE=1000
$EST METHOD=BAYES AUTO=1 NITER=1000 FILE=auto.txt PRINT=100

The settings of AUTO=1 for each method are as follows:

METHOD=DIRECT INTERACTION ISAMPLE=1000 CTYPE=3 NITER=500 STDOBJ=10
                                    ISAMPEND=10000 NOPRIOR=1 CITER=10 CINTERVAL=0 CALPHA=0.05
                                    EONLY=0

METHOD=BAYES INTERACTION CTYPE=3 NITER=10000 NBURN=4000
                                    NOPRIOR=0 CITER=10 CINTERVAL=0 CALPHA=0.05
                                    IACCEPT=0.4 ISCALE_MIN=1.0E-06 ISCALE_MAX=1.0E+06
                                    PACCEPT=0.5 PSCALE_MIN=0.01 PSCALE_MAX=1000
                                    PSAMPLE_M1=1 PSAMPLE_M2=-1 PSAMPLE_M3=1 OSAMPLE_M1=-1
                                    OSAMPLE_M2=-1 OACCEPT=0.5 ISAMPLE_M1=2 ISAMPLE_M1A=0
                                    ISAMPLE_M2=2 ISAMPLE_M3=3 ISAMPLE_M1B=2 MCETA=0

METHOD=SAEM INTERACTION CTYPE=3 NITER=1000 NBURN=4000
                                    ISAMPEND=10 NOPRIOR=1 CITER=10 CINTERVAL=0 CALPHA=0.05
                                    IACCEPT=0.4 ISCALE_MIN=1.0E-06 ISCALE_MAX=1.0E+06
                                    ISAMPLE_M1=2 ISAMPLE_M1A=0 ISAMPLE_M1B=2 ISAMPLE_M2=2
                                    ISAMPLE_M3=2 CONSTRAIN=1 EONLY=0 ISAMPLE=2 MCETA=0

METHOD=ITS INTERACTION CTYPE=3 NITER=500
                                    NOPRIOR=1 CITER=10 CINTERVAL=1 CALPHA=0.05

METHOD=IMP INTERACTION CTYPE=3 NITER=500 ISAMPLE=300 STDOBJ=10
                                    ISAMPEND=10000 NOPRIOR=1 CITER=10 CINTERVAL=1 CALPHA=0.05
                                    IACCEPT=0.0 ISCALE_MIN=0.1 ISCALE_MAX=10 DF=0 MCETA=3
                                    EONLY=0 MAPITER=1 MAPINTER=-1

METHOD=IMPMAP INTERACTION CTYPE=3 NITER=500 ISAMPLE=300 STDOBJ=10
                                    ISAMPEND=10000 NOPRIOR=1 CITER=10 CINTERVAL=1 CALPHA=0.05
                                    IACCEPT=0.0 ISCALE_MIN=0.1 ISCALE_MAX=10 DF=0 MCETA=3
                                    EONLY=0

As of nm74, for IMP estimation a second auto value, AUTO=2 is available. Same settings as AUTO=1, with additional:

GRDQ=-1.0 DERCONT=1 RANMETHOD=3S2P

As of nm75, for SAEM estimation, a second auto value, AUTO=2 is available. Same settings as AUTO=1, with additional:

MAPITERS=1
MCETA=100

to turn MAP assessment assist on the first iteration.

As of nm75, for BAYES estimation, a second auto value, AUTO=2 is available. Same settings as AUTO=1, with additional:

MAPITERS=1
MCETA=100

to turn on MAP assessment assist on the first iteration.

As of nm74, an AUTO=1 feature is available for NUTS algorithm

METHOD=NUTS INTERACTION CTYPE=0 NITER=2000 NBURN=10000
              NOPRIOR=0 NUTS_STEPITER=1 NUTS_STEPINTER=0 NUTS_TEST=0
              NUTS_INIT=75 NUTS_BASE=-3 NUTS_TERM=50 NUTS_GAMMA=0.05
              NUTS_DELTA=0.8 KAPPA=1.0 IKAPPA=1.0
              NUTS_REG=0.0 MADAPT=-1
              NUTS_EPARAM=0 NUTS_OPARAM=1 NUTS_SPARAM=1 NUTS_MASS=B
              NUTS_TRANSFORM=0 NUTS_MAXDEPTH=10

A second auto value, AUTO=2, may be used with NUTS estimation to setup the alternative sampling strategy, “Matt trick” (options that differ from AUTO=1 are shown in bold):

METHOD=NUTS INTERACTION CTYPE=0 NITER=2000 NBURN=10000
              NOPRIOR=0 NUTS_STEPITER=1 NUTS_STEPINTER=0 NUTS_TEST=0
              NUTS_INIT=75 NUTS_BASE=-3 NUTS_TERM=50 NUTS_GAMMA=0.05
              NUTS_DELTA=0.8 KAPPA=1.0 IKAPPA=1.0
              NUTS_REG=0.0 MADAPT=-1
              NUTS_EPARAM=2 NUTS_OPARAM=1 NUTS_SPARAM=1 NUTS_MASS=BD
              NUTS_TRANSFORM=0 NUTS_MAXDEPTH=10

A third auto value, AUTO=3, may be user with NUTS estimation to setup the alternative sampling strategy of eta sampling (options that differ from AUTO=1 are shown in bold):

METHOD=NUTS INTERACTION CTYPE=0 NITER=2000 NBURN=10000
              NOPRIOR=0 NUTS_STEPITER=1 NUTS_STEPINTER=0 NUTS_TEST=0
              NUTS_INIT=75 NUTS_BASE=-3 NUTS_TERM=50 NUTS_GAMMA=0.05
              NUTS_DELTA=0.8 KAPPA=1.0 IKAPPA=1.0
              NUTS_REG=0.0 MADAPT=-1
              NUTS_EPARAM=1 NUTS_OPARAM=1 NUTS_SPARAM=1 NUTS_MASS=D
              NUTS_TRANSFORM=0 NUTS_MAXDEPTH=10

The AUTO option is ignored by the FO/FOCE/Laplace methods. The AUTO setting itself transfers to the next $EST within the same $PROB, just like any other option settings explicitly set by the user in the control stream file, so AUTO remains on or off until then next AUTO option specified. For example, in the following example:

$EST METHOD=ITS AUTO=1 PRINT=10
$EST METHOD=SAEM AUTO=1 PRINT=50
$EST METHOD=IMP PRINT=1 EONLY=1 NITER=5 ISAMPLE=1000
$EST METHOD=BAYES AUTO=1 FILE=auto.txt PRINT=100 NITER=1000

the IMP statement also has AUTO=1. However, for the following example:

$EST METHOD=ITS AUTO=1 PRINT=10
$EST METHOD=SAEM AUTO=1 PRINT=50
$EST METHOD=IMP PRINT=1 EONLY=1 NITER=5 ISAMPLE=1000 AUTO=0
$EST METHOD=BAYES AUTO=1 FILE=auto.txt PRINT=100 NITER=1000

the AUTO setting is turned off for IMP, and turned back on for BAYES. Any option settings implicitly set by the AUTO feature does not transfer to the next $EST statement. Also, when using AUTO=1, the transfer of any options settings explicitly set by the user from previous $EST statements may or may not occur for those options set by the AUTO option, depending on the situation.

10.2 S-ADAPT / NONMEM parameter mapping

The mapping of parameters between S-ADAPT and NONMEM is as follows:

S-ADAPT	NONMEM
Pmethod=4	IMPMAP
Pmethod=8	IMP
Pmethod=1	ITS
Pmethod=6	DIRECT
Npopiter	NITER
Npopc	ISAMPLE
Npop	MCETA
optmethod	OPTMAP
covest	ETADER
Gefficiency	IACCEPT
Gamma_min	ISCALE_MIN
Gamma_max	ISACLE_MAX
DFRAN	DF
Popconv_test	CTYPE
Popconv_rows	CITER
Popconv_alpha	CALPHA
Ndelpar	MAPINTER
Poperr_type=3	$COV MATRIX=S
Poperr_type=8	$COV MATRIX=R
Poperr_type=9	$COV
POPFINAL subroutine	CONSTRAINT subroutine may be user modified to provide any constraining pattern on any population parameters
RANMETHOD	RANMETHOD
SEED	SEED

Hoffman, Matthew D., and Andrew Gelman. 2014. “The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.” Journal of Machine Learning Research 15: 1593–623.

Hong, Hee Sun, and Fred J. Hickernell. 2003. “Algorithm 823: Implementing Scrambled Digital Sequences.” ACM Transactions on Mathematical Software 29 (2): 95–109. https://doi.org/10.1145/779359.779360.

Lavielle, Marc. 2007. MONOLIX Users Manual, Version 2.1. Laboratoire de Mathématiques, Université Paris-Sud.

Lavielle, Marc, and Benjamin Ribba. 2016. “Enhanced Method for Diagnosing Pharmacometric Models: Random Sampling from Conditional Distributions.” Pharmaceutical Research 33 (12): 2979–88. https://doi.org/10.1007/s11095-016-2020-3.

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.

Stan Development Team. 2015. Stan Modeling Language Users Guide and Reference Manual, Version 2.6.3. https://mc-stan.org/docs/.