Hi Mats, Nick and all.
Two questions on the proposed model for simulation using the uncertainty in the
covariance matrix that was discussed last summer.
1. I would like to perform simulation with uncertainty in parameter estimates
followed by estimation with alternative models. Would it be possible to use the
model Mats proposed as the simulation model using the SSE module in PsN?
2. I did run the model proposed by Mats.
When I look at the output files 52 (estimated thetas) and 53 (similated
thetas), the mean of the estimated thetas are only around half of the simulated
values. Also there is no correlation between simulated and estimated thetas.
How can this be? Is there something I am misunderstanding.
I used a datafile with 20 subjects with two observations/subject as explained
below.
The model I used was:
$PROB
$DATA data1.csv IGNORE=@
$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 .001
$SUB ADVAN1
$PK
CL=THETA(1)*EXP(ETA(1)) ; with uncertainty for CL
V =THETA(2)*EXP(ETA(2)) ; with uncertainty for V
TCLE=THETA(1); estimated typical value
TVE=THETA(2); estimated typical value
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
$TABLE ID DV TIME TVCL TVV TCLE TVE CL V ETA1 ETA2 ETA3 ETA4
ONEHEADER NOPRINT FILE=patab37
--------------------------------------------------------------------------
Confidentiality Notice: This message is private and may contain confidential
and proprietary information. If you have received this message in error, please
notify us and remove it from your system and note that you must not copy,
distribute or take any action in reliance on it. Any unauthorized use or
disclosure of the contents of this message is not permitted and may be unlawful.
-----Original Message-----
From: owner-nmus...@globomaxnm.com [mailto:owner-nmus...@globomaxnm.com] On
Behalf Of Nick Holford
Sent: den 27 juli 2009 09:42
To: nmusers
Subject: Re: AW: [NMusers] 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
>
>
> -----Original Message-----
> From: owner-nmus...@globomaxnm.com [mailto:owner-nmus...@globomaxnm.com] 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
n.holf...@auckland.ac.nz 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
fort52
Description: fort52
fort53
Description: fort53
