Hello NMUsers
I am trying to reanalyse some data using the M3 method (Ahn et al, Likelihood
based approaches to handling data below the quantification limit using NONMEM
VI, J Pharmacokinet Pharmacodyn. 2008). I had previously used M6 and wanted to
try something potentially more robust.
I'm having problems with a "Hessian of posterior density is non positive
definite during search" error occurring as soon as NONMEM is started. It is
suggested that this may be caused by including time zero records in the
analysis (when MDV or EVID = 0) but I've double-checked and I can't find any
incidence of this.
My control stream is pasted below and I can send an extract of data if that
helps?
Any ideas very gratefully received!
Many thanks :)
Ann
$PROBLEM PK 3comp LAG BLQ M3
$INPUT ID DOSE AMT RATE DUR TIME DV EVID ART FLG AGE WGT MDV
$DATA step1&2PKBLQ.csv IGNORE=#
$SUBROUTINES ADVAN9 TOL=9
$MODEL
COMP(CENTRAL, DEFOBS) ;1
COMP(PERIPH1) ;2
COMP(PERIPH2) ;3
$PK
CL=THETA(1)*EXP(ETA(1))
Q2=THETA(2)*EXP(ETA(2))
Q3=THETA(3)*EXP(ETA(3))
V1=THETA(4)*EXP(ETA(4))
V2=THETA(5)*EXP(ETA(5))
V3=THETA(6)*EXP(ETA(6))
K10=CL/V1
K12=Q2/V1
K13=Q3/V1
K21=Q2/V2
K31=Q3/V3
S1=V1
IF (DUR.EQ.10) THEN
ALAG1=THETA(7)*EXP(ETA(7))
ELSE
ALAG1=THETA(8)*EXP(ETA(7))
ENDIF
$DES
DADT(1)=A(2)*K21 + A(3)*K31 - A(1)*(K10+K12+K13)
DADT(2)=A(1)*K12 - A(2)*K21
DADT(3)=A(1)*K13 - A(3)*K31
$ERROR
SIG1=THETA(9)
LOQ = 0.4 ; nm/L
IPRED = F
DUM = (LOQ-IPRED)/(SIG1*IPRED)
CUMD=PHI(DUM)
IF(FLG.EQ.0)THEN ; FLG=0, MDV=0 non BQL values
F_FLAG=0
Y = F*(1+SIG1*ERR(1))
ENDIF
IF(FLG.EQ.1)THEN ; FLG=1, MDV=0 BQL values
F_FLAG=1
Y=CUMD
ENDIF
$THETA (0, 696) ;CL
$THETA (0, 213) ;Q2
$THETA (0, 2600) ;Q3
$THETA (0, 5600) ;V1
$THETA (0, 50500) ;V2
$THETA (0, 38000) ;V3
$THETA (0.2167, 0.548,1) ;ALAG 10MIN
$THETA (0.00833, 0.101, 1) ;ALAG 1MIN
$THETA (0, 0.2) ;SIG1
$OMEGA (0.0121) ; ETA CL
$OMEGA (0.455) ; ETA Q2
$OMEGA (0.118) ; ETA Q3
$OMEGA (0.0919) ; ETA V1
$OMEGA (0.324) ; ETA V2
$OMEGA (0.0664) ; ETA V3
$OMEGA (1.02) ; ETA LAG SHARED
$SIGMA (0.0622)
$ESTIMATION METHOD=1 LAPLACIAN NUMERICAL SLOW PRINT=1 MAX=9999 NOABORT SIG=3
$COVA
$TABLE ID EVID AMT TIME IPRED
NOPRINT FILE=AllRecords.txt
$TABLE ID CL Q2 Q3 V1 V2 V3
ETA1 ETA2 ETA3 ETA4 ETA5 ETA6 ETA7
FIRSTONLY NOPRINT NOAPPEND FILE=FirstRecords.txt
_______________________________________________________________________
Ann Rigby-Jones PhD MRSC
Research Fellow in Pharmacokinetics & Pharmacodynamics
Peninsula College of Medicine & Dentistry, UK
_______________________________________________________________________
Hessian of posterior density is non positive definite during search
Ann,
This error message usually means try NONMEM thinks your OMEGA initial estimates are weird. A leading suspect would be the OMEGA estimates for the two lag time parameters. Try fixing both of them to zero and see if that gets your model going. If that doesn't work then you could try fixing all OMEGA values to zero then adding them in one at a time.
You don't need to use $DES for your linear model. $SUB ADVAN11 is a closed form 3 cpt model with bolus/zero order input (or ADVAN12 with first order absorption as well).
An efficiency tip -- there is no point calling the time consuming PHI() function until you really need it and also no need for a FLG data item. You can simplify your $ERROR code like this:
$ERROR
SD=F*THETA(9) ; any residual error model you like assuming $SIGMA 1 FIX
LLOQ = 0.4 ; nm/L
IF (DV.GE.LLOQ) THEN ; non BQL values
F_FLAG=0
Y = F + SD*ERR(1))
ELSE ; BQL values
F_FLAG=1
Y = PHI((LLOQ-F)/SD)
ENDIF
Best wishes,
Nick
Quoted reply history
On 29/10/2010 9:50 p.m., Ann Rigby-Jones wrote:
> Hello NMUsers
>
> I am trying to reanalyse some data using the M3 method (Ahn et al, Likelihood based approaches to handling data below the quantification limit using NONMEM VI, J Pharmacokinet Pharmacodyn. 2008). I had previously used M6 and wanted to try something potentially more robust.
>
> I'm having problems with a "Hessian of posterior density is non positive definite during search" error occurring as soon as NONMEM is started. It is suggested that this may be caused by including time zero records in the analysis (when MDV or EVID = 0) but I've double-checked and I can't find any incidence of this.
>
> My control stream is pasted below and I can sendan extract of data if that helps?
>
> Any ideas very gratefully received!
>
> Many thanks J
>
> Ann
>
> $PROBLEM PK 3comp LAG BLQ M3
> $INPUT ID DOSE AMT RATE DUR TIME DV EVID ART FLG AGE WGT MDV
> $DATA step1&2PKBLQ.csv IGNORE=#
> $SUBROUTINES ADVAN9 TOL=9
> $MODEL
> COMP(CENTRAL, DEFOBS) ;1
> COMP(PERIPH1) ;2
> COMP(PERIPH2) ;3
>
> $PK
>
> CL=THETA(1)*EXP(ETA(1))
> Q2=THETA(2)*EXP(ETA(2))
> Q3=THETA(3)*EXP(ETA(3))
> V1=THETA(4)*EXP(ETA(4))
> V2=THETA(5)*EXP(ETA(5))
> V3=THETA(6)*EXP(ETA(6))
>
> K10=CL/V1
> K12=Q2/V1
> K13=Q3/V1
> K21=Q2/V2
> K31=Q3/V3
>
> S1=V1
>
> IF (DUR.EQ.10) THEN
> ALAG1=THETA(7)*EXP(ETA(7))
> ELSE
> ALAG1=THETA(8)*EXP(ETA(7))
> ENDIF
>
> $DES
> DADT(1)=A(2)*K21 + A(3)*K31 - A(1)*(K10+K12+K13)
> DADT(2)=A(1)*K12 - A(2)*K21
> DADT(3)=A(1)*K13 - A(3)*K31
>
> $ERROR
>
> SIG1=THETA(9)
> LOQ = 0.4 ; nm/L
> IPRED = F
> DUM = (LOQ-IPRED)/(SIG1*IPRED)
> CUMD=PHI(DUM)
> IF(FLG.EQ.0)THEN ; FLG=0, MDV=0 non BQL values
> F_FLAG=0
> Y = F*(1+SIG1*ERR(1))
> ENDIF
> IF(FLG.EQ.1)THEN ; FLG=1, MDV=0 BQL values
> F_FLAG=1
> Y=CUMD
> ENDIF
>
> $THETA (0, 696) ;CL
> $THETA (0, 213) ;Q2
> $THETA (0, 2600) ;Q3
> $THETA (0, 5600) ;V1
> $THETA (0, 50500) ;V2
> $THETA (0, 38000) ;V3
> $THETA (0.2167, 0.548,1) ;ALAG 10MIN
> $THETA (0.00833, 0.101, 1) ;ALAG 1MIN
> $THETA (0, 0.2) ;SIG1
>
> $OMEGA (0.0121) ; ETA CL
> $OMEGA (0.455) ; ETA Q2
> $OMEGA (0.118) ; ETA Q3
> $OMEGA (0.0919) ; ETA V1
> $OMEGA (0.324) ; ETA V2
> $OMEGA (0.0664) ; ETA V3
> $OMEGA (1.02) ; ETA LAG SHARED
>
> $SIGMA (0.0622)
>
> $ESTIMATION METHOD=1 LAPLACIAN NUMERICAL SLOW PRINT=1 MAX=9999 NOABORT SIG=3
>
> $COVA
> $TABLE ID EVID AMT TIME IPRED
> NOPRINT FILE=AllRecords.txt
> $TABLE ID CL Q2 Q3 V1 V2 V3
> ETA1 ETA2 ETA3 ETA4 ETA5 ETA6 ETA7
> FIRSTONLY NOPRINT NOAPPEND FILE=FirstRecords.txt
>
> _______________________________________________________________________
>
> *Ann Rigby-Jones PhD MRSC*
> Research Fellow in Pharmacokinetics & Pharmacodynamics
>
> Peninsula College of Medicine & Dentistry, UK
>
> _______________________________________________________________________
--
Nick Holford, Professor Clinical Pharmacology
Dept Pharmacology& Clinical Pharmacology
University of Auckland,85 Park Rd,Private Bag 92019,Auckland,New Zealand
tel:+64(9)923-6730 fax:+64(9)373-7090 mobile:+64(21)46 23 53
email: [email protected]
http://www.fmhs.auckland.ac.nz/sms/pharmacology/holford
Ann,
I had run into similar problems. Here're some suggestions for you to
consider:
1. I had included only one post-Cmax BLQ value. It worked for my case.
2. One consultant proposed some boundaries on the parameters to constrain
the search. He also suggest me to stabilize my model first. It's work in
progress.
3. ALAG may have caused problem, especially with differential equation,
along with ETA.
Best regards,
Shelley
========================================
Xiao Hu (Shelley), Ph.D.
Scientist,
Development Pharmacokinetics & Disposition
Biogen Idec, Inc.
14 Cambridge Center
Cambridge, MA 02142
Ann Rigby-Jones <[email protected]>
Sent by: [email protected]
29-Oct-2010 10:03 AM
Message Size: 32.5 KB
To
"'[email protected]'" <[email protected]>
cc
Subject
RE: [NMusers] Hessian of posterior density is non positive definite during
search
Dear Nick
Many thanks for your reply J The reason I am using $DES is because as my
next step will evaluate adding transit compartments to the model. I?m
currently trying to sort the late kinetics before I tackle the early
disposition. I?ve already investigated transit models with this data set
having just deleted BQL data but the WRES vs time plots show bias at the
late time points (FOCE).
I?ve been trying out your suggestions (and others suggested by Jeroen
Elassaiss, many thanks Jeroen);
If I fix all ETAs to zero, I then get problems with:
ERROR: NONMEM terminated with the error message:
SUM OF "SQUARED" WEIGHTED INDIVIDUAL RESIDUALS IS INFINITE
If I allow an ETA on CL only (initial estimate 0.01) I get
0PRED EXIT CODE = 1
0INDIVIDUAL NO. 1 ID= 1.00000000000000E+00 (WITHIN-INDIVIDUAL)
DATA REC NO. 2
THETA=
6.96E+02 5.13E+02 2.60E+02 5.60E+03 2.05E+04 3.80E+04 5.48E-01
1.01E-01 2.00E-01
OCCURS DURING SEARCH FOR ETA AT A NONZERO VALUE OF ETA
ERROR IN LSODI1: CODE 21
ERROR OCURRED WHILE ATTEMPTING TO OBTAIN INITIAL VALUES FOR ATOL
Elapsed estimation time in seconds: 0.05
This time point is my first occurrence of BQL (FLG=1), not sure if this is
relevant. I don?t think the initial estimates of THETA are too far out,
as above, this data set has been modelling before using FOCE and omitted
BQL data points.
I?m wondering whether log transforming the data may help, but unsure of
the implementation with M3.
Many thanks for your input, much appreciated J
All best wishes
Ann
Quoted reply history
From: [email protected] [mailto:[email protected]]
On Behalf Of Nick Holford
Sent: Friday, 29 October, 2010 11:50
To: '[email protected]'
Subject: Re: [NMusers] Hessian of posterior density is non positive
definite during search
Ann,
This error message usually means try NONMEM thinks your OMEGA initial
estimates are weird. A leading suspect would be the OMEGA estimates for
the two lag time parameters. Try fixing both of them to zero and see if
that gets your model going. If that doesn't work then you could try fixing
all OMEGA values to zero then adding them in one at a time.
You don't need to use $DES for your linear model. $SUB ADVAN11 is a closed
form 3 cpt model with bolus/zero order input (or ADVAN12 with first order
absorption as well).
An efficiency tip -- there is no point calling the time consuming PHI()
function until you really need it and also no need for a FLG data item.
You can simplify your $ERROR code like this:
$ERROR
SD=F*THETA(9) ; any residual error model you like assuming $SIGMA 1 FIX
LLOQ = 0.4 ; nm/L
IF (DV.GE.LLOQ) THEN ; non BQL values
F_FLAG=0
Y = F + SD*ERR(1))
ELSE ; BQL values
F_FLAG=1
Y = PHI((LLOQ-F)/SD)
ENDIF
Best wishes,
Nick
On 29/10/2010 9:50 p.m., Ann Rigby-Jones wrote:
Hello NMUsers
I am trying to reanalyse some data using the M3 method (Ahn et al,
Likelihood based approaches to handling data below the quantification
limit using NONMEM VI, J Pharmacokinet Pharmacodyn. 2008). I had
previously used M6 and wanted to try something potentially more robust.
I'm having problems with a "Hessian of posterior density is non positive
definite during search" error occurring as soon as NONMEM is started. It
is suggested that this may be caused by including time zero records in the
analysis (when MDV or EVID = 0) but I've double-checked and I can't find
any incidence of this.
My control stream is pasted below and I can send an extract of data if
that helps?
Any ideas very gratefully received!
Many thanks J
Ann
$PROBLEM PK 3comp LAG BLQ M3
$INPUT ID DOSE AMT RATE DUR TIME DV EVID ART FLG AGE WGT MDV
$DATA step1&2PKBLQ.csv IGNORE=#
$SUBROUTINES ADVAN9 TOL=9
$MODEL
COMP(CENTRAL, DEFOBS) ;1
COMP(PERIPH1) ;2
COMP(PERIPH2) ;3
$PK
CL=THETA(1)*EXP(ETA(1))
Q2=THETA(2)*EXP(ETA(2))
Q3=THETA(3)*EXP(ETA(3))
V1=THETA(4)*EXP(ETA(4))
V2=THETA(5)*EXP(ETA(5))
V3=THETA(6)*EXP(ETA(6))
K10=CL/V1
K12=Q2/V1
K13=Q3/V1
K21=Q2/V2
K31=Q3/V3
S1=V1
IF (DUR.EQ.10) THEN
ALAG1=THETA(7)*EXP(ETA(7))
ELSE
ALAG1=THETA(8)*EXP(ETA(7))
ENDIF
$DES
DADT(1)=A(2)*K21 + A(3)*K31 - A(1)*(K10+K12+K13)
DADT(2)=A(1)*K12 - A(2)*K21
DADT(3)=A(1)*K13 - A(3)*K31
$ERROR
SIG1=THETA(9)
LOQ = 0.4 ; nm/L
IPRED = F
DUM = (LOQ-IPRED)/(SIG1*IPRED)
CUMD=PHI(DUM)
IF(FLG.EQ.0)THEN ; FLG=0, MDV=0 non BQL values
F_FLAG=0
Y = F*(1+SIG1*ERR(1))
ENDIF
IF(FLG.EQ.1)THEN ; FLG=1, MDV=0 BQL values
F_FLAG=1
Y=CUMD
ENDIF
$THETA (0, 696) ;CL
$THETA (0, 213) ;Q2
$THETA (0, 2600) ;Q3
$THETA (0, 5600) ;V1
$THETA (0, 50500) ;V2
$THETA (0, 38000) ;V3
$THETA (0.2167, 0.548,1) ;ALAG 10MIN
$THETA (0.00833, 0.101, 1) ;ALAG 1MIN
$THETA (0, 0.2) ;SIG1
$OMEGA (0.0121) ; ETA CL
$OMEGA (0.455) ; ETA Q2
$OMEGA (0.118) ; ETA Q3
$OMEGA (0.0919) ; ETA V1
$OMEGA (0.324) ; ETA V2
$OMEGA (0.0664) ; ETA V3
$OMEGA (1.02) ; ETA LAG SHARED
$SIGMA (0.0622)
$ESTIMATION METHOD=1 LAPLACIAN NUMERICAL SLOW PRINT=1 MAX=9999 NOABORT
SIG=3
$COVA
$TABLE ID EVID AMT TIME IPRED
NOPRINT FILE=AllRecords.txt
$TABLE ID CL Q2 Q3 V1 V2 V3
ETA1 ETA2 ETA3 ETA4 ETA5 ETA6 ETA7
FIRSTONLY NOPRINT NOAPPEND FILE=FirstRecords.txt
_______________________________________________________________________
Ann Rigby-Jones PhD MRSC
Research Fellow in Pharmacokinetics & Pharmacodynamics
Peninsula College of Medicine & Dentistry, UK
_______________________________________________________________________
--
Nick Holford, Professor Clinical Pharmacology
Dept Pharmacology & Clinical Pharmacology
University of Auckland,85 Park Rd,Private Bag 92019,Auckland,New Zealand
tel:+64(9)923-6730 fax:+64(9)373-7090 mobile:+64(21)46 23 53
email: [email protected]
http://www.fmhs.auckland.ac.nz/sms/pharmacology/holford
This message and any attachments are solely for the intended recipient. If
you are not the intended recipient, disclosure, copying, use or
distribution of the information included in this message is prohibited ---
Please immediately and permanently delete.
Dear Matt & All
I'm still battling I'm afraid...
Can't seem to get past the SUM OF "SQUARED" WEIGHTED INDIVIDUAL RESIDUALS IS
INFINITE problem.
I've now log-transformed my data, and deleted all but the first occurrence of a
BQL data point from each individual. I'm confident that the initial estimates
of the structural parameters aren't too far out as I have already analysed the
data set having deleted all BQL points.
Any other ideas very gratefully received :)
Many thanks
Ann
$PROBLEM Step1&2 PK 3comp LAG BLQ M3
$INPUT ID DOSE AMT RATE DUR TIME LNDV DV EVID ART FLG AGE WGT MDV
$DATA step1&2PKBLQnoABSLNv2.csv IGNORE=#
$SUBROUTINES ADVAN9 TOL=9
$MODEL
COMP(CENTRAL, DEFOBS, DEFDOSE) ;1
COMP(PERIPH1) ;2
COMP(PERIPH2) ;3
$PK
CL=THETA(1)*EXP(ETA(1))
Q2=THETA(2)*EXP(ETA(2))
Q3=THETA(3)*EXP(ETA(3))
V1=THETA(4)*EXP(ETA(4))
V2=THETA(5)*EXP(ETA(5))
V3=THETA(6)*EXP(ETA(6))
K10=CL/V1
K12=Q2/V1
K13=Q3/V1
K21=Q2/V2
K31=Q3/V3
S1=V1
IF (DUR.EQ.10) THEN
ALAG1=THETA(7)*EXP(ETA(7))
ELSE
ALAG1=THETA(8)*EXP(ETA(7))
ENDIF
$DES
DADT(1)=A(2)*K21 + A(3)*K31 - A(1)*(K10+K12+K13)
DADT(2)=A(1)*K12 - A(2)*K21
DADT(3)=A(1)*K13 - A(3)*K31
;DADT(1)=A(4)*KTR + A(2)*K21 + A(3)*K31 - A(1)*(K10+K12+K13)
;DADT(2)=A(1)*K12 - A(2)*K21
;DADT(3)=A(1)*K13 - A(3)*K31
;DADT(4)=-A(4)*KTR
;DADT(5)=A(4)*KTR - A(5)*KTR
;DADT(6)=A(5)*KTR - A(6)*KTR
;DADT(7)=A(6)*KTR - A(7)*KTR
;DADT(8)=A(7)*KTR - A(8)*KTR
;DADT(9)=A(8)*KTR - A(9)*KTR
$ERROR (ONLY OBSERVATIONS)
SIG1=THETA(9)
;LOQ = 0.4 ; nm/L
;IPRED = F
LOQ = LOG(0.4)
IPRED = LOG(F)
DUM = (LOQ-IPRED)/(SIG1*IPRED)
CUMD=PHI(DUM)
IF(FLG.EQ.0)THEN ; FLG=0, MDV=0 non BQL values
F_FLAG=0
;Y = F*(1+SIG1*ERR(1))
Y=LOG(F) + SIG1*ERR(1)
ENDIF
IF(FLG.EQ.1)THEN ; FLG=1, MDV=0 BQL values
F_FLAG=1
Y=CUMD
ENDIF
$THETA (0, 696) ;CL
$THETA (0, 513) ;Q2
$THETA (0, 260) ;Q3
$THETA (0, 5600) ;V1
$THETA (0, 20500) ;V2
$THETA (0, 38000) ;V3
$THETA (0.2167, 0.548,1) ;ALAG 10MIN
$THETA (0.00833, 0.101, 1) ;ALAG 1MIN
$THETA (0, 0.5) ;SIG1
$OMEGA (0.01) ; ETA CL
$OMEGA (0.01) ; ETA Q2
$OMEGA (0 FIX) ; ETA Q3
$OMEGA (0.01) ; ETA V1
$OMEGA (0 FIX) ; ETA V2
$OMEGA (0 FIX) ; ETA V3
$OMEGA (0 FIX) ; ETA LAG SHARED
$SIGMA (1 FIXED)
$ESTIMATION METHOD=1 LAPLACIAN NUMERICAL SLOW PRINT=1 MAX=9999 NOABORT SIG=3
MSFO=msfo.outputfile
$COVA
$TABLE ID EVID AMT TIME IPRED
NOPRINT FILE=AllRecords.txt
$TABLE ID CL Q2 Q3 V1 V2 V3
ETA1 ETA2 ETA3 ETA4 ETA5 ETA6 ETA7
FIRSTONLY NOPRINT NOAPPEND FILE=FirstRecords.txt
Ann,
Several ideas:
1. Fix ALAGs to zero, will it help?
2. Try to shield IPRED from being zero
IPRED=LOG(0.0001)
IF(F.GT.0.0001) IPRED=LOG(F)
Y=IPRED + SIG1*ERR(1)
(this could be a bad style, but will allow you to test the code)
3. Try ADVAN6, 8, 13
4. Remove BQL part, fit it as a 3-comp linear model
5. Check whether you have any FLG values NE 0 or 1 in the dataset
6. Try (to guard against problem with FLG values):
F_FLAG=0
Y=LOG(F) + SIG1*ERR(1)
IF(FLG.EQ.1)THEN
F_FLAG=1
Y=CUMD
ENDIF
---------------------
At least, this will help you to understand where the problem is
Leonid
--------------------------------------
Leonid Gibiansky, Ph.D.
President, QuantPharm LLC
web: www.quantpharm.com
e-mail: LGibiansky at quantpharm.com
tel: (301) 767 5566
Quoted reply history
On 11/1/2010 11:42 AM, Ann Rigby-Jones wrote:
> Dear Matt & All
>
> I’m still battling I’m afraid…
>
> Can’t seem to get past theSUM OF "SQUARED" WEIGHTED INDIVIDUAL RESIDUALS
> IS INFINITE problem.
>
> I’ve now log-transformed my data, and deleted all but the first
> occurrence of a BQL data point from each individual. I’m confident that
> the initial estimates of the structural parameters aren’t too far out as
> I have already analysed the data set having deleted all BQL points.
>
> Any other ideas very gratefully received J
>
> Many thanks
>
> Ann
>
> $PROBLEM Step1&2 PK 3comp LAG BLQ M3
>
> $INPUT ID DOSE AMT RATE DUR TIME LNDV DV EVID ART FLG AGE WGT MDV
>
> $DATA step1&2PKBLQnoABSLNv2.csv IGNORE=#
>
> $SUBROUTINES ADVAN9 TOL=9
>
> $MODEL
>
> COMP(CENTRAL, DEFOBS, DEFDOSE) ;1
>
> COMP(PERIPH1) ;2
>
> COMP(PERIPH2) ;3
>
> $PK
>
> CL=THETA(1)*EXP(ETA(1))
>
> Q2=THETA(2)*EXP(ETA(2))
>
> Q3=THETA(3)*EXP(ETA(3))
>
> V1=THETA(4)*EXP(ETA(4))
>
> V2=THETA(5)*EXP(ETA(5))
>
> V3=THETA(6)*EXP(ETA(6))
>
> K10=CL/V1
>
> K12=Q2/V1
>
> K13=Q3/V1
>
> K21=Q2/V2
>
> K31=Q3/V3
>
> S1=V1
>
> IF (DUR.EQ.10) THEN
>
> ALAG1=THETA(7)*EXP(ETA(7))
>
> ELSE
>
> ALAG1=THETA(8)*EXP(ETA(7))
>
> ENDIF
>
> $DES
>
> DADT(1)=A(2)*K21 + A(3)*K31 - A(1)*(K10+K12+K13)
>
> DADT(2)=A(1)*K12 - A(2)*K21
>
> DADT(3)=A(1)*K13 - A(3)*K31
>
> ;DADT(1)=A(4)*KTR + A(2)*K21 + A(3)*K31 - A(1)*(K10+K12+K13)
>
> ;DADT(2)=A(1)*K12 - A(2)*K21
>
> ;DADT(3)=A(1)*K13 - A(3)*K31
>
> ;DADT(4)=-A(4)*KTR
>
> ;DADT(5)=A(4)*KTR - A(5)*KTR
>
> ;DADT(6)=A(5)*KTR - A(6)*KTR
>
> ;DADT(7)=A(6)*KTR - A(7)*KTR
>
> ;DADT(8)=A(7)*KTR - A(8)*KTR
>
> ;DADT(9)=A(8)*KTR - A(9)*KTR
>
> $ERROR (ONLY OBSERVATIONS)
>
> SIG1=THETA(9)
>
> ;LOQ = 0.4 ; nm/L
>
> ;IPRED = F
>
> LOQ = LOG(0.4)
>
> IPRED = LOG(F)
>
> DUM = (LOQ-IPRED)/(SIG1*IPRED)
>
> CUMD=PHI(DUM)
>
> IF(FLG.EQ.0)THEN ; FLG=0, MDV=0 non BQL values
>
> F_FLAG=0
>
> ;Y = F*(1+SIG1*ERR(1))
>
> Y=LOG(F) + SIG1*ERR(1)
>
> ENDIF
>
> IF(FLG.EQ.1)THEN ; FLG=1, MDV=0 BQL values
>
> F_FLAG=1
>
> Y=CUMD
>
> ENDIF
>
> $THETA (0, 696) ;CL
>
> $THETA (0, 513) ;Q2
>
> $THETA (0, 260) ;Q3
>
> $THETA (0, 5600) ;V1
>
> $THETA (0, 20500) ;V2
>
> $THETA (0, 38000) ;V3
>
> $THETA (0.2167, 0.548,1) ;ALAG 10MIN
>
> $THETA (0.00833, 0.101, 1) ;ALAG 1MIN
>
> $THETA (0, 0.5) ;SIG1
>
> $OMEGA (0.01) ; ETA CL
>
> $OMEGA (0.01) ; ETA Q2
>
> $OMEGA (0 FIX) ; ETA Q3
>
> $OMEGA (0.01) ; ETA V1
>
> $OMEGA (0 FIX) ; ETA V2
>
> $OMEGA (0 FIX) ; ETA V3
>
> $OMEGA (0 FIX) ; ETA LAG SHARED
>
> $SIGMA (1 FIXED)
>
> $ESTIMATION METHOD=1 LAPLACIAN NUMERICAL SLOW PRINT=1 MAX=9999 NOABORT
> SIG=3
>
> MSFO=msfo.outputfile
>
> $COVA
>
> $TABLE ID EVID AMT TIME IPRED
>
> NOPRINT FILE=AllRecords.txt
>
> $TABLE ID CL Q2 Q3 V1 V2 V3
>
> ETA1 ETA2 ETA3 ETA4 ETA5 ETA6 ETA7
>
> FIRSTONLY NOPRINT NOAPPEND FILE=FirstRecords.txt
Ann,
Because you already log-tranformed your data, you may want to change
SIG1*IPRED to SIG1 for DUM.
I'm guessing the combination of ALAG and SIG1*IPRED is causing problem.
DUM = (LOQ-IPRED)/(SIG1)
Best regards,
Shelley
========================================
Xiao Hu (Shelley), Ph.D.
Scientist,
Development Pharmacokinetics & Disposition
Biogen Idec, Inc.
14 Cambridge Center
Cambridge, MA 02142
Ann Rigby-Jones <[email protected]>
Sent by: [email protected]
01-Nov-2010 11:42 AM
Message Size: 32.2 KB
To
"'[email protected]'" <[email protected]>
cc
Subject
RE: [NMusers] Hessian of posterior density is non positive definite during
search
Dear Matt & All
I?m still battling I?m afraid?
Can?t seem to get past the SUM OF "SQUARED" WEIGHTED INDIVIDUAL RESIDUALS
IS INFINITE problem.
I?ve now log-transformed my data, and deleted all but the first occurrence
of a BQL data point from each individual. I?m confident that the initial
estimates of the structural parameters aren?t too far out as I have
already analysed the data set having deleted all BQL points.
Any other ideas very gratefully received J
Many thanks
Ann
$PROBLEM Step1&2 PK 3comp LAG BLQ M3
$INPUT ID DOSE AMT RATE DUR TIME LNDV DV EVID ART FLG AGE WGT MDV
$DATA step1&2PKBLQnoABSLNv2.csv IGNORE=#
$SUBROUTINES ADVAN9 TOL=9
$MODEL
COMP(CENTRAL, DEFOBS, DEFDOSE) ;1
COMP(PERIPH1) ;2
COMP(PERIPH2) ;3
$PK
CL=THETA(1)*EXP(ETA(1))
Q2=THETA(2)*EXP(ETA(2))
Q3=THETA(3)*EXP(ETA(3))
V1=THETA(4)*EXP(ETA(4))
V2=THETA(5)*EXP(ETA(5))
V3=THETA(6)*EXP(ETA(6))
K10=CL/V1
K12=Q2/V1
K13=Q3/V1
K21=Q2/V2
K31=Q3/V3
S1=V1
IF (DUR.EQ.10) THEN
ALAG1=THETA(7)*EXP(ETA(7))
ELSE
ALAG1=THETA(8)*EXP(ETA(7))
ENDIF
$DES
DADT(1)=A(2)*K21 + A(3)*K31 - A(1)*(K10+K12+K13)
DADT(2)=A(1)*K12 - A(2)*K21
DADT(3)=A(1)*K13 - A(3)*K31
;DADT(1)=A(4)*KTR + A(2)*K21 + A(3)*K31 - A(1)*(K10+K12+K13)
;DADT(2)=A(1)*K12 - A(2)*K21
;DADT(3)=A(1)*K13 - A(3)*K31
;DADT(4)=-A(4)*KTR
;DADT(5)=A(4)*KTR - A(5)*KTR
;DADT(6)=A(5)*KTR - A(6)*KTR
;DADT(7)=A(6)*KTR - A(7)*KTR
;DADT(8)=A(7)*KTR - A(8)*KTR
;DADT(9)=A(8)*KTR - A(9)*KTR
$ERROR (ONLY OBSERVATIONS)
SIG1=THETA(9)
;LOQ = 0.4 ; nm/L
;IPRED = F
LOQ = LOG(0.4)
IPRED = LOG(F)
DUM = (LOQ-IPRED)/(SIG1*IPRED)
CUMD=PHI(DUM)
IF(FLG.EQ.0)THEN ; FLG=0, MDV=0 non BQL values
F_FLAG=0
;Y = F*(1+SIG1*ERR(1))
Y=LOG(F) + SIG1*ERR(1)
ENDIF
IF(FLG.EQ.1)THEN ; FLG=1, MDV=0 BQL values
F_FLAG=1
Y=CUMD
ENDIF
$THETA (0, 696) ;CL
$THETA (0, 513) ;Q2
$THETA (0, 260) ;Q3
$THETA (0, 5600) ;V1
$THETA (0, 20500) ;V2
$THETA (0, 38000) ;V3
$THETA (0.2167, 0.548,1) ;ALAG 10MIN
$THETA (0.00833, 0.101, 1) ;ALAG 1MIN
$THETA (0, 0.5) ;SIG1
$OMEGA (0.01) ; ETA CL
$OMEGA (0.01) ; ETA Q2
$OMEGA (0 FIX) ; ETA Q3
$OMEGA (0.01) ; ETA V1
$OMEGA (0 FIX) ; ETA V2
$OMEGA (0 FIX) ; ETA V3
$OMEGA (0 FIX) ; ETA LAG SHARED
$SIGMA (1 FIXED)
$ESTIMATION METHOD=1 LAPLACIAN NUMERICAL SLOW PRINT=1 MAX=9999 NOABORT
SIG=3
MSFO=msfo.outputfile
$COVA
$TABLE ID EVID AMT TIME IPRED
NOPRINT FILE=AllRecords.txt
$TABLE ID CL Q2 Q3 V1 V2 V3
ETA1 ETA2 ETA3 ETA4 ETA5 ETA6 ETA7
FIRSTONLY NOPRINT NOAPPEND FILE=FirstRecords.txt
Hello Everyone
Huge thanks to Peter, Jeroen, Leonid, Jonathan, Nick, Matt, Marc, Andreas,
Kyun-Seop, Shelley and all, for taking the time to look at my problem and for
offering me such detailed and helpful advice, much appreciated. Apologies for
the delay in replying, teaching duties have kept me from my modelling over the
last few days.
I've not finished with the analysis yet but certainly the main problem seemed
to lie with my use of ALAG as many of you suspected. Once I removed these
parameters from the control stream, the model minimized successfully (the late
BLQ points seemed to bias the model fit but that's a separate issue!). My next
step is to try to implement the method with transit compartments instead of
ALAG...
Very many thanks - have a great weekend!
Best wishes
Ann
_______________________________________________________________________
Ann Rigby-Jones PhD MRSC
Research Fellow in Pharmacokinetics & Pharmacodynamics
Peninsula College of Medicine & Dentistry, UK
_______________________________________________________________________