94  Clinical Trial Design Evaluation and Optimization (NM75)

Introduction to NONMEM 7.6.0

Author

Robert J. Bauer

Published

June 1, 2025

The design process can help you in designing or evaluating a clinical trial. It may be desired to evaluate specified time points, or find the optimal time points, dose levels, and number of time points appropriate for a particular sub-design, etc. The design algorithms have been modeled after POPED by Hooker et al., and PFIM by Mentre et al. (see references (Foracchia et al. 2004; Ogungbenro et al. 2005; Nyberg, Ueckert, et al. 2012; Nyberg et al. 2014; Lestini et al. 2015; Dumont et al. 2018)). The following features are available. Please also see the tutorial on clinical trial design (Bauer et al. 2021), which provide an extended tutorial on the most essential options and examples shown here, as well as additional examples.

$DESIGN

followed by a set of options.

1 APPROX=FO (default)

The nlme approximation method is specified here. First order (no interaction) is the default, and the appropriate type of covariance matrix is evaluated or used in optimization of the design. Other options are:

  • FOI: First order interaction with interaction.
  • FOCE: first order conditional estimation.
  • FOCEI: first order conditional estimation with interaction.
  • LAPLACE: Laplace conditional estimation.
  • LAPLACEI: Laplace conditional estimation with interaction.

2 OFVTYPE=1 (default)

The objective function types are comparable to those of PopED:

  • 0,1,3,4,5: design type: D-optimality, -log(det(FIM)), where FIM=Fisher Information Matrix (inverse of variance-covariance).
  • 2: design type: A-optimality, -1/tr(1/FIM).
  • 6: design type: DS-optimality, -log(det(FIM))+log(det(FIMuninteresting)).

To identify parameters as uninteresting, place UNINT at the parameter in the same manner you would place FIX. For example:

$THETA 0.1 (0.5 UNINT) 0.4

$OMEGA BLOCK(2)
0.1 UNINT ; Sets entire block to uninteresting
0.01 0.1
  • 7: design type: R-optimality (relative standard error), -1/tr(sqrt(1/FIM)/Parameter)).
  • 8: optimal design type: Individual Bayesian FIM, -log(det(Bayes FIM)), as described in the PFIM 4.0 manual (Lestini et al. 2015).

In this case, the ETC() are optimized, as listed in the .phi table. The ETC() are the conditional variance-covariances of the individual’s parameters. If just one subject (or subject type) is in the data file, then the criterion is that one subject’s Bayesian FIM. If more than one subject (or subject type), then the Bayesian FIM of all the subjects, averaged together, is the criterion for design.

The change in variances with each printed iteration are listed in the additional raw file (.ext extension by default). The shrinkage information is reported in the main report file, and the individual empirical Bayes variances (ETC()) are reported in the .phi file. For shrinkage information, only the following have meaning:

  • EBVSHRINKSD(%)
  • EBVSHRINKVR(%)
  • EPSSHRINKSD(%)
  • EPSSHRINKVR(%)

When OFVTYPE=8, the change of average individual conditional variances (average empirical Bayes conditional variance) with each iteration are shown in the root.bfm file (where root is name of control stream file), and the final one on the -1000000000 line, during an optimization. The RAW file (by default named root.ext) shows only the starting, and final values of the standard errors of population parameters, as extra information.

  • 9: Same as option 2, and using the UNINT filter.
  • 10: Same as option 7, and using the UNINT filter.

3 OFVDIMSCALED=0 (default)

If OFVDIMSCALED=1, then the objective function is scaled to be “per parameter”, the dimension of the FIM, as follows:

  • 0,1,3,4,5: design type: D-optimality, -log(det(FIM))/D, where D is the dimension of the FIM.
  • 2: design type: A-optimality, -1/(tr(1/FIM)/D).
  • 6: design type: DS-optimality, (-log(det(FIM))+log(det(FIMuninteresting)))/(D-U) where U=number of uninteresting parameters.
  • 7: design type: R-optimality (relative standard error), -1/(tr(sqrt(1/FIM)/Parameter))/D).
  • 8: -log(det(Bayes FIM))/D.
  • 9: design type: A-optimality, -1/(tr(1/FIM)/(D-U)).
  • 10: -1/(tr(sqrt(1/FIM)/Parameter))/(D-U)).

4 GROUPSIZE=1 (default)

The GROUPSIZE is comparable to that of POPED, in which the FIM is multiplied by this number to provide the subject number size of the dataset template. For a template of one subject, GROUPSIZE would then offer the variance-covariance expected from GROUPSIZE number of subjects.

5 FIMTYPE or FIMDIAG=0

FIMTYPE or FIMDIAG may be set to the following, and corresponds to fimcalc.type in POPED:

0 (default): Full information matrix, using finite difference assessment for Theta, Omega, and Sigma variances and covariances.

No assumptions are made about whether the approximated variance (C=hSigmah+gOMEGAg, as described in Guide I, section E.2) used in FO is dependent on any THETAs via h and g.

1: Create a block diagonal information matrix of the estimates that assumes there is neglible correlation between THETAs and (OMEGAs,SIGMAs). Also, independence of variance C with respect to THETAs is assumed (that is, h and g’s dependence on theta is ignored). This is similar to the fim.calc.type=1 option in POPED or the diagonal option in PFIM, and performs finite difference at the y=E(f) position on the thetaxtheta portion of the information matrix (for THETAs that are mu referenced, an analytical evaluation is made, and increases speed considerably), analytical evaluation of (SIGMAS OMEGAS)x(SIGMAS OMEGAS) portion of information matrix, and THETAs x (SIGMAS OMEGAS) cross variances are set to 0.

2: Create a block information matrix, by performing finite difference on the thetaxtheta portion of the information matrix, independence of variance C with respect to THETAs is not assumed (that is, h and g’s dependence on theta is not ignored), performs analytical evaluation of (SIGMAS OMEGAS)x(SIGMAS OMEGAS) portion of information matrix, and THETAs x (SIGMAS OMEGAS) are cross variances set to 0.

