From:Paul HutsonSubject:[NMusers] Storing simulations Date:Thu, 11 Jul 2002 13:57:19 -0500 Dear NM Users, In the archive (99feb152000.html), Dr. DeJongh asked how to automate MonteCarlo sims and store the output for later processing. Dr. Banken's reply mentioned the $SUPER loop, but it is not apparent that the simulated "individual" data sets are being stored for later fitting as a generated set of population data. Are there other control streams that others have used to generate (and store) sets of "population" data? Thanks in advance. Paul Paul Hutson, Pharm.D. Associate Professor (CHS) UW School of Pharmacy 777 Highland Avenue Madison, WI 53705-2222 Tel: (608) 263-2496 FAX: (608) 265-5421 Pager: (608) 265-7000, #7856 ******* From:Paul Hutson Subject:Re: [NMusers] Storing simulations Date:Thu, 11 Jul 2002 14:50:13 -0500 Before I get a flood of "Read the Manual" regarding my request, a clarification of a poorly written question... In seeking randomness from initial etas and sigmas in a PK model that is used to generate a new set of "population data" using SIMULATION, is it useful or necessary to use different SIMETA values to affect the generated values of ETA(1..i) and EPS(1..j)? Thank you ******* From:Nick Holford Subject:Re: [NMusers] Storing simulations Date:Fri, 12 Jul 2002 08:14:54 +1200 I am not sure what you are trying to do but let me guess. If you want to simulate a data set using NONMEM then estimate the parameters of that data set using NONMEM you have 2 choices: 1. One Stage Stage 1: Add the $SIMULATION (1234567) SUBPROBLEMS=100 to an existing control stream. This will use the dataset as a template, simulate a new DV which replaces the original DV value in a temporary (hidden from the user) dataset, then estimate the parameters using the simulated data. With 100 subproblems this will be repeateed 100 times. The parameter estimates will appear in the usual NONMEM output listing but instead of one problem you will find the results of 100 problems listed one after another. Sounds too easy? 'Fraid so. NONMEM has a bug/feature that means that if one of the subproblems encounters certain errors (which arise commonly when dealing with simulated data and using the same initial estimates for all runs) it will not go onto the next problem but just stops. This unreliability makes it worthless for serious work. 2. Two Stage Stage 1: Add the $SIMULATION (1234567) ONLSIMULATION to an existing control stream. This will simulate a new data set as in the one stage method but you can then save the simulated results (and ID and other covariates) in a TABLE file in the usual way. The generated "population data" will be saved in the TABLE file. Stage 2: Change the dataset in your control stream to the name of the table file generated in Stage 1. (Remove any $SIMULATION record!). This will then run NONMEM in the usual way. the But if you are doing Monte Carlo simulation and need to do this 100 times or more then you need to find a way to automate Stages 1 and 2. With both methods you will need to find an automated way of extracting things such as objective function values and parameter estimates in order to perform a meta-analysis. Wings for NONMEM ( http://wfn.sourceforge.net) will handle this for you if you use the One Stage method. I have also made programs for doing the Two Stage method but they are a bit more clumsy and not currently supported directly in the distributed version of WFN. The basic idea is to run Stage 1 with ONLYSIMULATION and SUBPROBLEMS=100. Then extract, one at time, the 100 simulated data sets saved in the TABLE file output and run each of these extracted data sets with a standard NONMEM control stream. -- Nick Holford, Divn Pharmacology & Clinical Pharmacology University of Auckland, 85 Park Rd, Private Bag 92019, Auckland, New Zealand email:n.holford@auckland.ac.nz tel:+64(9)373-7599x6730 fax:373-7556 http://www.health.auckland.ac.nz/pharmacology/staff/nholford/ ******* From: Nick Holford Subject: Re: [NMusers] Storing simulations Date:Fri, 12 Jul 2002 08:39:24 +1200 I doubt if the Users Guides will be of much use. In the simple case (see my previous email) you do not need to use SIMETA() or SIMEPS() when using $SIMULATION. But typically one needs to be able to truncate the distribution of ETA or EPS values e.g. If you use an additive residual error model then sooner or later you will generate a large enough negative EPS to give you a negative simulated concentration. Same problem with additive models for parameter variability -- you will generate negative clearances. Even if you use an exponential model which cannot generate negative values you may still wish to impose some reasonable constraints on the simulated DV values or simulated parameters. This is when you use SIMETA() and SIMEPS(). e.g. this is a simple pharmacodynamic problem that first of all ensures that the DV is positive (the DOWHILE loop) then makes the DV missing if the simulated value is less than 200. $SIGMA 100 ; SD $PRED E0=E0*EXP(CVE0) EMAX=EMAX*EXP(CVEMAX) EC50=EC50*EXP(CVEC50) EFFECT=E0 + EMAX*THEO/(THEO+EC50) Y = EFFECT + SD IF (ICALL.EQ.4) THEN ; simulation DOWHILE (Y.LT.0) CALL SIMEPS(EPS) Y=EFFECT + SD ENDDO IF (Y.LT.200) THEN MDV=1 ENDIF ENDIF If you want to control the simulated parameters then it might look like the following code which will throw away parameters which are less likely than 1 in a 1000 to have arisen from the distribution. Please note that you must sample all the ETAs at the same time if they are correlated via an OMEGA BLOCK. You will need to be more creative in making the DOWHILE condition fit in the 80 char line limit if you have more than 3 parameters. $THETA (0,150.,) ; E0 $THETA (0,200.,) ; EMAX $THETA (.001,10,) ; EC50 $OMEGA 0.5 ; CVE0 $OMEGA 0.5 ; CVEMAX $OMEGA 0.5 ; CVEC50 $PRED E0=E0*EXP(CVE0) EMAX=EMAX*EXP(CVEMAX) EC50=EC50*EXP(CVEC50) IF (ICALL.EQ.4) THEN ; simulation TRUNC=3.27 ; Z 2tailed alpha=0.01 i.e. include 99.9% TVE0=E0 E0=TVE0*EXP(CVE0) LNMU=LOG(TVE0) DLTA=TRUNC*0.717 ; MUST BE 3.27*SQRT(CVE0)! L0=EXP(LNMU-DLTA) H0=EXP(LNMU+DLTA) TVEMAX=EMAX EM=TVEMAX*EXP(CVEMAX) LNMU=LOG(TVEMAX) DLTA=TRUNC*0.717 ; MUST BE 3.27*SQRT(CVEMAX)! LM=EXP(LNMU-DLTA) HM=EXP(LNMU+DLTA) TVEC50=EC50 C5=TVEC50*EXP(CVEC50) LNMU=LOG(TVEC50) DLTA=TRUNC*0.717 ; MUST BE 3.27*SQRT(CVEC50)! L5=EXP(LNMU-DLTA) H5=EXP(LNMU+DLTA) DOWHILE(E0.LT.L0.OR.E0.GT.H0.OR.EM.LT.LM.OR.EM.GT.HM.OR.C5.LT.L5.OR.C5.GT.H5) CALL SIMETA(ETA) EMAX=TVEMAX*EXP(CVEMAX) E0=TVE0*EXP(CVE0) EC50=TVEC50*EXP(CVEC50) ENDDO ENDIF EFFECT=E0 + EM*THEO/(THEO+C5) EFFECT=E0 + EMAX*THEO/(THEO+EC50) -- Nick Holford, Divn Pharmacology & Clinical Pharmacology University of Auckland, 85 Park Rd, Private Bag 92019, Auckland, New Zealand email:n.holford@auckland.ac.nz tel:+64(9)373-7599x6730 fax:373-7556 http://www.health.auckland.ac.nz/pharmacology/staff/nholford/