Re: AW: Simulations with/without residual error
Mats,
Thanks for extending the code fragment I proposed earlier to add some extra useful features.
I don't know what was incorrect about the code I suggested but resampling the uncertainty ETA values is a pragmatic way to avoid PREDPP errors caused by randomly generated non-positive CL or V values. The code below produces the same results as your code but it shows how to check on reasonable values for CL and V that do not depend on the value of CL or V defined by THETA or other things such as the effect of covariates on the group parameter value.
I note that your example using simulation and estimation does not use the same residual error model for simulation and estimation. The simulation residual error is exp(EPS(1)) but METHOD=CONDITIONAL INTERACTION would approximate this by (1+EPS(1)). You may remember the commentary on this point that Stuart Beal wrote in 2002 :-)
Best wishes,
Nick
Beal SL. Commentary on Significance Levels for Covariate Effects in NONMEM. Journal of Pharmacokinetics & Pharmacodynamics. 2002;29(4):403-10.
$SIM (20090726 NEW) ONLYSIM SUBPROBLEMS=100
; estimates of THETA and OMEGA from previous run
$THETA
1 ; POP_CL theta1
10 ; POP_V theta2
$OMEGA
0.5 ; PPV_CL eta1
0.5 ; PPV_V eta2
;variance-covariance matrix of the THETA estimates from previous run
$OMEGA BLOCK(2)
0.2 ; UNC_POP_CL eta3
0.1 3 ; UNC_POP_V eta 4
$SIGMA .01
$SUB ADVAN1
$PK
; get CL and V uncertainties
IF (ICALL.EQ.4.AND.NEWIND.EQ.0) THEN ; do this just once per subproblem
UNCCL=THETA(1)+ETA(3)
UNCV=THETA(2)+ETA(4)
DO WHILE (UNCCL.LE.0.OR.UNCV.LE.0) ;resample if values unreasonable
CALL SIMETA(ETA)
UNCCL=THETA(1)+ETA(3)
UNCV=THETA(2)+ETA(4)
END DO
ENDIF
CL=UNCCL*EXP(ETA(1)) ; with uncertainty for CL
V =UNCV*EXP(ETA(2)) ; with uncertainty for V
K=CL/V
$ERROR
Y =F*EXP(EPS(1))
Mats Karlsson wrote:
> Dear Nick, Marc and all,
>
> An even later addition to this thread. Nick pointed out the possibility to do simulations in NONMEM that incorporated error in simulation parameters. The code provided gave a good hint at how to do that but would work incorrectly unless some additional NONMEM switches were flicked. Also, as pointed out by Marc and others, one might like to truncate the parameter distribution for simulation (avoid negative CL etc). Also it may be convenient to do simulation followed by re-estimation in a single run. For appropriate post-processing this requires that simulation parameters are output somewhere. A model file, extending Nick’s example, that does these things are provided below.
>
> $PROB
>
> $DATA data1
>
> $INPUT ID DV TIME TVCL TVV AMT
>
> ;data1=
>
> ; 1 0 0 0 0 100
>
> ; 1 0 1 0 0 0
>
> ; 1 0 21 0 0 0
>
> ; 2 0 0 0 0 100
>
> ; etc
>
> $SIM (20090726 NEW) SUBPROBLEMS=100
>
> $THETA ; estimates of THETA and OMEGA from previous run
>
> (0,1) ; POP_CL theta1
>
> (0,10) ; POP_V theta2
>
> $OMEGA
>
> 0.5 ; PPV_CL eta1
>
> 0.5 ; PPV_V eta2
>
> ;variance-covariance matrix of the THETA estimates from previous run
>
> $OMEGA BLOCK(2)
>
> 0.2 ; UNC_POP_CL eta3
>
> 0.1 3 FIX; UNC_POP_V eta 4
>
> $SIGMA .01
>
> $SUB ADVAN1
>
> $PK
>
> CL=THETA(1)*EXP(ETA(1)) ; with uncertainty for CL
>
> V =THETA(2)*EXP(ETA(2)) ; with uncertainty for V
>
> K=CL/V
>
> IF (ICALL.EQ.4.AND.NEWIND.EQ.0) THEN ; do this just once per subproblem
>
> DO WHILE (ETA(3).LT.-.99.OR.ETA(4).LT.-9.9) CALL SIMETA(ETA) ;resample if “bad” values
>
> END DO
>
> TVCL=THETA(1)+ETA(3)
>
> TVV =THETA(2)+ETA(4)
>
> ENDIF
>
> IF (ICALL.EQ.4) THEN ; do this for all simulated values
>
> CL=TVCL*EXP(ETA(1)) ; with uncertainty for CL
>
> V =TVV*EXP(ETA(2)) ; with uncertainty for V
>
> ENDIF
>
> $ERROR
>
> Y =F*EXP(EPS(1))
>
> $INFN
>
> IF (NEWIND.EQ.0) THEN
>
> WRITE (52,*) THETA ; Output estimated THETAs to file 52
>
> WRITE (53,*) TVCL, TVV ; Output simulated THETAs to file 53
>
> ENDIF
>
> $EST MAX=9999 METH=1 INTER
>
> Mats Karlsson, PhD
>
> Professor of Pharmacometrics
>
> Dept of Pharmaceutical Biosciences
>
> Uppsala University
>
> Box 591
>
> 751 24 Uppsala Sweden
>
> phone: +46 18 4714105
>
> fax: +46 18 471 4003
Nick Holford wrote:
> $SIM (20090709) ONLYSIM SUBPROBLEMS=100
> ; estimates of THETA and OMEGA from previous run
> $THETA
> 1 ; POP_CL theta1
> 10 ; POP_V theta2
> $OMEGA
> 0.5 ; PPV_CL eta1
> 0.5 ; PPV_V eta2
> ;variance-covariance matrix of the THETA estimates from previous run
> $OMEGA BLOCK(2)
> 0.2 ; UNC_POP_CL eta3
> 0.1 3 ; UNC_POP_V eta 4
>
> $PK
> ; get CL and V uncertainties
> IF (NEWIND.EQ.0) THEN ; do this just once per subproblem
> UNCCL=THETA(1)+ETA(3)
> UNCV=THETA(2)+ETA(4)
> ENDIF
> CL=UNCCL*EXP(ETA(1)) ; with uncertainty for CL
> V =UNCV*EXP(ETA(2)) ; with uncertainty for V
> ...
--
Nick Holford, Professor Clinical Pharmacology
Dept Pharmacology & Clinical Pharmacology
University of Auckland, 85 Park Rd, Private Bag 92019, Auckland, New Zealand
[email protected] tel:+64(9)923-6730 fax:+64(9)373-7090
mobile: +64 21 46 23 53
http://www.fmhs.auckland.ac.nz/sms/pharmacology/holford