Follow up: Model building advice
From: "Eleveld, DJ" d.j.eleveld@anest.umcg.nl
Subject: [NMusers] Follow up: Model building advice
Date: Thu, May 19, 2005 7:01 am
Thanks everyone for the model building advice. I implemented Katya's advice and
used COM(1) for communication between $ERROR and $PK to set the appropriate bioavailibility
and added dose records at the appropriate times in the datafile. This almost doubled
the amount of records to about 6100 and as expected also increased runtimes. The code
looks clean and simulations are visually quite reasonable. But I still cant get any
estimation methods to converge. Using the FO method to start here is a brief summary
my results so far:
Using ADVAN6:
-------------
0MINIMIZATION TERMINATED
DUE TO MAX. NO. OF FUNCTION EVALUATIONS EXCEEDED
NO. OF FUNCTION EVALUATIONS USED:10024
NO. OF SIG. DIGITS UNREPORTABLE
The resulting theta values are completly unreasonable so I cant use them to try and
improve the inital theta values.
Using ADVAN9:
-------------
ERROR IN LSODI1: CODE 204
0PROGRAM TERMINATED BY OBJ
MESSAGE ISSUED FROM ESTIMATION STEP
Using ADVAN8:
-------------
Nonmem has been running for more than 48 hours now and still no results.
This is too long for a FO method? Should I assume that something has gone wrong?
Here is my control file:
------------------------
$PROB Potentiation fitting
$DATA potent.da2
$INPUT ID TIME DV CFLG MDV AMT RATE CMT
$SUBROUTINES ADVAN9 TOL=6
$ABBREVIATED COMRES=1
$MODEL NCOMPARTMENTS=5
NPARAMETERS=12
COMP(CENTRAL DEFDOSE)
COMP(PERIF1 NOOFF NODOSE)
COMP(PERIF2 NOOFF NODOSE)
COMP(EFFECT NOOFF NODOSE)
COMP(POTENT NOOFF)
$PK
V1=THETA(1)*EXP(ETA(1))
V2=THETA(2)*EXP(ETA(2))
V3=THETA(3)*EXP(ETA(3))
CL=THETA(4)*EXP(ETA(4))
Q3=THETA(6)*EXP(ETA(6))
Q2=(THETA(2)*(THETA(6)/THETA(3) + THETA(5)))*EXP(ETA(5))
S1=V1
KEO=THETA(7)*EXP(ETA(7))
EC50=THETA(8)*EXP(ETA(8))
GAMM=THETA(9)*EXP(ETA(9))
POTR=ABS(THETA(10)+ETA(10)) ; r has normal dist
POTK=THETA(11)*EXP(ETA(11))
SCAL=THETA(12)+ETA(12) ; scale has normal dist
IF (NEWIND.EQ.0) COM(1)=0
F5=POTR*(1-COM(1))
$DES
C1=A(1)/V1 ; Conc in V1
C2=A(2)/V2 ; Conc in V2
C3=A(3)/V3 ; Conc in V3
DADT(1)=Q2*(C2-C1)+Q3*(C3-C1)-CL*C1 ; Change in V1
DADT(2)=Q2*(C1-C2) ; Change in V2
DADT(3)=Q3*(C1-C3) ; Change in V3
DADT(4)=(C1-A(4))*KEO ; Effect compartment conc
DADT(5)=-POTK*A(5) ; Decay potentiation
$ERROR
DPOT=A(5) ; The degree of potentiation
IF (DPOT.LT.0) DPOT=0
CEFF=A(4)
IF (CEFF.LT.0) CEFF=0
DPD1=CEFF**GAMM ; Degree of NMB
NMB=DPD1/(DPD1+(EC50**GAMM))
COM(1)=NMB
CPRE=A(1)/V1*EXP(ERR(1)) ; Conc prediction
RPRE=SCAL*(1+DPOT)*(1-NMB)+ERR(2) ; Twitch prediction
IF (DV.EQ.-999) RPRE=-999+ERR(2) ; Missing twitch sizes
Y=CFLG*CPRE+(1-CFLG)*RPRE
$THETA (0,0.042,0.07)(0,0.041,0.08)(0,0.089,0.15)
(0,0.003,0.03)(0,0.005,1)(0,0.001,0.01)
(0.05,0.15,0.3)(1100,1500,1800)(3,4.5,5.5)
(0,0.006,0.006)(0.01,0.10,0.35)(90,100,110)
$OMEGA 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.001 0.04 2.
$SIGMA 0.2 1.
$ESTIMATION MAX=9999 SIG=6
;$SIMULATION (434, 372) (2, 3) ONLY
;$ESTIMATION MAX=9999 SIG=3 METHOD=COND NOABORT POSTHOC
;$ESTIMATION MAX=9999 NOABORT POSTHOC
$TABLE TIME V1 V2 V3 CL Q2 Q3 NOPRINT FILE=potent1.txt
$TABLE TIME KEO EC50 GAMM POTR POTK SCAL NOPRINT FILE=potent2.txt
$TABLE TIME NMB CPRE RPRE DPOT NOPRINT FILE=potent3.txt
Thanks,
Doug Eleveld