3: Create a full information matrix by performing finite difference on the Theta x Theta portion of the information matrix, and on the THETAs x (SIGMAS OMEGAS) portion of the information matrix, and performs analytical evaluation of (SIGMAS OMEGAS)x(SIGMAS OMEGAS) portion of information matrix.

For the film.calc.type=4 equivalent, see VARCROSS below.

With FIMTYPE=1, it is worth while to mu reference all non-fixed THETAs. This will allow analytical evaluation of the thetaXtheta portion of the information matrix, and can provide considerable speed improvement. Even derivatives of mu-referenced THETAs with associated OMEGA() fixed to 0 will be analytically evaluated.

6 VARCROSS=0 (default)

By default (VARCROSS=0), NONMEM interprets residual variance such as the following:

$PK
…
IPRED=F
Y=IPRED+EPS(1)+IPRED*EPS(2)
…
$SIGMA 0.5 0.15

as residual variance var=SIGMA(1,1)+IPREDIPREDSIGMA(2,2)

To model variance as residual standard deviation squared:

var=(sqrt(SIGMA(1,1))+IPRED*sqrt(SIGMA(2,2))**2= SIGMA(1,1)+2sqrt(SIGMA(1,1)SIGMA(2,2))IPRED+IPREDIPRED*SIGMA(2,2)

as in the manner of the PFIM 4.0 software, one would normally model the residual variance using THETAs, such as:

$PK
…
IPRED=F
W=SQRT(THETA(4))+SQRT(THETA(5))*IPRED
Y=IPRED+W*EPS(1)
…
$SIGMA (1.0 FIXED)

However, when FIMTYPE=1 is used, the theta(4) and theta(5) are not identified as Sigma type parameters, and the appropriate contribution by theta(4) and theta(5) is not estimated. To provide an interpretation of residual variance incorporating the the cross-term

2sqrt(SIGMA(1,1)SIGMA(2,2))*IPRED

model using the SIGMAS as shown in the first example, and set VARCROSS to 1:

$PK
…
IPRED=F
Y=IPRED+EPS(1)+IPRED*EPS(2)
…
$SIGMA 0.5 0.15
…
$DESIGN … VARCROSS=1 FIMTYPE=1

FIMTYPE=1 VARCROSS=1 is equivalent to fim.calc.type=4 in POPED, and diagonal option in PFIM. Note that the variances are still modeled, rather than the standard deviation form. See .._pfim.ctl for an example of its use.

7 EOPTD=1

For each iteration, this creates a random sample of THETAs, OMEGAs, or SIGMAs, using the prior information. $PRIOR NWPRI prior information is required, and PLEV=0.999 must be specified. Best used with STGR. See example optdesign16.

8 SEED=223345

Set the starting seed for any random samples to be created, whether for EOPTD=1, or for FOCE type FIM in creating random etas (see below).

9 CLOCKSEED=0 (default)

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.

10 MODE=0 (default), 1, or 2

Used for specifying data and prediction value type when specifying APPROX=FOCEI.

MODE=0 means FOCEI with data at F(ETAsim), and predicted function evaluated at f(ETAsim), is to be used. This method works well.

MODE=2 means linearized FOCEI, with data at F(ETAsim)-G*ETAsim, and predicted function at F(ETAsim). Works well.

MODE=1 means FOCEI with data at F(ETAsim), and predicted function evaluated at the mode, f(ETAhat), is to be used (see below). The results are not satisfactory.

11 DATASIM=0 (default)

Normally, y-expectation evaluation of the FIM is performed. To actually simulate data, set DATASIM=1, and along with APPROX=FOCEI, this will produce simulated etas as well.

So,

$DESIGN APPROX=FOCEI MODE=1 NELDER FIMDIAG=0 DATASIM=1 GROUPSIZE=32 OFVTYPE=0

will produce the most empirical, “clinical trial simulation” (CTS) style covariances, complete with simulated etas and eps, and standard FIM is assessed. If FIMDIAG>0, then a y-expectation covariance will be evaluated, but mode will be evaluated with the simulated data.

12 Additional Control Options for $DESIGN

The following options may be set within the $DESIGN record, and they operate exactly as their equivalents in the $EST record. In fact, these options specified in $DESIGN are transferred to the ESTM record in the FCON file, and NONMEM treats these in the same manner as if they came from the $EST record, but with the interpretation that is appropriate for the optimal design or evaluation setting.

  • SIGL
  • SIGLO
  • ABORT/NOABORT/NOHABORT
  • MAXEVAL (=0 indicates design evaluation, which is default, whereas MAXEVAL>0 indicates design optimization)
  • PRINT (control iteration printing during optimal design)
  • NUMERICAL/NONUMERICAL
  • -2LL/LIKELIHOOD/LLIKELIHOOD
  • SLOW/NOSLOW/FAST
  • POSTHOC
  • NOPRIOR
  • FORMAT
  • FILE
  • FNLETA

Thus, the following are equivalent:

$DESIGN FIMDIAG=1 OFVTYPE=6 APPROX=FO
              MAXEVAL=9999 NOHABORT PRINT=20 SIGL=10 POSTHOC

$DESIGN FIMDIAG=1 OFVTYPE=6 APPROX=FO
$EST MAXEVAL=9999 NOHABORT PRINT=20 SIGL=10 POSTHOC

In addition, $DESIGN sets up the covariance step as $COV MATRIX=R UNCONDITIONAL without the user requiring this record entered in the control stream. If you wish to specify additional control for the covariance, you can add these in a $COV record, such as:

$DESIGN FIMDIAG=1 OFVTYPE=6 APPROX=FO
              MAXEVAL=9999 NOHABORT PRINT=20 SIGL=10 POSTHOC

$COV PRINT=E

13 An Example of a Design Evaluation

The options described above can be used to set up the design criteria, and can be used for finding the best design parameters (optimal design), or for evaluating design parameters proposed by the user (design evaluation). We will consider a simple example of design evaluation, before we consider additional options needed for optimal design. Using the example of .._pfim.ctl, note the following record:

$DESIGN GROUPSIZE=200 FIMTYPE=1 VARCROSS=1

The csv file contains only one subject (that is one subject type), but the subject number size may be specified with GROUPSIZE. Here, you request the variance-covariance one should expect from 200 subjects, all of the one particular subject type (equivalent time values, doses, covariates) defined in the csv file. You can have more than one subject in the csv, each with their own design (different doses, sample times, etc.). Also, by use of STRAT and STRATF (see below), one can desire a specific weighting of each subject type to contribute to the FIM.

FIMTYPE=1 means to evaluate the variance-covariance with the block-diagonal modality. VARCROSS=1 means to treat the residual variance model in the manner of PFIM 4.0, as described in the manual.

The intended result for this problem is essentially the FIM information and its derived components. In warfarin_pfim.coi, is the inverse covariance matrix (the FIM itself), for a sample size of 200 subjects. In the .res file, the -log(det(PFIM)) is given. Also, standard errors are given in report form here. The .ext file has essentially the same information is the .res file, but is in rows-column format. Much of these outputs (.res, .ext., .coi) are already a part of standard NONMEM outputs, they are now just recruited for optimal design and design evaluation.

Note the usual variance-covariance intended for an actual data set is not created, but rather the asymptotic, y-expectation evaluation version, since the $DESIGN record was specified. Much of the control stream setup is as is if one were performing an FO estimation analysis on a real data set, while the addition of the $DESIGN record then re-purposes the analysis for design evaluation (or optimization if MAXEVAL>0).

View the standard error and variance-covariance outputs in the NONMEM report file, .ext file, .cov, .cor, and .coi files. For optimal design, you can view the final optimal design in the final data set, which can be outputted with the $TABLE step (see Examples for Optimal Design, below).

14 Options for Setting up Types of Optimal Design

The additional options for $DESIGN listed below are for optimizing parts of the design components. For example, the DESEL, DESELSTRAT, DESELMIN, DESELMAX can be specified for all the various design elements that you want optimized. It might be TIME for time samples, AMT for dose, or some type of covariate specific for the problem. Certainly, any combination of covariates can be requested to be optimized.

14.1 NELDER

Use Nelder method to search for optimal continuous parameters.

14.2 FEDOROV

Use to find ideal set of discrete time points from a larger set of possible time points.

14.3 RS

Random Search method to find optimal continuous parameters.

14.4 STGR

Stochastic Gradient method to find optimal continuous parameters.

14.5 DISCRETE

Find optimal number of time points for each sub-design (subject template), and use NELDER method to find optimal continuous parameters.

14.6 DISCRETE_RS

Find optimal number of time points for each sub-design (subject template), and use RS method to find optimal continuous parameters.

14.7 DISCRETE_SG

Find optimal number of time points for each sub-design (subject template), and use STGR method to find optimal continuous parameters.

Next, data items must be specified to be optimized. This is done using the following options:

14.8 DESEL=data item

The data item (column name) that contains the design element (DESEL) values that are to be modified and optimized. For example, TIME column would indicate that you want times to be estimated.

14.9 DESELSTRAT=data item indicating stratification

The DESELSTRAT data item should contain integer indices to indicate what values are to be shared, and estimated together. If a record contains a value of 0 in the DESELSTRAT column, then this record is not included in the estimation process, and its value (say its time value in DESEL=TIME column) will not be changed. If the record contains a value >0 in DESELSTRAT, let us suppose a 1, then all records with the value of 1 in DESELSTRAT will share the same time value (or whatever DESEL selected), estimated together. Those records with value 2 will be another set of records which will share a time value, etc. Thus, within a subject, there may be a PK record and a PD record which should share the same time value. Also, a group of subjects may share the same time values. Within a subject, times will be automatically sorted as they are changed, so that NONMEM’s time ordering policy is not violated.

14.10 DESELMIN=data item containing the minimal value

14.11 DESELMAX=data item containing the maximal value

You must impose boundaries on the values that are being optimized. That is done with these two data items. Only those records with a stratification value >0 in DESELSTRAT column will require a min and max value, and only those records that define that stratification value for the first time.

DESEL, DESELMIN, DESELMAX, and DESELSTRAT may be repeated for as many design elements are to be optimized. For example for times and amounts:

DESEL=TIME DESELMMIN=TMIN DESELMAX=TMAX DESELSTRAT=TSTRAT
DESEL=AMT DESELMIN=AMTMIN DESELMAX=AMTMAX DESELSTRAT=AMTSTRAT

14.12 NMIN=data item containing minimal number of time points to the subject

14.13 NMAX=data item containing maximal number of time points to the subject

If NMIN<0 or NMAX<0, then most previous non-negative value is used. The NMIN and NMAX column are only used for the DISCRETE* analyses, to bound the number of time points that may be permitted for a given subject. With DISCRETE*, the total N of time points among all subjects is determined by the total number of time points whose MDV=0 in the starting data set.

14.14 STRAT=data item containing grouping or stratification number pertaining to that subject

14.15 STRATF=data item containing starting fraction representation for the stratification value in data item STRAT

If STRAT and STRATF are specified, and there is at least one STRAT value >0, then the STRATF values are optimized, and represent the weight to the contribution of that subject to the Information matrix. For STRAT<=0, then their STRATF values are not optimized, and remain fixed at their initial values, but are still used as weights to the information matrix. It is up to the user to ensure that the sum of STRATF values among unique STRAT values sum to 1. If value of STRATF<0.0, then that subject is not included in the assessment.

15 Examples of Optimal Design

Here are portions of some examples of scripts located in examples. Note that results are reported in the outputted tables. Study the data sets as well, as these establish part of the setup.

The following use FO by default.

Optdesign.ctl; FEDOROV search for ideal time points (fixed number), from a possible set of discrete time points given in the data set. The data set contains many time points, only 5 of which have MDV=0 (see optdesignw.csv). The FEDOROV algorithms tries different combinations of 5 time points among those given in the data set, and finds the best ones. Final (optimal) data set is listed in optdesign.tab.

$DESIGN FEDOROV FIMTYPE=1 MAXEVAL=9999 SIGL=12
$TABLE ID TIME EVID MDV DV NOPRINT NOAPPEND EXCLUDE_BY XCLD

            FILE=optdesign.tab

Optdesign2.ctl; Nelder search for ideal time points (fixed number). The data set has initial times, and these times are varied continuously, until an optimal set are found. The algorithm assures that the new times are sorted appropriately so as not to violate NONMEM rules. Notice in this example that the Omega values are UNINTed, and only the THETAs and the residual proportional error are of interest.

$THETA                                 2.37436E+00
 1.68338E+00 1.58812E+00 8.12710E-01

;$OMEGA BLOCK(4) VALUES(0.0225,0.001)
$OMEGA (0.0225 UNINT)X4
$SIGMA
( 0.0225) ( 0.0001 UNINT)

$DESIGN NELDER FIMDIAG=1 OFVTYPE=6
              DESEL=TIME DESELSTRAT=TSTRAT DESELMIN=TMIN DESELMAX=TMAX
              APPROX=FO MAXEVAL=9999 NOHABORT PRINT=20 SIGL=10 POSTHOC

$TABLE ID TIME EVID MDV DV IPRED NOPRINT NOAPPEND FILE=optdesign2.tab
            FORMAT=S1PE23.16

Optdesign5.ctl; Search of best times, using stochastic gradient method of search.

$DESIGN FIMDIAG=0 STGR
              DESEL=TIME DESELSTRAT=TSTRAT DESELMIN=TMIN DESELMAX=TMAX APPROX=FO
              MAXEVAL=1000 SIGL=10 SIGLO=12 nohabort PRINT=50

$TABLE ID TIME EVID MDV DV NOPRINT NOAPPEND FILE=optdesign5.tab
FORMAT=S1PE23.16

Optdesign6.ctl; Search of best times, using Random Search method of search.

$DESIGN RS DESEL=TIME DESELSTRAT=TSTRAT DESELMIN=TMIN DESELMAX=TMAX
                    MAXEVAL=400 SIGL=12 nohabort PRINT=100

$TABLE ID TIME EVID MDV DV NOPRINT NOAPPEND FILE=optdesign6.tab
FORMAT=S1PE23.16

Optdesign11.ctl; Nelder search of best times, as well as optimizing for number of time points in each subject. Note that the MDV are set to 0 or 1, to add or remove data points to a subject. To select for STGR, use DISCRETE_SG (optdesign12.ctl), and for RS, use DISCRETE_RS (optdesign13.ctl).

$DESIGN DISCRETE FIMTYPE=1
              NMIN=NMIN NMAX=NMAX DESEL=TIME DESELSTRAT=TSTRAT DESELMIN=TMIN

DESELMAX=TMAX
              MAXEVAL=400 SIGL=10 nohabort PRINT=100

$TABLE ID TIME EVID MDV DV NOPRINT NOAPPEND FILE=optdesign11.tab
FORMAT=S1PE23.16

Optdesign14.ctl; Nelder search of best times, as well as optimizing the weight contributions to the FIM for each subject.

$DESIGN NELDER FIMTYPE=1 STRAT=STRAT STRATF=STRATF DESEL=TIME
DESELSTRAT=TSTRAT DESELMIN=TMIN DESELMAX=TMAX

                    MAXEVAL=400 SIGL=12 nohabort PRINT=10
$COV PRINT=E
$TABLE ID STRAT STRATF TIME EVID MDV DV NOPRINT NOAPPEND FILE=optdesign14.tab
FORMAT=S1PE23.16

Optdesign15.ctl; Nelder search of best times, DISCRETE optimization of number of data points for a subject, as well as optimizing the weight contributions to the FIM for each subject.

$DESIGN DISCRETE FIMTYPE=1 NMIN=NMIN NMAX=NMAX STRAT=STRAT
              STRATF=STRATF DESEL=TIME DESELSTRAT=TSTRAT
              DESELMIN=TMIN DESELMAX=TMAX
              MAXEVAL=400 SIGL=10 nohabort PRINT=10

$TABLE ID STRAT STRATF TIME EVID MDV DV NOPRINT NOAPPEND FILE=optdesign15.tab
FORMAT=S1PE23.16

Optdesign16.ctl; STGR search of best times, and random samples of THETA are generated at each iteration, for an “expected value” of FIM, since EOPTD=1. Note the need for $PRIOR information.

$PRIOR NWPRI PLEV=0.9999
$THETAP (1.68338E+00 FIXED) (1.58812E+00 FIXED) (8.12710E-01 FIXED)
(2.37436E+00 FIXED)
$THETAPV BLOCK(4) FIXED VALUES(0.03,0.001)

$DESIGN STGR EOPTD=1 FIMTYPE=1 SEED=1222345 DESEL=TIME DESELSTRAT=TSTRAT
DESELMIN=TMIN DESELMAX=TMAX

                    MAXEVAL=9999 SIGL=12 nohabort PRINT=10 NOPRIOR=1
$TABLE ID CLN V1N QN V2N TIME EVID MDV DV NOPRINT NOAPPEND
FILE=optdesign16.tab FORMAT=S1PE23.16

Optdesign17.ctl; The following evaluates FO with interaction for the FIM.

$DESIGN NELDER FIMTYPE=1
              DESEL=TIME DESELSTRAT=TSTRAT DESELMIN=TMIN DESELMAX=TMAX
              APPROX=FOI MAXEVAL=1000 SIGL=10 nohabort PRINT=100

$TABLE ID TIME EVID MDV DV NOPRINT NOAPPEND FILE=optdesign17.tab
FORMAT=S1PE23.16

Optdesign21.ctl; The following evaluates FOCE with interaction, non-MODE, for the FIM. FOCE is best executed with a number of replicate subjects, perhaps 100, each with their own random set of etas. You may use $SIM record to introduce seed values. Because of the new REPL option in $DATA, you can have a template data file that contains just one example subject for each design type desired, and ask for these to be replicated REPL times. For example, in optdesign21.csv, there is only one subject type, but one can want 100 replicates of these:

$DATA optdesign21.csv IGNORE=C REPL=100
…
$DESIGN NELDER MODE=2 FIMTYPE=1 GROUPSIZE=0.01 DESEL=TIME DESELSTRAT=TSTRAT
DESELMIN=TMIN DESELMAX=TMAX APPROX=FOCEI

                    MAXEVAL=400 SIGL=10 nohabort FNLETA=0
$SIM (1222345 NORMAL) (2933012 UNIFORM)
$TABLE ID ETA1 ETA2 ETA3 ETA4 TIME EVID MDV DV NOPRINT NOAPPEND
FILE=optdesign21.tab FORMAT=S1PE23.16

Also, notice that GROUPSIZE=0.01. In this way, 100 subjects are in the data set, but you may want the FIM to represent the FIM per subject, not the FIM of 100 subjects. Thus GROUPSIZE is set to 1/100.

Optdesign22.ctl; The following evaluates FOCE with interaction (modal) for the FIM. In this case, a SEED for source 1 may be entered at the $DESIGN record. If SEED is given at the $DESIGN record, as well as on a $SIM record, the SEED given at $DESIGN SEED option will have precedence for seed source 1.

$DATA optdesign21.csv IGNORE=C REPL=100
…
$DESIGN NELDER MODE=1 FIMTYPE=1 GROUPSIZE=0.01 SEED=1222345 APPROX=FOCEI

                    DESEL=TIME DESELSTRAT=TSTRAT DESELMIN=TMIN DESELMAX=TMAX
                    MAXEVAL=400 SIGL=10 nohabort FNLETA=0 PRINT=10
$TABLE ID ETA1 ETA2 ETA3 ETA4 TIME EVID MDV DV NOPRINT NOAPPEND
FILE=optdesign22.tab FORMAT=S1PE23.16

Optdesign25.ctl; The following evaluates FO for the FIM, for best times, and best doses.

$DESIGN NELDER FIMTYPE=1 DESEL=TIME DESELSTRAT=TSTRAT DESELMIN=TMIN
DESELMAX=TMAX

                    DESEL=AMT DESELSTRAT=AMTSTRAT DESELMIN=AMTMIN DESELMAX=AMTMAX
                    MAXEVAL=400 SIGL=10 nohabort
$TABLE ID TIME DOSE EVID MDV DV NOPRINT NOAPPEND FILE=optdesign25.tab
FORMAT=S1PE23.16

See example optex6d17.ctl for an A-optimality Nelder Mead search of best time points for a TMDD problem.

See example optex6d17_8.ctl for an example of Bayes FIM (OFVTYPE=8). The table optex6d17_8.tab reports the optimal data set, and optex6d17_8.bfm records the intermediate results of conditional variances during optimization. Final conditional variances are reported in optex6d17_8.phi.

You can also run several problems in one control stream, to use one optimizer followed by another. To do this across two problems note the following setup in warfarin_optimize.ctl:

$PROB 1
…
$DESIGN RS FIMDIAG=1 GROUPSIZE=32 OFVTYPE=0 DESEL=TIME DESELSTRAT=TSTRAT

                    DESELMIN=TMIN DESELMAX=TMAX
                    MAXEVAL=4000 SIGL=10 nohabort PRINT=100

$TABLE ID TIME EVID MDV DV IPRED NOPRINT NOAPPEND FILE=warfarin_optimize.tab
FORMAT=S1PE23.16
$TABLE ID TIME AMT RATE EVID MDV DV TSTRAT TMIN TMAX NOPRINT NOHEADER
NOAPPEND FILE=warfarin_1.CSV FORMAT=S1PE23.16
…
$PROB 2
$INPUT ID TIME AMT RATE EVID MDV DV TSTRAT TMIN TMAX
$DATA warfarin_1.csv NOOPEN
(1X,10(1X1PE23.16))

…
$DESIGN NELDER FIMDIAG=1 GROUPSIZE=32 OFVTYPE=0

              DESEL=TIME DESELSTRAT=TSTRAT DESELMIN=TMIN DESELMAX=TMAX
              MAXEVAL=4000 SIGL=10 nohabort PRINT=100
$TABLE ID TIME EVID MDV DV IPRED NOPRINT NOAPPEND FILE=warfarin_optimize2.tab
FORMAT=S1PE23.16
$TABLE ID TIME AMT RATE EVID MDV DV TSTRAT TMIN TMAX NOPRINT NOAPPEND
FILE=warfarin_nelder.CSV FORMAT=S1PE23.16

Notice the way in which $DATA of the second problem is used to pick up the file created in the first problem, to resume where the first problem left off regarding the modified time points.

You can also have one problem stack two or more $DESIGN records, as shown in warfarin_optimize3.ctl:

$DESIGN RS FIMDIAG=1 GROUPSIZE=32 OFVTYPE=0
                    APPROX=FOI SIGL=10 NOHABORT PRINT=100 MAXEVAL=400
                    DESEL=TIME DESELSTRAT=TSTRAT DESELMIN=TMIN DESELMAX=TMAX

$DESIGN STGR DESEL=TIME DESELSTRAT=TSTRAT DESELMIN=TMIN DESELMAX=TMAX
$DESIGN NELDER DESEL=TIME DESELSTRAT=TSTRAT DESELMIN=TMIN DESELMAX=TMAX

The only disadvantage to this method is that while the .ext and .phi files will record the intermittent change and results of each optimal design step, the main report file (.res or .lst) will only report the results of the final analysis, and tables will be outputted only after the last analysis, so you won’t be able to observe intermediate observation times.

The $DESIGN record will use the values of the options of the previous $DESIGN record for most options, except those of DESEL, DESELSTRAT, DESELMIN, DESELMAX. These options need to be specified for each $OPT record.

You can use the $SIML TRUE=PRIOR with prior information to randomly generate THETAs (and OMEGAs and SIGMAs, if desired), and optimal time points for each of these theta sets will be obtained (priortrue.ctl):

…
$PRIOR NWPRI PLEV=0.99
$THETAP (0.15 FIXED) (0.80 FIXED) (1.0 FIXED)
$THETAPV BLOCK(3) ; 30% CV
0.0025 FIX
0.0 0.07
0.0 0.0 0.11

$SIM (4442223) SUBPROB=100 TRUE=PRIOR
$DESIGN NELDER FIMDIAG=1 GROUPSIZE=32 OFVTYPE=0

              DESEL=TIME DESELSTRAT=TSTRAT DESELMIN=TMIN DESELMAX=TMAX
              MAXEVAL=4000 SIGL=10 nohabort PRINT=100 NOPRIOR=1
$TABLE ID TVCL TVV TVKA TIME EVID MDV DV IPRED NOPRINT NOAPPEND
FILE=priortrue.tab

In this example, 100 sub-problems are generated, each with their own theta sets. The variance on the THETAs is about 30%CV, and by observing the optimal times listed in priortrue.tab, you will get a sense of what kind of variability of optimal time points will be obtained.

As of nm751, you can add a prior from a previous analysis, or design, to a problem. Problem warfarin_desprior_opt4.ctl is a typical design evaluation, generating a .coi file. Then, warfarin_desprior_opt5.ctl considers another design, and adds the prior information in warfarin_desprior_opt4.coi, to the FIM:

…
$RCOVI FILE=warfarin_desprior_opt4.coi

$DESIGN GROUPSIZE=50 FIMTYPE=1 VARCROSS=1 APPROX=FO
              SIGL=12 MAXEVAL=9999 PRINT=100
              DESEL=TIME DESELSTRAT=TSTRAT DESELMIN=TMIN DESELMAX=TMAX

$TABLE ID TIME TSTRAT TMIN TMAX EVID MDV IPRED NOAPPEND NOPRINT
FILE=warfarin_desprior_opt5.tab
…

Such an example is given in reference (Dumont et al. 2018).

The same problem is done without $RCOVI in warfarin_desprior_opt.ctl. Note that the first subject represents the design of warfarin_desprior_opt5, and the other 4 subjects represents the elementary designs in warfarin_desprior_opt4, with their TSTRAT 0 so that they do not change with the optimization, but remain as unchanged, prior information.

You can also add prior information via $PRIOR NWPRI. Only the dispersion information is used, as centrality is assumed at the set population parameters. For example, $THETAPV, $OMEGAPD, and $SIGMAPD information is used, along with $SIGMA and $OMEGA, to construct the prior information. Example warfarin_pfim2p.ctl shows this:

$THETA
0.15 ;[CL]
8.0 ;[V]
1.0 ;[KA]

$OMEGA (0.07) (0.02) (0.6)
$SIGMA 0.01

$PRIOR NWPRI
$THETAP
0.15 FIXED ;[CL]
8.0 FIXED ;[V]
1.0 FIXED ;[KA]

$THETAPV BLOCK(3)
0.00005 FIXED
0.0 0.06
0.0 0.0 0.015
$OMEGAP (0.07 FIXED) (0.02 FIXED) (0.6 FIXED)
$OMEGAPD (32 FIXED) (40 FIXED) (56 FIXED)
$SIGMAP 0.01 FIXED
$SIGMAPD (96 FIXED)

$DESIGN GROUPSIZE=32 FIMTYPE=1 VARCROSS=1 APPROX=FO NOPRIOR=0

Note that NOPRIOR can be used to turn on or off the prior contribution.

As of nm751, to include a theta that is only involved in the residual variance, and when using FIMDIAG=1, declare that variable “Sigma-like” with GRD=TS() option, as shown in warfarin_pfim23.ctl for the gamma parameter:

$PROBLEM WARFARIN
$INPUT ID TIME AMT RATE EVID MDV DV TSTRAT TMIN TMAX
$DATA warfarin_pfim2.csv ignore=C

$SUBROUTINES ADVAN2 TRANS2
$PK
MU_1=LOG(THETA(1))
MU_2=LOG(THETA(2))
MU_3=LOG(THETA(3))
;MU_4=LOG(THETA(4))
CL=EXP(MU_1+ETA(1))
V=EXP(MU_2+ETA(2))
KA=EXP(MU_3+ETA(3))
gamma=THETA(4)
S2=V
F1=1.0

$ERROR
IPRED=A(2)/V
Y=IPRED + (IPRED**gamma)*EPS(1)

$THETA
0.15 ;[CL]
8.0 ;[V]
1.0 ;[KA]
1.0 ;[gamma]

$OMEGA (0.07) (0.02) (0.6)
$SIGMA 0.01
$DESIGN GROUPSIZE=32 FIMTYPE=1 VARCROSS=1 APPROX=FO

$EST GRD=TS(4) SIGL=10

The standard errors to gamma and sigma, however are somewhat inaccurate (by about 25%, when compared with the equivalent clinical trial simulations warfarin_pfim23_cts.ctl (FOCEI) and warfarin_pfim23_ctsi.ctl (IMP)). Therefore, avoid incorporating a theta parameter in the residual variance assessment when performing optimal designs and evaluations.

Serial correlation models may also be implemented. According to Nyberg et . al (Nyberg, Höglund, et al. 2012), incorporating serial correlation can account for possible model mis-specification, and allow a wider distribution of times among the samples to be found, in order to make the sampling more robust. The following example, warfarin_fobc.ctl shows how this may be implemented (as of nm751, FIMDIAG=1 may be used with a serial correlation problem).

$PROBLEM WARFARIN
$INPUT ID TIME AMT RATE EVID MDV DV TSTRAT TMIN TMAX REPZ
$DATA warfarinz.csv ignore=@
$SUBROUTINES ADVAN2 TRANS2
$ABBR DECLARE T(NO), INTEGER I, DOWHILE J
$PK
MU_1=LOG(THETA(1))
MU_2=LOG(THETA(2))
MU_3=LOG(THETA(3))
CL=THETA(1)*EXP(ETA(1))
V=THETA(2)*EXP(ETA(2))
KA=THETA(3)*EXP(ETA(3))
S2=V
F1=1.0
$ERROR

 IF(NEWIND.NE.2) I=0
 IF (NEWL2 ==1 .AND. MDV==0) THEN

       I=I+1
       T(I)=TIME
       J=1
       DO WHILE (J<=I)
       CORRL2(J,1)=EXP(-THETA(4)*(TIME-T(J)))
       J=J+1
       ENDDO
 ENDIF
IPRED=A(2)/V
Y=IPRED*(1.0+EPS(1)) + EPS(2)
$THETA
0.15 ;[CL]
8.0 ;[V]
1.0 ;[KA]
(0.1) ; [COR]
$OMEGA (0.07 ) (0.02 ) (0.6 )
$SIGMA (0.01 ) (0.001 FIX )
$DESIGN GROUPSIZE=100 FIMDIAG=1 APPROX=FO

                MAXEVAL=9999 SIGL=12 DESEL=TIME DESELSTRAT=TSTRAT DESELMIN=TMIN
                DESELMAX=TMAX
$EST GRD=TS(4)
$TABLE ID TIME AMT EVID MDV DV IPRED NOPRINT NOAPPEND FILE=warfarin_fobc.tab

This results in 7-8 distinct time points (warfarin_fobc.tab). The analysis of real data itself may not require the CORRL2 algorithm, it is mainly designed to create a greater spread of distinct time samples for robustness, as mentioned earlier.

This above example requires some decision about the de-correlation rate THETA(4), which could be arbitrary. Also, it might be more desirable to consider that the farther apart that data are in terms of value (rather than in time), the less correlated they are, to provide a separation of samples based on the degree of change in concentration. In other words, in regions where the concentration changes rapidly, the samples may be more closely spaced, and in regions where concentration changes less rapidly, samples would be spaced farther apart. Furthermore, relative change in concentration could be most appropriate. For this purpose, the following control stream (.._fobc5.ctl) would be suitable:

$PROBLEM WARFARIN
$INPUT ID TIME AMT RATE EVID MDV DV TSTRAT TMIN TMAX REPZ
$DATA warfarinz2.csv ignore=@
$SUBROUTINES ADVAN2 TRANS2
$ABBR DECLARE T(NO), INTEGER I, DOWHILE J
$PK
MU_1=LOG(THETA(1))
MU_2=LOG(THETA(2))
MU_3=LOG(THETA(3))
CL=THETA(1)*EXP(ETA(1))
V=THETA(2)*EXP(ETA(2))
KA=THETA(3)*EXP(ETA(3))
S2=V
F1=1.0
$ERROR

 IPRED=A(2)/V
 MEASURE=DV
; MEASURE=LOG(DV) ; Alternative measure
 IF(NEWIND.NE.2) I=0
 IF (NEWL2 ==1 .AND. MDV==0) THEN

       I=I+1
       T(I)=MEASURE
       CORRL2(I,1)=1.0
       J=1
       DO WHILE (J<I)
       DVAL1=MEASURE
       DVAL2=T(J)
; 0.999 needed to make sure correlation between two different data
; points is less than 1, even if they are of identical value
       CORRL2(J,1)=0.999*EXP(-THETA(4)* &
SQRT( (DVAL1-DVAL2)*(DVAL1-DVAL2)/(DVAL1*DVAL1+DVAL2*DVAL2) ) )
; Alternative correlation algorithm
; CORRL2(J,1)=0.999*EXP(-THETA(4)*ABS((DVAL1-
DVAL2))/(ABS(DVAL1)+ABS(DVAL2)) )
       J=J+1
       ENDDO
 ENDIF
Y=IPRED*(1.0+EPS(1)) + EPS(2)
$THETA
0.15 ;[CL]
8.0 ;[V]
1.0 ;[KA]
(1.0 FIXED) ; [COR]
$OMEGA (0.07 ) (0.02 ) (0.6 )
$SIGMA (0.01 ) (0.001 FIXED)
$DESIGN GROUPSIZE=100 FIMDIAG=1 APPROX=FOI

                MAXEVAL=99999 PRINT=200 SIGL=9 DESEL=TIME DESELSTRAT=TSTRAT
DESELMIN=TMIN DESELMAX=TMAX
$EST GRD=TS(4)
$TABLE ID TIME AMT EVID MDV DV IPRED NOPRINT NOAPPEND FILE=warfarin_fobc5.tab

The measure now used is not TIME, but DV. The IPRED should not be used as a measure, as NONMEM cannot handle random variables (those with eta associated with them) within the CORRL2 structure. Note that the THETA(4) is fixed to 1, which is very reasonable for any problem, and we have removed the arbitrariness of having to select its value.

Table warfarin_fobc5.tab shows final times of samples that provide a reasonable distribution to cover incremental changes in the predictive values of the PK curve. What is also of value is that the standard errors approximate what one would get with a clinical trial simulation, not using the CORRL2 algorithm (cts_fobc5).

Another way of obtaining a robust range of optimal times is to take into consideration that there is a variation of possible parameters by providing a distribution of possible population PK parameters, and assess the distribution of optimal times that result from this. We may want to determine how differences in sample times may occur, given uncertainty in the population parameters. Thus we use the $SIML TRUE=PRIOR with prior information to randomly generate THETAs (and OMEGAs and SIGMAs, if desired), and optimal time points for each of these theta sets will be obtained. The prior information may be selected based on some understanding of likely variability in THETAs, OMEGAs, and SIGMAs, or may have come from some firm experience of previous studies. We add prior information, and simulation instructions (complete control stream priortrue.ctl):

$PRIOR NWPRI PLEV=0.99
; THETAs are in log units, to ensure 30% variation does not
; produce negative parameters

$THETAP (-1.897 FIXED) (2.079 FIXED) (0.0001 FIXED)

$THETAPV BLOCK(3) ; 30% CV=0.3*0.3=0.09 variance
0.09 FIX
0.0 0.09
0.0 0.0 0.09

$SIM (4442223) SUBPROB=1000 TRUE=PRIOR

$DESIGN NELDER FIMDIAG=1 GROUPSIZE=32 DESEL=TIME DESELSTRAT=TSTRAT
DESELMIN=TMIN DESELMAX=TMAX

              MAXEVAL=4000 SIGL=10 nohabort PRINT=100 NOPRIOR=1

$TABLE ID LTVCL LTVV LTVKA TIME EVID MDV DV IPRED NOPRINT NOAPPEND

FILE=priortrue.tab

In this example, 1000 sub-problems are generated, each with their own theta sets. The variance on the THETAs is about 30%CV, and by observing the optimal times listed in priortrue.tab, we get a sense of what kind of differences of optimal time points will be obtained, based on the likely variability of the THETAs, in this example. A FORTRAN program, ...f90, is available, to provide a convenient summarization of the results in priortrue.tab (STD=standard deviation, RSTD=relative standard deviation, or %STD/MEAN):

TIME:       Mean      STD          RSTD              2.50%      97.50%
Row          0.0000   0.0000       0.0000         0.0000        0.0000
            0.13934   5.09463E-02  36.563         6.32440E-02  0.26127
         1   7.1222   1.7673       24.814         4.0921        11.134
         2   7.1222   1.7673       24.814         4.0921        11.134
         3   7.1222   1.7673       24.814         4.0921        11.134
         4   155.81   37.671       24.177         82.548        199.97
         5
         6                                              2.50%   97.50%
                                                      0.0000    0.0000
IPRED:       Mean      STD         RSTD             0.84534     1.3868
Row           0.0000    0.0000     0.0000             4.5215    12.773
              1.0839  0.13800      12.732             4.5215    12.773
         1    7.9587    2.1256     26.708             4.5215    12.773
         2    7.9587    2.1256     26.708          0.37273      1.1902
         3    7.9587    2.1256     26.708
         4  0.52453   0.20825      39.702
         5
         6

The full output from summary.exe contains many percentiles, to allow for a large selection of dispersed time samples. The syntax of executing this program is as follows (the third argument, a FORTRAN style format, is optional):

summary priortrue.tab summary.tab "1X,1PG15.8"

For this example, we need only 2 time samples for dispersing the clustered set near 7.0 hours, so we use 2.5 and 97.5 percentile positions, to obtain coverage of the possible THETAs that may occur. Of particular interest for us is to take some sense of lower and upper bounds around the 7 hour time, such as the 2.5% and 97.5% values of the cluster points (whose mean is 7.122, thus representing the 7 hours). In this case the values are 4.09 and 11.13 that we obtain from the above summary table, and we can use the floor of 4.09 (->4) and ceiling of 11.13 (->12) and use these as recommended times around the 7.0 hour clustered time to improve robustness.

We may also wish to make the theta variance equal to the OMEGA inter-subject variance, to take into account likely fluctuations of PK profiles among the subjects, and make the sampling robust for between subject variability. We would thus make the $THETAPV variances equivalent to the $OMEGA variances (priortrue2.ctl):

$THETAPV BLOCK(3) ; make variance same as that of $OMEGA
0.07 FIX
0.0 0.02
0.0 0.0 0.6

This process of determining what variability of THETAs could provide what degree of variability of ideal sample times (or other design elements) provides a robust means of assessing sensitivity between parameters and design elements.

16 Modifying the Weighting of Components of LOG(DET(FIM)) for use in Robust Design

As of nm760, you can modify the weighting of components of the LOG(DET(FIM)), for OFVTYPE=0,1,6, and 8, for use in robust design. Please read reference (Fayette et al. 2023) for details about robust design. Example ..5.ctl codes four possible models, in accordance with (Fayette et al. 2023), and the subroutine .._wfims.f90 provides the appropriate weights for the various parts of the (DET(FIM)) in accordance with equations 5 and 9 of (Fayette et al. 2023). It is worth while to read the comments in wfims.f90, and (Fayette et al. 2023), to understand how to set the weight information for various stages of the robust design process. The template ...f90 is available to modify for robust design problems in general. If WFIMS.f90 is used unmodified (with no executable code, as a dummy plug-in), the weights are all set to 1 for the default, non-robust design condition.

The relevant files for this robust design example are as follows:

Robust5.ctl: Stage 1, robust design optimization for 4 models of doses, dose group proportions, and sample times of each dose group. The 4 models are linear, log-linear, ED50 (Emax) model, and ED50s model (Emax with sigmoid coefficient=2).

Robust5_sim.ctl: The result in robust5.tab is used to manually create a simulation template file called robust5_sim.dat, from which simulated data are created using ed50s model as the “true” model. Simulated data are placed in ed50s_fit.dat.

linear_fit.ctl: The linear model fitting ed50s_fit.dat as the data from a “clinical trial”.

llinear_fit.ctl: The log-linear model fitting ed50s_fit.dat as the data from a “clinical trial”.

Ed50_fit.ctl: The Emax model fitting ed50s_fit.dat as the data from a “clinical trial”.

Ed50s_fit.ctl: The Emax model with sigmoid fitting ed50s_fit.dat as the data from a “clinical trial”.

Weights are calculated from the OFV’s of the four model fits, using exp(-(OFV/2+p)) as the weights, where p=number of fitted parameters to the model, in accordance with (Fayette et al. 2023).

Robust5b.ctl: Stage 2, the weights are used to set WT() array in robust5b_wfims.f90, used by robust5b.ctl in the relative weighting of the various model sections in the FIM. In this second stage of robust design, the weights favor the ed50 models, in that these gave the best fits to the “clinical trial” data. The final doses (DOSE), fractional representation of dose groups (STRATF), and optimized sample times of each dose can be used for the design of the second clinical trial.