Dear NONMEM users,
A few month ago, a post by Li Li ("Question regarding Calculation Process in $DES BLOCK") generated some interesting discussions about how to implement abrupt changes in a model parameter in $DES. The original problem proposed by Li Li was to implement a step function that turns on and off between two time points reset on a 24h clock time basis. Alison Boeckmann strongly recommended NOT to use IF statements or create a 24-h time variable using the INT function in $DES. Instead she suggested to use model time events (MTIME's). My question relates to the implementation of such a solution.
The attachments to this post are simulation control stream illustrating some attempts to use MTIME's to switch on and off a FLAG variable every 24h. On a 24h clock time, the flag switches from 0 to 1 at time SHIFT (the so-called shift from time 0, bound between 0 and 24) and turns back to 0 at time SHIFT+DUR (DUR being the duration of the step function, bound between 0 and 24). Two MTIME's variables define the "switch" times, and the flag is set as follows: FLAG=MPAST(1)-MPAST(2). This would work nice and well if I did not need to reset the MTIME every 24h.
In the attached model, MTIME's depend on the TX variable which is incremented by 24 based upon a comparison of the TIME variable with various equations of TX and DUR. As documented in the manual, the MTDIFF reserved variable is set to 1 to reflect the changes in MTIME variables. Furthermore, $DES defines 3 states to monitor MPAST(1), MPAST(2) and FLAG. The integration of the third state over time is expected to increase by DUR every 24h. But, this is not always the case depending on the dataset records, the values of SHIFT, DUR, and the methods of MTIME update:
circadian1.ctl: MTIME's are coded to be reset when MTIME(1)+24 is reached. Results: MTIME's are never reset if SHIFT+DUR<24 and, apparently, if there is no dataset record between MTIME(2) and 0.5*(MTIME(2) + MTIME(1)+24) circadian2.ctl: MTIME's are coded to be reset when MTIME(1)+DUR is reached. Results: MTIME's are never properly set if SHIFT+DUR<24 circadian3.ctl: a third MTIME is used and MTIME's are coded to be reset at the midpoint between MTIME(2) and MTIME(1)+24 (ie in between the steps). Results: if SHIFT+DUR<24, MTIME's are never properly reset if there is no dataset record between MTIME(2) and the midpoint.
For all methods, the integration is "wrong" when there is a gap of more than 24h between dataset records (try running any models with exclusion of the EVID=2 records). Additionally, because SHIFT and DUR will be estimated for each subject in my real model, adding EVID=2 records in the dataset does not seem to be a reasonable option to apply the "circadian3.ctl" method.
Basically, I am stuck and would greatly appreciate any code or dataset suggestions to fix this integration problem. Obviously, any alternative approach to implement this circadian switch function would also be welcome.
Thank you for your time
Sebastien
circadian.csv
Description:
MS-Excel spreadsheet
$PROBLEM circadian1
;; -----------------------------------------------------------------------------
;; PURPOSE: Circadian rhythm - update of MTIME's at MTIME(1)+24
;; -----------------------------------------------------------------------------
$INPUT ID TIME MYTIME AMT CMT EVID DV MDV
$DATA circadian.csv IGNORE=@
$THETA
1
1
$OMEGA
0 FIX
$SIGMA
0.1
$SUBROUTINE ADVAN13 TOL=9
$MODEL NCOMP=3
COMP=(MP1)
COMP=(MP2)
COMP=(FLAG)
$PK
SHIFT=THETA(1)*EXP(ETA(1))
DUR=THETA(2)
IF (ID.GT.0.5) THEN
SHIFT = 4.2
DUR = 7.1
END IF
IF (ID.GT.1.5) THEN
SHIFT = 3.3
DUR = 17.6
END IF
IF (ID.GT.2.5) THEN
SHIFT = 22.7
DUR = 8.3
END IF
IF (ID.GT.3.5) THEN
SHIFT = 14.5
DUR = 17.9
END IF
;Circadian rhythm
IF ((SHIFT+DUR).GT.24) THEN
OFFSET = 24
ELSE
OFFSET = 0
END IF
OFFSETI = SHIFT - OFFSET
IF (NEWIND.LE.1) TX=0
IF (TIME.EQ.0) THEN
TX = SHIFT-OFFSET
END IF
TX = TX
IF (TIME.GE.(TX+24)) TX=TX+24
MTIME(1)=TX
MTIME(2)=TX+DUR
MTIME(3)=(TX+DUR+24+TX)/2
MT1 = MTIME(1)
MT2 = MTIME(2)
IF(NEWIND.LE.1) MTDIFF=0
IF(NEWIND.GT.1) MTDIFF=1
$DES
FLG=MPAST(1)-MPAST(2)
DADT(1) = MPAST(1)
DADT(2) = MPAST(2)
DADT(3) = FLG
$ERROR
IPRED = A(3)
MPAST1 = A(2)
MPAST2 = A(3)
Y = IPRED*(1+EPS(1))
$SIMULATION(12345) NSUB=1 ONLYSIM
$TABLE ID TIME SHIFT DUR TX IPRED MPAST1 MPAST2 NOPRINT ONEHEADER
FILE=circadian1.tbl
$PROBLEM circadian2
;; -----------------------------------------------------------------------------
;; PURPOSE: Circadian rhythm - update of MTIME's at MTIME(1)+DUR
;; -----------------------------------------------------------------------------
$INPUT ID TIME MYTIME AMT CMT EVID DV MDV
$DATA circadian.csv IGNORE=@
$THETA
1
1
$OMEGA
0 FIX
$SIGMA
0.1
$SUBROUTINE ADVAN13 TOL=9
$MODEL NCOMP=3
COMP=(MP1)
COMP=(MP2)
COMP=(FLAG)
$PK
SHIFT=THETA(1)*EXP(ETA(1))
DUR=THETA(2)
IF (ID.GT.0.5) THEN
SHIFT = 4.2
DUR = 7.1
END IF
IF (ID.GT.1.5) THEN
SHIFT = 3.3
DUR = 17.6
END IF
IF (ID.GT.2.5) THEN
SHIFT = 22.7
DUR = 8.3
END IF
IF (ID.GT.3.5) THEN
SHIFT = 14.5
DUR = 17.9
END IF
;Circadian rhythm
IF ((SHIFT+DUR).GT.24) THEN
OFFSET = 24
ELSE
OFFSET = 0
END IF
OFFSETI = SHIFT - OFFSET
IF (NEWIND.LE.1) TX=0
IF (TIME.EQ.0) THEN
TX = SHIFT-OFFSET
END IF
TX = TX
IF (TIME.GE.(TX+DUR)) TX=TX+24
MTIME(1)=TX
MTIME(2)=TX+DUR
MTIME(3)=(TX+DUR+24+TX)/2
MT1 = MTIME(1)
MT2 = MTIME(2)
IF(NEWIND.LE.1) MTDIFF=0
IF(NEWIND.GT.1) MTDIFF=1
$DES
FLG=MPAST(1)-MPAST(2)
DADT(1) = MPAST(1)
DADT(2) = MPAST(2)
DADT(3) = FLG
$ERROR
IPRED = A(3)
MPAST1 = A(2)
MPAST2 = A(3)
Y = IPRED*(1+EPS(1))
$SIMULATION(12345) NSUB=1 ONLYSIM
$TABLE ID TIME SHIFT DUR TX IPRED MPAST1 MPAST2 NOPRINT ONEHEADER
FILE=circadian2.tbl
$PROBLEM circadian3
;; -----------------------------------------------------------------------------
;; PURPOSE: Circadian rhythm - update of MTIME's at mid point between MTIME(2)
and MTIME(1)+24
;; -----------------------------------------------------------------------------
$INPUT ID TIME MYTIME AMT CMT EVID DV MDV
$DATA circadian.csv IGNORE=@
$THETA
1
1
$OMEGA
0 FIX
$SIGMA
0.1
$SUBROUTINE ADVAN13 TOL=9
$MODEL NCOMP=3
COMP=(MP1)
COMP=(MP2)
COMP=(FLAG)
$PK
SHIFT=THETA(1)*EXP(ETA(1))
DUR=THETA(2)
IF (ID.GT.0.5) THEN
SHIFT = 4.2
DUR = 7.1
END IF
IF (ID.GT.1.5) THEN
SHIFT = 3.3
DUR = 17.6
END IF
IF (ID.GT.2.5) THEN
SHIFT = 22.7
DUR = 8.3
END IF
IF (ID.GT.3.5) THEN
SHIFT = 14.5
DUR = 17.9
END IF
;Circadian rhythm
IF ((SHIFT+DUR).GT.24) THEN
OFFSET = 24
ELSE
OFFSET = 0
END IF
OFFSETI = SHIFT - OFFSET
IF (NEWIND.LE.1) TX=0
IF (TIME.EQ.0) THEN
TX = SHIFT-OFFSET
END IF
TX = TX
IF (TIME.GE.((TX+DUR+24+TX)/2)) TX=TX+24
MTIME(1)=TX
MTIME(2)=TX+DUR
MTIME(3)=(TX+DUR+24+TX)/2
MT1 = MTIME(1)
MT2 = MTIME(2)
IF(NEWIND.LE.1) MTDIFF=0
IF(NEWIND.GT.1) MTDIFF=1
$DES
FLG=MPAST(1)-MPAST(2)
DADT(1) = MPAST(1)
DADT(2) = MPAST(2)
DADT(3) = FLG
$ERROR
IPRED = A(3)
MPAST1 = A(2)
MPAST2 = A(3)
Y = IPRED*(1+EPS(1))
$SIMULATION(12345) NSUB=1 ONLYSIM
$TABLE ID TIME SHIFT DUR TX IPRED MPAST1 MPAST2 NOPRINT ONEHEADER
FILE=circadian3.tbl