RE: output SE of population PK parameters in Nonmem
Dear All
I think you have to write your own INFN subroutine.
You may write something like this:
The config file SIZES must be on the same place as this code and must be the
same as the SIZES used by your NONMEM installation
c Begin of own subroutine
SUBROUTINE INFN (ICALL,THETA,DATREC,INDXS,NEWIND)
INCLUDE 'SIZES'
INTEGER PE,PT,PCOV,IERE,IERC,I,MODE
PARAMETER (PE=LVR)
PARAMETER (PT=LTH)
PARAMETER (PCOV=LPAR3)
DIMENSION THETA(*),DATREC(*),INDXS(*)
DOUBLE PRECISION THETA
COMMON /ROCM6/ THETAF(PT),OMEGAF(PE,PE),SIGMAF(PE,PE)
COMMON /ROCM7/ SETH(PT),SEOM(PE,PE),SESIG(PE,PE)
COMMON /ROCM8/ OBJECT
COMMON /ROCM9/ IERE,IERC
COMMON /CM12/ COV(PCOV)
COMMON /CM2/ NETA
DOUBLE PRECISION OBJECT, OBJECT1 DELTA
DOUBLE PRECISION THETAF,OMEGAF,SIGMAF
DOUBLE PRECISION INFTH(PT),INFOM(PE),INFSIG(PE)
DOUBLE PRECISION SUPTH(PT),SUPOM(PE),SUPSIG(PE)
DOUBLE PRECISION ETA(PE)
INTEGER NTH,NETA,NEPS
CHARACTER*16 FILECOV
1060 FORMAT(80F23.10)
IF (ICALL.EQ.0) THEN
OPEN(50,FILE=filecov,status='REPLACE')
IPROB=1
ENDIF
IF (ICALL.EQ.3) THEN
DO 20 I=1,PCOV
c PRINT *,'INFO',COV(I)
IF (COV(I).EQ.0) THEN
CLOSE(50)
GOTO 30
ENDIF
c PRINT *,'INFO',I,COV(I)
WRITE (50,1060) COV(I)
20 CONTINUE
30 CONTINUE
ENDIF
RETURN
END
C End of own subroutine
In the contol file
Write the following:
$SUBROUTINES ADVAN4 TRANS4 INFN=myown.f
The result will be a vector, which gives you the diagonal half matrix of the
covariance matrix of estimate. (loaded in R it may look like the following lines
1 9.497490e-02
2 1.000000e+10
3 1.000000e+10
4 1.825657e-01
5 1.000000e+10
6 1.332238e+00
....
228 -2.389190e-05
229 1.000000e+10
230 6.468500e-06
231 1.393760e-05
The number of lines in this case 231 lines, shows you that you have a 21x21
Matrix.
This is in the example the way to fill the results into the matrix
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[,15] [,16] [,17] [,18] [,19] [,20] [,21]
[1,] 1 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
[2,] 2 3 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
[3,] 4 5 6 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
[4,] 7 8 9 10 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
[5,] 11 12 13 14 15 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
[6,] 16 17 18 19 20 21 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
[7,] 22 23 24 25 26 27 28 0 0 0 0 0 0
0 0 0 0 0 0 0 0
[8,] 29 30 31 32 33 34 35 36 0 0 0 0 0
0 0 0 0 0 0 0 0
[9,] 37 38 39 40 41 42 43 44 45 0 0 0 0
0 0 0 0 0 0 0 0
[10,] 46 47 48 49 50 51 52 53 54 55 0 0 0
0 0 0 0 0 0 0 0
[11,] 56 57 58 59 60 61 62 63 64 65 66 0 0
0 0 0 0 0 0 0 0
[12,] 67 68 69 70 71 72 73 74 75 76 77 78 0
0 0 0 0 0 0 0 0
[13,] 79 80 81 82 83 84 85 86 87 88 89 90 91
0 0 0 0 0 0 0 0
[14,] 92 93 94 95 96 97 98 99 100 101 102 103 104
105 0 0 0 0 0 0 0
[15,] 106 107 108 109 110 111 112 113 114 115 116 117 118
119 120 0 0 0 0 0 0
[16,] 121 122 123 124 125 126 127 128 129 130 131 132 133
134 135 136 0 0 0 0 0
[17,] 137 138 139 140 141 142 143 144 145 146 147 148 149
150 151 152 153 0 0 0 0
[18,] 154 155 156 157 158 159 160 161 162 163 164 165 166
167 168 169 170 171 0 0 0
[19,] 172 173 174 175 176 177 178 179 180 181 182 183 184
185 186 187 188 189 190 0 0
[20,] 191 192 193 194 195 196 197 198 199 200 201 202 203
204 205 206 207 208 209 210 0
[21,] 211 212 213 214 215 216 217 218 219 220 221 222 223
224 225 226 227 228 229 230 231
Some coding in R and you have the result of the covariance matrix of estimate
in R
Than you can calculate in R with cov2cor the CORRELATION MATRIX OF ESTIMATE
Bevor you should delete all rows and columns equal to 1.000000e+10 (INF in
NONMEM) and as result you should have a covariance matrix of estimates in R
with your calculated estimates.
I hope this helps.
Mit freundlichen Grüßen / Best regards / Cordialement
Dipl. Ing. Heiner Speth
Sanofi-Aventis Deutschland GmbH
GMPK
Industriepark Hoechst
Bldg. H831, Room 364
D-65926 Frankfurt am Main
t: +49 69 305 2148
f: +49 69 305 942 052
w: www.sanofi-aventis.de
**************************************************************************************************************************************************************
Sanofi-Aventis Deutschland GmbH · Sitz der Gesellschaft: Frankfurt am Main ·
Handelsregister: Frankfurt am Main, Abt. B Nr. 40661
Vorsitzender des Aufsichtsrats: Hanspeter Spek - Geschäftsführer: Dr.
Heinz-Werner Meier (Vorsitzender),
Dr. Matthias Braun, Hervé Gisserot, Prof. Dr. Dr. Werner Kramer, Dr. Klaus
Menken, Dr. Martin Siewert
**************************************************************************************************************************************************************
Quoted reply history
-----Original Message-----
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Rik Schoemaker
Sent: Wednesday, February 06, 2008 2:31 PM
To: [email protected]
Subject: RE: [NMusers] output SE of population PK parameters in Nonmem
Dear all,
It's a neat trick and works very well, but can anyone tell me the keyword (like
for instance SETHET below) for exporting the correlation matrix of the
estimates and for the non-parametric estimates: 'expected value of ETA' and
'Covariance matrix of ETA'....
Thanks in advance,
Rik
-----Original Message-----
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of GIRARD PASCAL
Sent: 01 February 2008 18:08
To: Jin, Bo; [email protected]
Subject: RE : [NMusers] output SE of population PK parameters in Nonmem
Hi Jin,
See in NONMEM help:
___________________________________________________________________
| |
| FINALIZATION EXAMPLE ($PRED AND $INFN) |
|_________________________________________________________________|
This example contains abbreviated code which can be inserted in a $PRED
or $INFN block. This code outputs final parameter estimates, standard
errors, minimum value of the objective function, and condi- tional estimates
of etas to various user files. The return codes from Estimation and
Covariance steps (zero for normal termination) are also output.
IF (ICALL.EQ.3) THEN
DO WHILE(DATA)
IF (NEWIND.LE.1) WRITE (50,*) ETA
ENDDO
WRITE (51,*) OBJECT
WRITE (52,*) THETA
WRITE (53,*) SETHET
WRITE (54,*) OMEGA(BLOCK)
WRITE (55,*) SEOMEG(BLOCK)
WRITE (56,*) SIGMA(BLOCK)
WRITE (57,*) SESIGM(BLOCK)
WRITE (58,*) IERE,IERC
ENDIF
Alternatively, you take advantage of the fact that all those estimmates are
stored in Vectors and matrices of various COMMON:
COMMON /ROCM6/ THETAF(40),OMEGAF(30,30),SIGMAF(30,30)
COMMON /ROCM7/ SETH(40),SEOM(30,30),SESIG(30,30)
COMMON /ROCM8/ OBJECT
COMMON /ROCM9/ IERE,IERC
COMMON /CM12/ COV(2850)
and write your own INFN subroutine (see INFN help).
Best regards,
Pascal Girard, PhD
EA 3738, CTO
Fac Medecine Lyon-Sud, BP12
69921 OULLINS Cedex, France
[EMAIL PROTECTED]
Tel +33 (0)4 26 23 59 54 / Fax +33 (0)4 26 23 59 76
Master Recherche Lyon 1 Santé et Populations, Spécialité PhIT
http://master-sante-pop.univ-lyon1.fr/
>-----Message d'origine-----
>De : [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
>De la part de Jin, Bo Envoyé : vendredi 1 février 2008 17:48 À :
>[email protected] Objet : [NMusers] output SE of population PK
>parameters in Nonmem
>
>Dear all: I have a question. Is there any way to output the SE or SD
>of a population PK parameter estimate in Nonmem (say typical value of
>Cl) into some files which can be read by SAS or S-Plus etc.? (e.g.
>some excel or TAB files?)
>
>thanks,
>
>- Bo
>
>_____________________________________
>
>BO JIN
>Clinical Pharmacology Statistics
>Merck Research Labs
>Phone: 267-305-7876
>
>
>
>-----------------------------------------------------------------------
>----
>---
>Notice: This e-mail message, together with any attachments, contains
>information of Merck & Co., Inc. (One Merck Drive, Whitehouse Station,
>New Jersey, USA 08889), and/or its affiliates (which may be known
>outside the United States as Merck Frosst, Merck Sharp & Dohme or MSD
>and in Japan, as Banyu - direct contact information for affiliates is
>available at http://www.merck.com/contact/contacts.html) that may be
>confidential, proprietary copyrighted and/or legally privileged. It is
>intended solely for the use of the individual or entity named on this
>message. If you are not the intended recipient, and have received this
>message in error, please notify us immediately by reply e-mail and then
>delete it from your system.
>
>-----------------------------------------------------------------------
>----
>---