Dear R-users,
I have recently observed a very puzzling difference of behavior between NM 6.2.0 and NM 7.1.2, while trying to reproduce the simulation examples described in Martin Bergstrand's paper 'Handling data below the limit of quantification in mixed effect models' (AAPS Journal 2009, 11-2). My model corresponds to Martins' model B (2-cmt linear model with single IV dosing) and implement my interpretation of the residual variability described in his equation 3 (see code below). When I run the code in NONMEM 7.1.2,, all the PRED values are reported as NaN and RES values either as NaN or 0. NM6 returns non-NaN values for PRED and RES, but there does not seem to be any differences in simulated DV.
Did anybody experienced the same differences? In my case, could they be explained by an improper implementation of the RV model?
As a side note, I would also be interested to know if a reference in the literature would describe the properties of this particular RUV model.
Thank you
Sebastien Bihorel
--------
$PROBLEM base-2cmt-sim
$DATA basedata.csv IGNORE=@
$INPUT ID TIME AMT RATE CMT EVID DV MDV STDY
$THETA
(0.,5.) ;1- clearance
(0.,20.) ;2- central compartment volume
(0.,5.) ;3- distribution clearance
(0.,100.) ;4- peripheral compartment volume
(0.,0.1) ;5- 'additive' log RV
(0.,0.1) ;6- 'second' log RV term
$OMEGA
0.3 ;1- IIV in clearance
0.3 ;2- IIV in central compartment volume
0.3 ;3- IIV in distribution clearance
0.3 ;4- IIV in peripheral compartment volume
$SIGMA
1 FIX ;1- 'additive' log RV
1 FIX ;2- 'second' log RV term
$SUBROUTINES ADVAN3 TRANS4
$PK
; Model parameter assignment
TVCL=THETA(1)
TVV1=THETA(2)
TVQ =THETA(3)
TVV2=THETA(4)
ECL=EXP(ETA(1))
EV1=EXP(ETA(2))
EQ =EXP(ETA(3))
EV2=EXP(ETA(4))
; PREDPP required variables
F1=1.0 ; bioavailability in central compartment
CL=TVCL*ECL ; elimination clearance
V1=TVV1*EV1 ; volume of the central compartment
Q =TVQ*EQ ; inter-compartment distribution clearance
V2=TVV2*EV2 ; volume of the peripheral compartment
S1=V1
$ERROR
;set up Dose flag
DFLG=0
IF(AMT.GT.0)DFLG=1
IPRED=LOG(F+DFLG)
W1=THETA(5)*EPS(1)
W2=THETA(6)*EPS(2)
W=SQRT(W1*W1+(W2*(1-F/(0.5+F)))**2)
Y=IPRED+W
$SIM (123456) ONLYSIM
$TABLE ID TIME AMT RATE CMT EVID DV MDV STDY W NOPRINT ONEHEADER
FILE=base-2cmt-sim.tbl
basedata.csv
Description:
MS-Excel spreadsheet
Unexpected differences in predictions between NM 6.2.0 and NM 7.1.2
Dear Sebastien,
I can't answer your primary question but I want to make a correction on the
way that I specified the residual error model when using it (see updated
code below). Your implementation will in my opinion cause bias since it will
make all residuals positive (-x*-x = x^2).
The reason that I used this model for simulations was to avoid the problems
with the estimation error model, i.e. [W = SQRT(W1**2+(W2/F)**2]. These
problems have been mentioned before at NMusers
( http://www.mail-archive.com/[email protected]/msg02445.html).
The error model that you are referring wasn't used for estimation since it
depends on tree parameters (W1,W2 and WH) and these are in most cases not
all identifiable. I didn't get this error model out of the literature
(haven't seen it described), it was simply my own solution to the problems I
was experiencing.
;; --- Altered error model code ---------
IPRED=LOG(F+DFLG)
W1=THETA(5)
W2=THETA(6)
WH=THETA(7) ; Equal to 0.5 in examples out of the publication)
W=SQRT(W1*W1+(W2*(1-F/(WH+F)))**2)
Y=IPRED+W*EPS(1)
$SIMGA 1 FIX
;; --------------------------------------
Best regards,
Martin Bergstrand, MSc, PhD student
-----------------------------------------------
Pharmacometrics Research Group,
Department of Pharmaceutical Biosciences,
Uppsala University
-----------------------------------------------
[email protected]
-----------------------------------------------
Work: +46 18 471 4639
Mobile: +46 709 994 396
Quoted reply history
-----Original Message-----
From: [email protected] [mailto:[email protected]] On
Behalf Of Sebastien Bihorel
Sent: Thursday, May 13, 2010 6:06 PM
To: nmusers; Martin Bergstrand
Subject: [NMusers] Unexpected differences in predictions between NM 6.2.0
and NM 7.1.2
Dear R-users,
I have recently observed a very puzzling difference of behavior between NM
6.2.0 and NM 7.1.2, while trying to reproduce the simulation examples
described in Martin Bergstrand's paper 'Handling data below the limit of
quantification in mixed effect models' (AAPS Journal 2009, 11-2). My model
corresponds to Martins' model B (2-cmt linear model with single IV
dosing) and implement my interpretation of the residual variability
described in his equation 3 (see code below). When I run the code in NONMEM
7.1.2,, all the PRED values are reported as NaN and RES values either as NaN
or 0. NM6 returns non-NaN values for PRED and RES, but there does not seem
to be any differences in simulated DV.
Did anybody experienced the same differences? In my case, could they be
explained by an improper implementation of the RV model?
As a side note, I would also be interested to know if a reference in the
literature would describe the properties of this particular RUV model.
Thank you
Sebastien Bihorel
--------
$PROBLEM base-2cmt-sim
$DATA basedata.csv IGNORE=@
$INPUT ID TIME AMT RATE CMT EVID DV MDV STDY
$THETA
(0.,5.) ;1- clearance
(0.,20.) ;2- central compartment volume
(0.,5.) ;3- distribution clearance
(0.,100.) ;4- peripheral compartment volume
(0.,0.1) ;5- 'additive' log RV
(0.,0.1) ;6- 'second' log RV term
$OMEGA
0.3 ;1- IIV in clearance
0.3 ;2- IIV in central compartment volume
0.3 ;3- IIV in distribution clearance
0.3 ;4- IIV in peripheral compartment volume
$SIGMA
1 FIX ;1- 'additive' log RV
1 FIX ;2- 'second' log RV term
$SUBROUTINES ADVAN3 TRANS4
$PK
; Model parameter assignment
TVCL=THETA(1)
TVV1=THETA(2)
TVQ =THETA(3)
TVV2=THETA(4)
ECL=EXP(ETA(1))
EV1=EXP(ETA(2))
EQ =EXP(ETA(3))
EV2=EXP(ETA(4))
; PREDPP required variables
F1=1.0 ; bioavailability in central compartment
CL=TVCL*ECL ; elimination clearance
V1=TVV1*EV1 ; volume of the central compartment
Q =TVQ*EQ ; inter-compartment distribution clearance
V2=TVV2*EV2 ; volume of the peripheral compartment
S1=V1
$ERROR
;set up Dose flag
DFLG=0
IF(AMT.GT.0)DFLG=1
IPRED=LOG(F+DFLG)
W1=THETA(5)*EPS(1)
W2=THETA(6)*EPS(2)
W=SQRT(W1*W1+(W2*(1-F/(0.5+F)))**2)
Y=IPRED+W
$SIM (123456) ONLYSIM
$TABLE ID TIME AMT RATE CMT EVID DV MDV STDY W NOPRINT ONEHEADER
FILE=base-2cmt-sim.tbl
Thanks Martin,
Running the model in NM7 using your parameterization returns non-NaN values for DV, PRED and RES. I agree with you that my parametrization would have force the residues to be positive but I'm still puzzled by NM7 results: W is within 0.001 and 1, which should not create numerical problems... As a side note, I am a bit concerned with the use of a single sigma parameter for simulation purposes. I believe this will induce a correlation between the additive and the pseudo-proportional parts of the RV. Shouldn't we use the following parameterization for simulation?
DFLG=0
IF(AMT.GT.0)DFLG=1
IPRED=LOG(F+DFLG)
W1=THETA(5)*EPS(1)
W2=THETA(6)*EPS(2)*(1-F/(THETA(7)+F))
Y=IPRED+W1+W2
This should theoretically give the same variance for Y as your parameterization and enable the two components of the RV model to be independent during the simulation.
Sebastien
Martin Bergstrand wrote:
> Dear Sebastien,
>
> I can't answer your primary question but I want to make a correction on the
> way that I specified the residual error model when using it (see updated
> code below). Your implementation will in my opinion cause bias since it will
>
> make all residuals positive (-x*-x = x^2).
>
> The reason that I used this model for simulations was to avoid the problems
> with the estimation error model, i.e. [W = SQRT(W1**2+(W2/F)**2]. These
> problems have been mentioned before at NMusers
>
> ( http://www.mail-archive.com/ [email protected] /msg02445.html ).
>
> The error model that you are referring wasn't used for estimation since it
> depends on tree parameters (W1,W2 and WH) and these are in most cases not
> all identifiable. I didn't get this error model out of the literature
> (haven't seen it described), it was simply my own solution to the problems I
> was experiencing.
>
> ;; --- Altered error model code ---------
> IPRED=LOG(F+DFLG)
> W1=THETA(5)
> W2=THETA(6)
>
> WH=THETA(7) ; Equal to 0.5 in examples out of the publication)
>
> W=SQRT(W1*W1+(W2*(1-F/(WH+F)))**2)
>
> Y=IPRED+W*EPS(1)
>
> $SIMGA 1 FIX
> ;; --------------------------------------
>
> Best regards,
>
> Martin Bergstrand, MSc, PhD student
> -----------------------------------------------
> Pharmacometrics Research Group,
> Department of Pharmaceutical Biosciences,
> Uppsala University
> -----------------------------------------------
> [email protected]
> -----------------------------------------------
> Work: +46 18 471 4639
> Mobile: +46 709 994 396
>
Quoted reply history
> -----Original Message-----
> From: [email protected] [mailto:[email protected]] On
> Behalf Of Sebastien Bihorel
> Sent: Thursday, May 13, 2010 6:06 PM
> To: nmusers; Martin Bergstrand
> Subject: [NMusers] Unexpected differences in predictions between NM 6.2.0
> and NM 7.1.2
>
> Dear R-users,
>
> I have recently observed a very puzzling difference of behavior between NM
> 6.2.0 and NM 7.1.2, while trying to reproduce the simulation examples
> described in Martin Bergstrand's paper 'Handling data below the limit of
> quantification in mixed effect models' (AAPS Journal 2009, 11-2). My model
> corresponds to Martins' model B (2-cmt linear model with single IV
> dosing) and implement my interpretation of the residual variability
> described in his equation 3 (see code below). When I run the code in NONMEM
> 7.1.2,, all the PRED values are reported as NaN and RES values either as NaN
> or 0. NM6 returns non-NaN values for PRED and RES, but there does not seem
> to be any differences in simulated DV.
>
> Did anybody experienced the same differences? In my case, could they be
> explained by an improper implementation of the RV model?
>
> As a side note, I would also be interested to know if a reference in the
> literature would describe the properties of this particular RUV model.
>
> Thank you
>
> Sebastien Bihorel
>
> --------
>
> $PROBLEM base-2cmt-sim
>
> $DATA basedata.csv IGNORE=@
>
> $INPUT ID TIME AMT RATE CMT EVID DV MDV STDY
>
> $THETA
> (0.,5.) ;1- clearance
> (0.,20.) ;2- central compartment volume
> (0.,5.) ;3- distribution clearance
> (0.,100.) ;4- peripheral compartment volume
> (0.,0.1) ;5- 'additive' log RV
> (0.,0.1) ;6- 'second' log RV term
>
> $OMEGA
> 0.3 ;1- IIV in clearance
> 0.3 ;2- IIV in central compartment volume
> 0.3 ;3- IIV in distribution clearance
> 0.3 ;4- IIV in peripheral compartment volume
>
> $SIGMA
> 1 FIX ;1- 'additive' log RV
> 1 FIX ;2- 'second' log RV term
>
> $SUBROUTINES ADVAN3 TRANS4
>
> $PK
>
> ; Model parameter assignment
>
> TVCL=THETA(1)
> TVV1=THETA(2)
> TVQ =THETA(3)
> TVV2=THETA(4)
>
> ECL=EXP(ETA(1))
> EV1=EXP(ETA(2))
> EQ =EXP(ETA(3))
> EV2=EXP(ETA(4))
>
> ; PREDPP required variables
>
> F1=1.0 ; bioavailability in central compartment
> CL=TVCL*ECL ; elimination clearance
> V1=TVV1*EV1 ; volume of the central compartment
> Q =TVQ*EQ ; inter-compartment distribution clearance
> V2=TVV2*EV2 ; volume of the peripheral compartment
>
> S1=V1
>
> $ERROR
>
> ;set up Dose flag
> DFLG=0
> IF(AMT.GT.0)DFLG=1
>
> IPRED=LOG(F+DFLG)
> W1=THETA(5)*EPS(1)
> W2=THETA(6)*EPS(2)
> W=SQRT(W1*W1+(W2*(1-F/(0.5+F)))**2)
>
> Y=IPRED+W
>
> $SIM (123456) ONLYSIM
>
> $TABLE ID TIME AMT RATE CMT EVID DV MDV STDY W NOPRINT ONEHEADER
> FILE=base-2cmt-sim.tbl