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

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 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