Re: AW: Simulations with/without residual error
Mats,
Thanks for pointing out the error in my original code fragement -- it was, of course, an intentional model misspecification :-)
The code as you write it below won't run for me because of PREDPP errors so I dont know how you get any output at all but I do accept that it does need to have ICALL.EQ.4 in the block that calculates the initial values of UNCCL and UNCV.
IF (ICALL.EQ.4.AND.NEWIND.EQ.0) THEN ; do this just once per subproblem
UNCCL=THETA(1)*EXP(ETA(3))
UNCV=THETA(2)*EXP(ETA(4))
ENDIF
I don't understand why ICALL.EQ.4 is needed because this is a $SIM ONLYSIM problem. I expected the code to run as if it was all in an ICALL.EQ.4 block but it seems that once again NONMEM does the unexpected thing -- no doubt consist with some deeply hidden rule in the documentation.
When ICALL.EQ.4 is included and ETA is resampled to get reasonable values for CL and V then this is the output with the same value of UNCCL and UNCV for all records in each subproblem.
TABLE NO. 1
REP ID UNCCL CL UNCV V
1 1 1.3889 0.50415 11.591 6.1054
1 1 1.3889 0.50415 11.591 6.1054
1 1 1.3889 0.50415 11.591 6.1054
1 2 1.3889 3.1992 11.591 10.021
1 2 1.3889 3.1992 11.591 10.021
1 2 1.3889 3.1992 11.591 10.021
TABLE NO. 1
REP ID UNCCL CL UNCV V
2 1 1.5026 2.6093 11.656 2.4064
2 1 1.5026 2.6093 11.656 2.4064
2 1 1.5026 2.6093 11.656 2.4064
2 2 1.5026 0.66622 11.656 4.6768
2 2 1.5026 0.66622 11.656 4.6768
2 2 1.5026 0.66622 11.656 4.6768
Best wishes,
Nick
Mats Karlsson wrote:
> Nick,
>
> I said the code was incorrect because it actually simulates values dependent
> on the uncertainty for the very first record in each problem only. Here is
> output from your code:
>
> LINE NO. ID TIME DV UNCC UNCV CL V
> DV PRED RES WRES
>
> 1
>
> + 1.00E+00 0.00E+00 0.00E+00 1.36E+00 9.38E+00 2.36E+00
> 3.23E+00 0.00E+00 1.00E+02 0.00E+00 0.00E+00
>
> 2
>
> + 1.00E+00 1.00E+00 1.#RE+00 0.00E+00 0.00E+00 0.00E+00
> 0.00E+00 1.#RE+00 1.#RE+00 1.#RE+00 0.00E+00
>
> 3
>
> + 1.00E+00 2.10E+01 1.#RE+00 0.00E+00 0.00E+00 0.00E+00
> 0.00E+00 1.#RE+00 1.#RE+00 1.#RE+00 0.00E+00
>
> 4
>
> + 2.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00
> 0.00E+00 0.00E+00 1.00E+02 0.00E+00 0.00E+00
>
> And here is that code:
> $PROB
> $DATA data2.csv
> $INPUT ID DV TIME DROP DROP AMT
>
> $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
> $SIGMA .1
> $SUB ADVAN1
> $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
>
> K=CL/V
> $ERROR
> Y =F*EXP(EPS(1))
> $TABLE ID TIME DV UNCCL UNCV CL V
>
> As for the error model this was of course intentional - checking impact of
> model misspecification is one of the reasons we do these kinds of simulation
> studies. :)
>
> Best regards,
> Mats
>
> 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
>
Quoted reply history
> -----Original Message-----
> From: [email protected] [mailto:[email protected]] On
> Behalf Of Nick Holford
> Sent: Sunday, July 26, 2009 10:42 PM
> To: nmusers
> Subject: Re: AW: [NMusers] 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