69 Bootstrap, Selecting a Random Method, and Other Options for Simulation (NM73, NM74)
Introduction to NONMEM 7.6.0
1 BOOTSTRAP (NM73)
$SIML BOOTSTRAP=-1 SUBP=100
$EST METHOD=1 INTERACTIONThe above example requests a bootstrap rearrangement (with replacement) of an existing data set, followed by analysis of that data set. The BOOTSTRAP number refers to how many subjects are to be randomly selected from the data set. Setting -1 or to a value larger than the number of subjects in the data set means to randomly select as many subjects as are in the data set. For example, if 400 subjects are in the simulation template data set, then 400 subjects are randomly selected (with replacement, so some are selected more than once, others not at all). In this case, NONMEM’s simulator does not perform the usual activity of randomly creating DV values for a new data set, but rather selects a random set of subjects of an existing data set (which must already have legitimate DV values), uniformly selected (using source 1) with replacement. This results in some subjects not being selected at all, and some subjects selected more than once.
2 NOREPLACE (NM73)
$SIML BOOTSTRAP=50 SUBP=100 NOREPLACE
$EST METHOD=1 INTERACTIONIn the above example, 50 unique subjects are to be randomly selected from the simulation template data set. The NOREPLACE feature is reasonable if there are many more than 50 subjects to choose from template set (for example, 1000 subjects in the template, and for each sub-problem, 50 of them are randomly chosen without replacement, that is, without repeating a subject).
3 STRAT (NM73)
$SIML BOOTSTRAP=50 SUBP=100 NOREPLACE STRAT=CATA single stratification data item may be entered. In the above example, the data item CAT serves as the stratification. This splits the data set into distinct sub-sets, guaranteeing a specific number of subjects will be selected from each category. For example, if in the base data set CAT has values of 1 or 2, with 33 subjects in group 1 and 67 subjects in group 2 out of 100 total subjects, then exactly 33% of of subjects from group 1 will be randomly selected out of 50 total (16), and exactly 67% of subjects will be randomly selected from group 2 (34). This has value when desiring that a bootstrap analysis maintain the same proportion of subjects belonging to certain categories, such as gender, or age bracket. To stratify by both age bracket and gender, create a stratification data item that would be, for example, valued 1 for subjects who are male under 30, 2 for subjects that are female under 30, 3 for subjects who are male over 30, 4 for subjects who are female over 30. Any discrete numerical values will do, as long as the stratifier is not a continuous variable, and the subjects need not be sorted according to the stratification data item.
4 STRATF (NM73)
$SIML BOOTSTRAP=50 SUBP=100 NOREPLACE STRAT=CAT STRATF=FCATThe option STRATF points to a data item that contains the fraction that should represent a category in the bootstrapped data set. Without STRATF, the number of subjects to be taken from a given category is proportional to the number of subjects in the base data set. If you want the category to be represented at a different proportion, then specify a STRATF data item, in this example, FCAT. Suppose FCAT=0.5 for CAT=1 and 0.5 for CAT=2 as well. Even though only 33% of subjects in the base data set belong to category 1, exactly 50% of subjects from group 1 will be randomly selected out of 50 total (25), and exactly 50% of subjects will be randomly selected from group 2 (25) in the formation of each bootstrap data set. This allows you to alter the proportions in each category from what is in the original data set.
5 BOOTDATA=0 (Default) (NM75)
By default (BOOTDATA=0), when data are selected based on $SIML BOOSTRAP, the randomly selected subjects are analyzed during the subsequent estimation method. If BOOTDATA=1, then the subjects not selected are analyzed. Usually, you may wish to simply perform an evaluation, not an estimation, on these non-selected data (MAXEVAL=0, or NITER=0), at the population parameter values based on estimation of the selected subjects. For example:
$SIML BOOTSTRAP=-1 SUBP=100
$EST METHOD=1 MAXEVAL=9999 BOOTDATA=0 ; Estimated using randomly selected subjects
$EST METHOD=1 MAXEVAL=0 BOOTDATA=1 ; Evaluate non-selected subjects at estimated parameters from selected subjects.Or for EM methods:
$SIML BOOTSTRAP=-1 SUBP=100
$EST METHOD=ITS NITER=50 BOOTDATA=0 ; Estimated using randomly selected subjects
$EST METHOD=ITS NITER=0 BOOTDATA=1 ; Evaluate at estimation of parameters6 Bootstrapping Single Subject Data (NM74)
The BOOTSTRAP option in $SIML is most suitable for multi-subject data, in which there is an ID data column identifying the subjects. However, because NONMEM will parse data records belonging to a single subject in single subject data as if they belonged to separate subjects, the BOOTSTRAP option can be used to create a series of data sets with random samples of all records in the data set, and this may have some utility in obtaining summary estimates statistics for an individual curve fit. This is shown in the following control3boot.ctl example:
$PROBLEM THEOPHYLLINE SINGLE SUBJECT DATA
$INPUT DOSE=AMT TIME CP=DV CAT
$DATA DATA3B
$SUBROUTINES ADVAN2
$PK
CALLFL=1
KA=THETA(1)
K=THETA(2)
SC=THETA(3)
$ERROR
IPRED=F
W=1.0
; first observation after dose is part of "first subject". So, put in dummy record, CAT=3,
; and give it a residual variance that is very large, so it does not influence the fit.
IF(CAT==3.0) W=1.0E+10
Y=F+W*ERR(1)
$THETA (0,1.7) (0,.102) (0,29)
$OMEGA 0.2
$SIML (567666 NORMAL) (33012 UNIFORM) BOOTSTRAP=-1 STRAT=CAT SUBP=100
$ESTIMATION MAXEVAL=9999 PRINT=2
$COVR
$TABLE TIME CAT AMT CP IPRED W NOAPPEND NOPRINT file=control3boot.tabThe data file is as follows:
320 .0 .1
. .1 3.0 3
. .27 1.71 2
. .52 7.91 2
. 1. 8.31 2
. 1.92 8.33 2
. 3.5 6.85 2
. 5.02 6.08 2
. 7.03 5.4 2
. 9. 4.55 2
. 12. 3.01 2
. 24.3 .90 2Note that to modify the original control3 problem, the CAT data item was added. Each dose record should have its own category number, so that it will be presented in every random dataset, rather than sometimes there and sometimes not, and sometimes duplicated. In this case, the single dose record is given category 1. The observed data are given category 2, from which a random set will be selected for each data set created. Furthermore, the dose record and the first data record after it is considered “one individual”, and therefore, this first data record (in this example the one for time=0.1) will always travel with the dose record. Since this means that the time=0.1 data record is not randomly selected, but always present in every random set, it was made to be a “dummy” record, one that was not originally in the control3 problem. To ensure that this dummy record does not contribute to the estimate, it is flagged to have a residual deviation of W=1.0E+10, so that it is given very little weight, and does not impact the estimate. One can readily give it a CAT=3, to distinguish it from the other data records during estimation, although its category number is not recognized as a separate category during the random data set creation, since the data record is associated with the dose record with category 1.
7 PARAFILE (NM74)
As of NONMEM 7.4, there is the option to perform simulations in a parallel computing setting, if parallel computing was requested, by the –parafile option on the command line, for example (see section I.75 Parallel Computing (NM72)). By default, simulations are not parallel computed because they are performed rapidly. There may be occasions, however, when the model is sufficiently complicated, and/or there are many subjects in the template data set, that simulation could benefit from parallel computing. If you wish to turn parallel computing on for the simulation step, then set PARAFILE to ON:
$SIM PARAFILE=ON …Alternatively, set –simparon on the nmfe75 command line to turn on parallelization during simulation step. The –simparon switch takes precedence.
When modeling with super-ID nested ETA levels ($LEVEL record is present), parallelization will not occur, since these etas are shared across individuals, and there is no guarantee that all subjects sharing the same etas will be simulated by the same process.
8 RANMETHOD=[n|S|m|P] (NM73)
As of NM73, the RANMETHOD option is available for the $SIM record, to use alternative random numbers generators (default is NONMEM’s traditional one, number 4):
$SIML RANMETHOD=[n|S|m|P]Where n is the random number generator type, S is Sobol sequence, and m is the Sobol scrambler. See the description of RANMETHOD under I.38 Monte Carlo Importance Sampling EM.
As of NM74, RANMETHOD will also act on the P modifier, which will retain separate random number sequences for each subject, so that the random variable patterns are retained regardless of whether the simulation is done in single computing or parallel computing mode. So, when parallel computing, you may select
$SIM PARAFILE=ON RANMETHOD=PBefore NM74, a P descriptor could be given without an error message, but it was ignored because the parallel random number generation for each individual was not yet developed for $SIM.
NONMEM’s default random number generator for the $SIM step is 4 (in contrast, default random number generator for $EST and $TABLE is 3). Number 4 is NONMEM’s classic random number generator. Whatever random number source is selected, it affects all seed1 sources, and all source seed2 if not also selecting P. The P descriptor sets initial seed2 to 0 when creating starting seeds for each individual.
The Sobol method is used only to generate normally distributed random vectors of etas and epsilons, when the S descriptor is selected, and SEED1 of source 1 is used to set the seed. Among the Sobol sequence methods, the S2 method appears to provide the least biased random samples, that is nearly uniform distribution, with good mixing in multi-dimensional spaces.
9 REWIND (NM74)
The REWIND feature in $SIM allows the original data set to be used for all sub-problems. By default, if any data item is changed by a sub-problem, those data items remain changed for the start of the next sub-problem. If you want that each sub-problem start with using the values from the original data set, use the REWIND feature of $SIM. So, any changes to the data set made during simulation (when ICALL=4) of a sub-problem are used for that sub-problem only, and are not preserved for the next sub-problem. Keep in mind that any transgeneration you may have performed on the data set when using an $INFN when ICALL=1 will be considered original data set. For example:
$INFN
IF (ICALL==1) THEN
DOWHILE(DATA)
..modifying statements here
ENDDO
ENDIF10 NOSUPRESET (NM74)
By default, (SUPRESET), with subsequent iterations of a super-problem, the simulation random number generator is re-initialized back to the seed that is listed in the $SIM record of the control stream file. It may be desirable that each iteration serves as a new random instance, so use NOSUPRESET.
11 Simulation Error Forgiveness (NM72)
As of NM72, if a simulation error occurs (such as “KA nonpositive”), or a user-implemented EXIT 1 is invoked during simulation, then another eta and eps sample set will be generated and tested. If ten such errors occur in the same subject, then it is supposed that the cause of the simulation error is not due to an occasional bad random sample, but is caused by a systematic error in the control stream file, and the simulation is terminated.
12 Extensions to Simulation Error Forgiveness (NM75)
You may wish to test a predictive value result, and if unacceptable, request a new set of etas and epsilons. For example,
$ERROR
…
IF(ICALL==4.and.IPRED<0.01 .and. TIME>20.0) EXIT 1 2300If the first value following the EXIT statement is 1 (out of a possible range of 0-2), and the second value is between 1000 and 9999 (out of a possible range of 0-9999), then upon the condition being true, NONMEM will try another random eta and eps sample. Be careful that the condition does not occur too often (causing wasteful computation) or that it always occurs due to logic error in the control stream, or NONMEM can be caught in an infinite loop.
Note that if you wish to simply create a truncated normal distribution that only requires testing the eta value directly, this can be done with:
$PK
…
DO WHILE (ETA(1).GT.5)
CALL SIMETA(ETA)
ENDDOas suggested in Guide VIII, in ABBREVIATED CODE section.
If it is desired that the simulation be immediately terminated, then generate an EXIT 2 code:
IF(ICALL==4.and.IPRED<0.1 .and. TIME>20.0) EXIT 2This will also over-ride the simulation error forgiveness.
13 Providing Separate Random Number Sequences for Etas and Epsilons (NM75)
By default, etas and epsilons are generated from the same random source 1, with epsilons generated after etas are generated, for each subject, in sequence. You may wish to retain the sequence of randomly generated etas (for variables) using one source, and retain a separate sequence of randomly generated epsilons (for observed data) using another source. To do this, specify the random source index with SOURCE_EPS, and provide starting seed values for both sources, as follows
$PROB RUN# Example 1 (from samp5l)
$INPUT C SET ID JID TIME DV=CONC AMT=DOSE RATE EVID MDV CMT
$DATA seedeps.csv IGNORE=C
…
$SIMULATION ONLYSIMULATION SUBPROBLEMS=5
SOURCE_EPS=2 ; specify that epsilons should be generated from source 2
(567811 NORMAL) ; starting seed for Source 1 for Etas
(2567811 NORMAL); stating seed for Source 2 for epsilons
(2933012 UNIFORM) ; starting seed for source 3
; for something else (for calling RANDOM(R,3)).This is useful if you generate data from two separate data file templates which differ only in the number of data points among subjects, but you have the same number of subjects, and you wish to retain the same etas between the two data sets generated across subproblems.
14 Simulating THETAs with t-Distribution (NM75)
By setting TTDF to a non-zero integer value, THETAs can be simulated with TTDF degrees of freedom t-distribution, when setting also TRUE=PRIOR. For example:
$SIMULATION (567811 NORMAL) (2933012 UNIFORM) SUBPROBLEMS=10 ONLYSIM
TTDF=4 TRUE=PRIORMake sure to also have the priors appropriately set up, such as:
$PRIOR NWPRI PLEV=0.999
$THETAP (1.68 FIXED) (1.58 FIXED) (8.12E-01 FIXED) ( 2.37 FIXED)
$THETAPV BLOCK(4)
0.04 FIX
0.00 0.4
0.00 0.00 0.4
0.00 0.00 0.0 0.4
$OMEGAP BLOCK(4)
1.65276E-01 FIXED
4.62368E-03 1.33770E-01
6.35274E-03 1.69155E-02 2.13881E-01
-1.53098E-02 1.25855E-02 5.32524E-02 1.62771E-01
$OMEGAPD (50 FIX)
$SIGMAP BLOCK(1) 0.0544 FIXED
$SIGMAPD (300 FIXED)See see $TTDF 0 (default, NM75) to specify degrees of freedom for each theta. Specifying $SIML TTDF over-rides the $TTDF settings.
15 CLOCKSEED=0 (default) (nm75)
As of nm75, starting seed will be 10000*(seconds after midnight)+SEED (SEED may be set to 0 for this purpose), if CLOCKSEED=1, for each source seed (SEED1) or source seed pair (SEED1, SEED2). 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. If SEED2<=0, then the clock time will not be added to SEED2.