83  $ANNEAL to facilitate EM search methods (NM73)

Introduction to NONMEM 7.6.0

Syntax:

$ANNEAL number-list1:value1 number-list2:value2
etc. for as many lists that are needed.

Example:

$ANNEAL 1-3,5:0.3 6,7:1.0

Sets starting diagonal Omega values for purposes of simulated annealing. Thus, initial values of OMEGA(1,1), OMEGA(2,2), OMEGA(3,3), and OMEGA(5,5) are set to 0.3, while initial OMEGA(6,6) and OMEGA(7,7) are set to 1.0. When $EST CONSTRAIN>=4, an algorithm in constraint.f90 will initially set the OMEGAs to these values, and then shrink these OMEGA values more and more with each iteration, and eventually shrinks the OMEGAs to 0, the intended target value for that Omega. This is a technique that may be used especially with SAEM, to provide an annealing method for moving THETAs that have 0 omega values associated with them. The default is the use of gradient methods, which are good for problems starting near the solution, whereas the annealing method is more suitable for problems starting far from the solution.

An example is anneal.ctl, an EMAX model in which the Hill coefficient does not have inter-subject variance (that is, its omega variance is set to 0):

$PROB Emax model with hill=3
$INPUT ID DOSE DV
$DATA anneal.dat IGNORE=@
$PRED

MU_1 = THETA(1)
EMAX = EXP(MU_1+ETA(1))
MU_2 = THETA(2)
ED50 = EXP(MU_2+ETA(2))
MU_3 = THETA(4)
E0 = EXP(MU_3+ETA(3))

MU_4=THETA(3)
HILL = EXP(MU_4+ETA(4))

IPRED = E0+EMAX*DOSE**HILL/(ED50**HILL+DOSE**HILL)

Y  = IPRED + EPS(1)

$THETA 4.1 ; 1. Emax
$THETA 6.9 ; 2. ED50
$THETA 0.001 ; 3. Hill
$THETA 2.3 ; 4. E0

$OMEGA BLOCK(2) 0.1
                              0.01 0.1

$OMEGA 0.1
$OMEGA 0.0 FIXED

$ANNEAL 4:0.3

$SIGMA 1
$ESTIMATION METH=SAEM INTER NBURN=1000 NITER=500 ISAMPLE=5 IACCEPT=0.3 CINTERVAL=25 CTYPE=0
NOABORT PRINT=50 CONSTRAIN=5 SIGL=8
$ESTIMATION METH=IMP INTER PRINT=1 NITER=0 ISAMPLE=10000 EONLY=1 CONSTRAIN=0 MAPITER=0 DF=4
$COV MATRIX=R UNCONDITIONAL

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. The subroutine CONSTRAINT may also be used to provide any kind of constraint pattern on any parameters.

Another technique is to use an initial Monte Carlo search method using $EST METHOD=CHAIN ISAMPEND, and then use the standard gradient method for SAEM, as follows:

$PROB Emax model with hill=3
$INPUT ID DOSE DV
$DATA anneal.dat IGNORE=@
$PRED

MU_1 = THETA(1)
EMAX = EXP(MU_1+ETA(1))
MU_2 = THETA(2)
ED50 = EXP(MU_2+ETA(2))
MU_3 = THETA(4)
E0 = EXP(MU_3+ETA(3))

MU_4=THETA(3)
HILL = EXP(MU_4+ETA(4))

IPRED = E0+EMAX*DOSE**HILL/(ED50**HILL+DOSE**HILL)

Y       = IPRED + EPS(1)

$THETA  4.1 ; 1. Emax
$THETA  6.9 ; 2. ED50
$THETA  (-3.0,0.001,3.0) ; 3. Hill
$THETA  2.3 ; 4. E0

$OMEGA BLOCK(2) 0.1
                              0.01 0.1

$OMEGA 0.1
$OMEGA 0.0 FIXED

$SIGMA 1
$EST METHOD=CHAIN ISAMPLE=1 ISAMPEND=30 NSAMPLE=30 FILE=anneal2.chn
$ESTIMATION METH=SAEM INTER NBURN=4000 NITER=200 ISAMPLE=5 IACCEPT=0.3 CINTERVAL=25 CTYPE=3
NOABORT PRINT=100
$ESTIMATION METH=IMP INTER PRINT=1 NITER=0 ISAMPLE=10000 EONLY=1 MAPITER=0
$COV MATRIX=R UNCONDITIONAL

Notice that the range of Monte Carlo search for the Hill coefficient is from -3 to 3, the specified lower and upper bound values (note that theta(3) is actually the log of the Hill coefficient). See I.66 Method for creating several instances for a problem starting at different randomized initial positions: $EST METHOD=CHAIN and $CHAIN Records.

For importance sampling, the CONSTRAINT subroutine is set up so that when CONSTRAIN=4, it uses NBURN as the first NBURN iterations of the total NITER to be for the annealing process, and the remaining NITER-NBURN iterations will continue the problem with the annealing OMEGA being set to 0. See ..\examples\anneal4.ctl, and notice that NBURN=40, NITER=60, and CONSTRAIN=4.