97  Stochastic Differential Equation Plug-In (NM72)

Introduction to NONMEM 7.6.0

Author

Robert J. Bauer

Published

June 1, 2025

An alternative method to evaluating stochastic differential equation problems is to utilize the plug-in routine SDE.f90 in the NONMEM ..\examples\sde directory, which numerically evaluates the SDE equations, without requiring in-line coding into the control stream. An example control stream file is as follows (..\examples\sde\sde_ex2_impo.ctl):

; Based on sde_ex2_base.ctl, using SDE.f90 OTHER file, and .dat file
; modified. From Christoffer Tornoe, example 2

$PROBLEM PK ODE HANDS ON ONE

$ABBR DES=FULL ; Must have this, so DA array is simple-mapped.

$INPUT ID TIME DV AMT CMT FLAG MDV SDE ; Add SDE data item. Should have
; values between 0 and 5, as described in comments of SDE.f90
; The data item SDE must be added. It has values of 0,1,2,3,4, or 5:
; SDE (SDE data item)=0, BEGINNING OF SUBJECT, OR WHENEVER YOU WANT TO
; INITIALIZE SDE SYSTEM (SOMETIMES YOU WANT TO DO THIS
; EVEN WITHIN A SUBJECT, LIKE A NEW OCCASION)
; SDE=1, FIRST OBSERVATION OF PRESENT TIME.
; SDE=2 LAST OBSERVATION FOR PRESENT TIME.
; SDE=5 FIRST AND LAST OBSERVATION FOR PRESENT TIME.
; SDE=6, MIDDLE OBSERVAION FOR PRESENT TIME.
; THUS:
; DOSE RECORD, TIME=0, THEN SDE=0
; PK OBSERVATION, TIME=0.1, SDE=1
; PD OBSERVATION, TIME=0.1, SDE=6
; EFFICACY OBSERVATION, TIME=0.1, SDE=2 (SDE=2 SINCE NEXT
; RECORD HAS A NEW TIME)
; PK OBSERVATION, TIME=0.2, SDE=1
; PD OBSERVATION, TIME=0.2, SDE=2
; PK OBSERVATION, TIME=0.5, SDE=5
; PK OBSERVATION, TIME=1.0, SDE=5

$DATA sde_ex2o.dat
              IGNORE=@

$SUBROUTINE ADVAN6 TOL=9 DP OTHER=SDE.f90
; nde=number of base equations, ncmt=number of observation compartments
$ABBR DECLARE SGW(3) ; need at least nde of these

$MODEL
   COMP = (CENTRAL); there are nde base states from original sde_ex2_base.ctl
   COMP = (DFDX1) ; need to add ncmt observation compartments for SDE
   COMP = (DPDT11) ; Will need (nde+1)*nde/2 of these for SDE

$PK
   IF(NEWIND.NE.2) OT = 0
   MU_1 = THETA(1)
   CL = EXP(MU_1+ETA(1))
   MU_2 = THETA(2)
   VD = EXP(MU_2+ETA(2))
   SGW1 = THETA(4) ; Add estimable scalar for modeling the SDE noise.
   NCMT=1.0 ; number of compartments. Added for calls to SDE_DER and SDE_CADD

 NDE=1.0 ;Number of original, base ODEs. Added for calls to SDE_DER, SDE_CADD

$DES
 FIRSTEM=1 ; MAke sure FIRSTEM=1 so that DA arrays (Derivatives of DADT() wrt
; A()), are calculated, even when IMP is done.
 DADT(1) = - CL/VD*A(1) ; Original base derivative from sde_ex2_base.ctl

;NEXT DERIVATIVES ARE ACTUALLY PREDICTIVE VALUES FOR CMPTS
; 1 & 2, RESPECTIVELY
; Derivatives of these with respect to A() will be calculated
; symbolically by DES routine created by NMTRAN

 DADT(2) = A(1)/VD ; Add output equations, required for each CMT value.
