From:Paul Hutson
Subject:[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
Storing simulations
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/