NONMEM Users Network Archive

Hosted by Cognigen

Re: AW: Simulations with/without residual error

From: Nick Holford <n.holford>
Date: Mon, 27 Jul 2009 19:42:22 +1200

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
>
>
> -----Original Message-----
> From: owner-nmusers
> 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
n.holford
mobile: +64 21 46 23 53
http://www.fmhs.auckland.ac.nz/sms/pharmacology/holford
Received on Mon Jul 27 2009 - 03:42:22 EDT

The NONMEM Users Network is maintained by ICON plc. Requests to subscribe to the network should be sent to: nmusers-request@iconplc.com.

Once subscribed, you may contribute to the discussion by emailing: nmusers@globomaxnm.com.