Re: Unexpected differences in predictions between NM 6.2.0 and NM 7.1.2
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