; DUMMY PLACEMENT FOR DERIVATIVES OF THE STOCHASTIC ERROR SYSTEM. THESE ARE
; FILLED OUT BY SDE_DER
SGW(1)=SGW1 ; Specify SGW with appropriate index , for appropriate
; DES equation number
; the DA() array THEN contains all derivatives of DADT (=DXDT) with
; respect to A(=X).
; number of base model derivative equations (nde)=1,
; Number of compartments (ncmt)=1.
; DA is a reserved array, dimensioned DA(IR,*)

"LAST
"         CALL SDE_DER(DADT,A,DA,IR,SGW,NDE,NCMT)

$ERROR (OBS ONLY)
       IPRED = A(1)/VD
       IRES = DV - IPRED
       W    = THETA(3)
       IWRES = IRES/W
       WS=1000.0
; CENTRAL COMPARTMENT, PLASMA LEVELS
; EPS(1) = USER MODEL ERROR CONTRIBUTION
; EPS(2) = STOCHASTIC ERROR CONTRIBUTION. THE WS IS JUST A
; PLACEHOLDER COEFFICIENT. SDE_CADD WILL REPLACE THIS
; WITH THE CORRECT VALUE
       Y    = IPRED+W*EPS(1) + WS*EPS(2)
; SDE_CADD WILL EVALUATE THE TRUE COEFFICIENTS (WS) TO THE
; STOCHASTIC COMPONENTS.
; In general, if you have nmcmt observation compartments, then first ncmt
; EPS() will pertain to
; measurement error, and the second ncmt set of EPS() will pertain
; to stochastic errors.
; This means you cannot have L2 type correlations, and prop+additive
; should be packaged into a single EPS().
; For two observations, you may have:
; IF(CMT==1) THEN
; IPRED=A(1)/V
; W=SQRT(THETA((5)*THETA(5)*IPED*IPRED+THETA(6)*THETA(6))
; Y=IPRED+W*EPS(1)+WS*EPS(3)
; ENDIF
; IF(CMT==2) THEN
; IPRED=A(2)/V
; W=SQRT(THETA((7)*THETA(7)*IPED*IPRED+THETA(8)*THETA(8))
; Y=IPRED+W*EPS(2)+WS*EPS(4)
; ENDIF
; Number of compartments=1, number of base model derivative equations=1

"LAST
"         CALL SDE_CADD(A,HH,TIME,DV,CMT,NDE,NCMT,SDE)

$THETA (0,2.3)                ;1 CL
$THETA (0,3.5)                ;2 VD
$THETA (0, 2)                ;4 SIGMA
$THETA (0,1) ; SGW1

$OMEGA 0.1                   ;1 CL
$OMEGA 0.01                  ;2 VD
$SIGMA (1 FIX) (1 FIX)
                                       ; PK

$EST METHOD=ITS INTERACTION LAPLACE NUMERICAL SLOW NOABORT PRINT=1 CTYPE=3
SIGL=5
$EST METHOD=IMP INTERACTION NOABORT SIGL=5 PRINT=1 IACCEPT=1.0 CTYPE=3
$EST MAXEVAL=9999 METHOD=1 LAPLACE INTER NOABORT NUMERICAL SLOW NSIG=3
PRINT=1 SIGL=9

$COV MATRIX=R UNCONDITIONAL

$TABLE ID TIME FLAG AMT CMT IPRED IRES IWRES
            ONEHEADER NOPRINT FILE=sde_ex2_impo.tab

Please note that using this add-in requires $ABBR DES=FULL setting. This process works well with the methods such as importance sampling, SAEM, or BAYES, but works only partially for classical NONMEM methods or ITS. If using with classical NONMEM methods or ITS, it is better to set LAPLACE NUMERICAL, although it does not solve the problem perfectly. Classical methods rely on NMTRAN creating symbolic derivatives of the residual variance components with respect to eta, which is used to create the proper individual objective function. For this to occur, NMTRAN has to see all of the relevant equations in the control stream file, or the user must have the eta derivatives evaluated. This method has some of the SDE differential equations and RVAR components calculated in subroutines SDE_DER and SDE_CADD, “hidden” from NMTRAN. Despite this problem, classical NONMEM methods provide parameters using the SDE call routines that are similar, although not identical, to those when the SDE equations are placed in-line into the control stream file. To see how the SDE call routines work for each of the analysis methods, see sde_ex2_impo.res that uses SDE.f90, and compare the results with sde_ex2_imp_mu.res, which uses the in-line equations. The new methods (except ITS) do not need these NMTRAN constructed components, so they work with the SDE call routines quite well.

As of NM73, numerical eta derivatives are now available for FOCE/ITS, so that it is not necessary for NMTRAN to see all the code, or for the user to supply evaluation of the eta derivatives. In the following example, OPTMAP=1 is chosen to provide forward finite difference eta derivatives for the search, and ETADER=2 is chosen to provide numerically assessed central finite difference derivatives to the Hessian matrix of the posterior density (sde_ex2_impo2.ctl), allowing ITS and FOCE to obtain results similar to Importance sampling:

$EST METHOD=ITS INTERACTION NOABORT PRINT=1 CTYPE=3 OPTMAP=1 ETADER=2 SIGLO=6
 SIGL=6 MCETA=1

$EST METHOD=IMP INTERACTION NOABORT PRINT=1 IACCEPT=1.0 CTYPE=3 OPTMAP=0
 ETADER=0 SIGLO=6 SIGL=6 MCETA=1 MAPITER=0

$EST MAXEVAL=9999 METHOD=1 INTER NOABORT NSIG=1 PRINT=1
 OPTMAP=1 ETADER=2 SIGLO=6 SIGL=6 MCETA=1 SLOW

$COV MATRIX=R UNCONDITIONAL TOL=9 SIGL=8 SIGLO=8
$TABLE ID TIME FLAG AMT CMT IPRED IRES IWRES
            ONEHEADER NOPRINT FILE=sde_ex2_impo2.tab

Example 4 from Christoffer Tornoe’s workshop has also been processed in several ways. They are as follows:

sde_ex4_base.ctl: User defined control stream, without SDE components.

sde_ex4_foce.ctl: SDE components added to sde_ex4_base.ctl, by Tornoe’s R/Splus command, for example:

source("NONMEMSDEscript.R")
SDEmodel( CS="sde_ex4_base.ctl", datafile="sde_ex4_base.dat")

The output files created are

sde_ex4_base.ctlSDE
sde_ex4_base.datSDE

These were manually renamed to

sde_ex4_foce.ctl
sde_ex4.dat

for the present example.

sde_ex4_foceo.ctl: This problem was modified to use the SDE.f90 OTHER routine, and FOCE was used for estimation.

sde_ex4_foceo.ctl: This problem was modified to use the SDE.f90 OTHER routine, and FOCE was used for estimation.

sde_ex4_impo.ctl: This problem was modified to use the SDE.f90 OTHER routine, and IMP was used for estimation.

See ..\examples\sde_inline\readme.txt in using Tornoe’s inline code enhancer method. The choice to use Tornoe’s inline code enhancer (using NONMEMSDEscript.R) or whether to use the SDE.f90 OTHER routine method depends on the complexity of the model. The NONMEMSDEscript can only enhance NMTRAN code for which there are observations from only one output compartment. Its advantage is the NMTRAN sees all of the ODE code, and provides appropriate analytical derivative components for random variables (ETAS). Typically this results in faster execution, especially for FOCE. The SDE.f90 OTHER routine can handle multiple output compartments, but does not provide derivatives to the additional ODE components, so finite difference evaluation of random eta variables (through OPTMAP=1 and ETADER=2) must be used, which can be 2-4 times slower in execution.