Examples
Introduction to NONMEM 7.6.0
Ten worked examples demonstrating NONMEM 7 features: ADVAN compartment models, covariate modelling, population mixtures, receptor-mediated disposition, inter-occasion variability, MCMC Bayesian sampling, SAEM annealing, and categorical data analysis.
Two compartment Model, Using ADVAN3, TRANS4
;Model Desc: Two compartment Model, Using ADVAN3, TRANS4
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# Example 1 (from samp5l)
$INPUT C SET ID JID TIME DV=CONC AMT=DOSE RATE EVID MDV CMT CLX
V1X QX V2X SDIX SDSX
$DATA example1.csv IGNORE=C
$SUBROUTINES ADVAN3 TRANS4
$PK
; The THETAs are MU modeled.
; Best that there is a linear relationship between THETAs and Mus
; The linear MU modeling of THETAS allows them to be efficiently
; Gibbs sampled.
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
CL=DEXP(MU_1+ETA(1))
V1=DEXP(MU_2+ETA(2))
Q=DEXP(MU_3+ETA(3))
V2=DEXP(MU_4+ETA(4))
S1=V1
$ERROR
Y = F + F*EPS(1)
; Initial values of THETA
$THETA
(0.001, 2.0) ;[LN(CL)]
(0.001, 2.0) ;[LN(V1)]
(0.001, 2.0) ;[LN(Q)]
(0.001, 2.0) ;[LN(V2)]
;INITIAL values of OMEGA
$OMEGA BLOCK(4)
0.15 ;[P]
0.01 ;[F]
0.15 ;[P]
0.01 ;[F]
0.01 ;[F]
0.15 ;[P]
0.01 ;[F]
0.01 ;[F]
0.01 ;[F]
0.15 ;[P]
;Initial value of SIGMA
$SIGMA
(0.6 ) ;[P]
;Prior information is important for MCMC Bayesian analysis,
;not necessary for maximization methods
;Note the syntax used for defining priors that is available
;as of NONMEM 7.3
$PRIOR NWPRI
; Prior information of THETAS
$THETAP (2.0 FIX)X4
; Variance to prior information of THETAS.
; Because variances are very large, this means that the prior
; information to the THETAS is highly uninformative.
$THETAPV BLOCK(4) FIX VALUES(10000,0.0)
; Prior information to the OMEGAS.
$OMEGAP BLOCK(4) FIX VALUES(0.2,0.0)
; Degrees of freedom to prior OMEGA matrix.
; Because degrees of freedom is very low, equal to the
; the dimension of the prior OMEGA, this means that the
; prior information to the OMEGAS is highly uninformative
$OMEGAPD (4 FIX)
; Prior information to the SIGMAS
$SIGMAP 0.06 FIX
; Degrees of freedom to prior SIGMA matrix.
; Because degrees of freedom is very low, equal to the
; the dimension of the prior SIGMA, this means that the
; prior information to the SIGMA is highly uninformative
$SIGMAPD (1 FIX)
; The first analysis is iterative two-stage,
; maximum of 500 iterations (NITER), iteration results
; are printed every 5 iterations, gradient precision (SIGL) is 4.
; Termination is tested on all of
; the population parameters (CTYPE=3),
; and for less then 2 significant digits change (NSIG).
; Prior information is not necessary for ITS, so NOPRIOR=1.
; The intermediate and final results of the ITS method will be
; recoded in row/column format in example1.ext
$EST METHOD=ITS MAPITER=0 INTERACTION FILE=example1.ext NITER=500
PRINT=5 NOABORT SIGL=4 CTYPE=3 CITER=10
CALPHA=0.05 NOPRIOR=1 NSIG=2
; The results of ITS are used as the initial values for the
; SAEM method. A maximum of 3000 ; stochastic iterations (NBURN)
; is requested, but may end early if statistical test determines
; that variations in all parameters is stationary
; (note that any settings from the previous $EST
; carries over to the next $EST statement, within a $PROB).
; The SAEM is a Monte Carlo process,
; so setting the SEED assures repeatability of results.
; Each iteration obtains only 2 Monte Carlo samples ISAMPLE),
; so they are very fast.
; But many iterations are needed, so PRINT only
; every 100th iteration.
; After the stochastic phase, 500 accumulation iterations will be
; Performed (NITER), to obtain good parameters estimates with
; little stochastic noise.
; As a new FILE has not been given, the SAEM results will append to
; example1.ext.
$EST METHOD=SAEM INTERACTION NBURN=3000 NITER=500 PRINT=100
SEED=1556678 ISAMPLE=2
; After the SAEM method, obtain good estimates of the marginal
; density (objective function),
; along with good estimates of the standard errors.
; This is best done with importance sampling ; (IMP),
; performing the expectation step only (EONLY=1), so that
; final population parameters remain at the final SAEM result.
; Five iterations (NITER) should allow the importance sampling
; proposal density to become stationary.
; This is observed by the objective function settling
; to a particular value (with some stochastic noise).
; By using 3000 Monte Carlo samples
; (ISAMPLE), this assures a precise assessment of standard errors.
$EST METHOD=IMP INTERACTION EONLY=1 NITER=5 ISAMPLE=3000 PRINT=1
SIGL=8 NOPRIOR=1
; The Bayesian analysis is performed.
; While 10000 burn-in iterations are requested as a maximum,
; because the termination test is on (CTYPE<>0, set at the
; first $EST statement), and because the initial parameters are at
; the SAEM result, which is the maximum likelihood position,
; the analysis should settle down to a stationary distribution in
; several hundred iterations.
; Prior information is also used to facilitate Bayesian analysis.
; The individual Bayesian iteration results are important,
; and may be need for post-processing analysis.
; So specify a separate FILE for the Bayesian analysis.
$EST METHOD=BAYES INTERACTION FILE=example1.txt NBURN=10000
NITER=10000 PRINT=100 NOPRIOR=0
; Just for old-times sake, let's see what the traditional
; FOCE method will give us.
; And, remember to introduce a new FILE, so its results won't
; append to our Bayesian FILE.
; Appending to example1.ext with the EM methods is fine.
$EST METHOD=COND INTERACTION MAXEVAL=9999 NSIG=3 SIGL=10
PRINT=5 NOABORT NOPRIOR=1
FILE=example1.ext
; Time for the standard error results.
; You may request a more precise gradient precision (SIGL)
; that differed from that used during estimation.
$COV MATRIX=R PRINT=E UNCONDITIONAL SIGL=12
; Print out results in tables. Include some of the new weighted
; residual types
$TABLE ID TIME PRED RES WRES CPRED CWRES EPRED ERES EWRES NOAPPEND
ONEHEADER FILE=example1.TAB NOPRINT
$TABLE ID CL V1 Q V2 FIRSTONLY NOAPPEND NOPRINT FILE=example1.PAR
$TABLE ID ETA1 ETA2 ETA3 ETA4 FIRSTONLY NOAPPEND
NOPRINT FILE=example1.ETA2 Compartment model with Clearance and central volume modeled with covariates age and gender
;Model Desc: Two Compartment model with Clearance and
; central volume modeled with covariates age and gender
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# example2 (from sampc)
$INPUT C SET ID JID TIME DV=CONC AMT=DOSE RATE EVID MDV CMT GNDR AGE
$DATA example2.csv IGNORE=C
$SUBROUTINES ADVAN3 TRANS4
$PK
; LCLM=log transformed clearance, male
LCLM=THETA(1)
;LCLF=log transformed clearance, female.
LCLF=THETA(2)
; CLAM=CL age slope, male
CLAM=THETA(3)
; CLAF=CL age slope, female
CLAF=THETA(4)
; LV1M=log transformed V1, male
LV1M=THETA(5)
; LV1F=log transformed V1, female
LV1F=THETA(6)
; V1AM=V1 age slope, male
V1AM=THETA(7)
; V1AF=V1 age slope, female
V1AF=THETA(8)
; LAGE=log transformed age
LAGE=DLOG(AGE)
;Mean of ETA1, the inter-subject deviation of Clearance,
; is ultimately modeled as linear function of THETA(1) to THETA(4).
; Relating THETAs to Mus by linear functions is not essential for
; ITS, IMP, or IMPMAP methods, but is very helpful for MCMC methods
; such as SAEM and BAYES.
MU_1=(1.0-GNDR)*(LCLM+LAGE*CLAM) + GNDR*(LCLF+LAGE*CLAF)
; Mean of ETA2, the inter-subject deviation of V1,
; is ultimately modeled as linear function of THETA(5) to THETA(8)
MU_2=(1.0-GNDR)*(LV1M+LAGE*V1AM) + GNDR*(LV1F+LAGE*V1AF)
MU_3=THETA(9)
MU_4=THETA(10)
CL=DEXP(MU_1+ETA(1))
V1=DEXP(MU_2+ETA(2))
Q=DEXP(MU_3+ETA(3))
V2=DEXP(MU_4+ETA(4))
S1=V1
$ERROR
CALLFL=0
; Option to model the residual error coefficient in THETA(11),
; rather than in SIGMA.
SDSL=THETA(11)
W=F*SDSL
Y = F + W*EPS(1)
IPRED=F
IWRES=(DV-F)/W
;Initial THETAs
$THETA
( 0.7 ) ;[LCLM]
( 0.7 ) ;[LCLF]
( 2 ) ;[CLAM]
( 2.0);[CLAF]
( 0.7 ) ;[LV1M]
( 0.7 ) ;[LV1F]
( 2.0 ) ;[V1AM]
( 2.0 ) ;[V1AF]
( 0.7 ) ;[MU_3]
( 0.7 );[MU_4]
( 0.3 ) ;[SDSL]
;Initial OMEGAs
$OMEGA BLOCK(4)
0.5 ;[p]
0.001 ;[f]
0.5 ;[p]
0.001 ;[f]
0.001 ;[f]
0.5 ;[p]
0.001 ;[f]
0.001 ;[f]
0.001 ;[f]
0.5 ;[p]
; SIGMA is 1.0 fixed, serves as unscaled variance for EPS(1).
; THETA(11) takes up the residual error scaling.
$SIGMA
(1.0 FIXED)
;Prior information is important for MCMC Bayesian analysis,
; not necessary for maximization methods
; In this example, only the OMEGAs have a prior distribution,
; the THETAS do not.
; For Bayesian methods, it is most important for at least the
; OMEGAs to have a prior, even an uninformative one,
; to stabilize the analysis. Only if the number of subjects
; exceeds the OMEGA dimension number by at least 100,
; then you may get away without priors on OMEGA for BAYES analysis.
$PRIOR NWPRI
; Prior OMEGA matrix
$OMEGAP BLOCK(4) FIX VALUES(0.01,0.0)
; Degrees of freedom to OMEGA prior matrix:
$OMEGAPD 4 FIX
; The first analysis is iterative two-stage.
; Note that the GRD specification is THETA(11) is a
; Sigma-like parameter. This will allow NONMEM to make
; efficient gradient evaluations for THETA(11), which is useful
; for later IMP,IMPMAP, and SAEM methods, but has no impact on
; ITS and BAYES methods.
$EST METHOD=ITS INTERACTION FILE=example2.ext NITER=1000 NSIG=2
PRINT=5 NOABORT SIGL=8 NOPRIOR=1 CTYPE=3 GRD=TS(11)
; Results of ITS serve as initial parameters for the IMP method.
$EST METHOD=IMP INTERACTION EONLY=0 MAPITER=0 NITER=100 ISAMPLE=300
PRINT=1 SIGL=8
; The results of IMP are used as the initial values for the SAEM method.
$EST METHOD=SAEM NBURN=3000 NITER=2000 PRINT=10 ISAMPLE=2
CTYPE=3 CITER=10 CALPHA=0.05
; After the SAEM method, obtain good estimates of the marginal density
; (objective function),
; along with good estimates of the standard errors.
$EST METHOD=IMP INTERACTION EONLY=1 NITER=5 ISAMPLE=3000
PRINT=1 SIGL=8 SEED=123334
CTYPE=3 CITER=10 CALPHA=0.05
; The Bayesian analysis is performed.
$EST METHOD=BAYES INTERACTION FILE=example2.TXT NBURN=10000
NITER=3000 PRINT=100 NOPRIOR=0
CTYPE=3 CITER=10 CALPHA=0.05
; Just for old-times sake, lets see what the traditional
; FOCE method will give us.
; And, remember to introduce a new FILE, so its results wont
; append to our Bayesian FILE.
$EST METHOD=COND INTERACTION MAXEVAL=9999 FILE=example2.ext NSIG=2
SIGL=14 PRINT=5 NOABORT NOPRIOR=1
$COV MATRIX=R UNCONDITIONALPopulation Mixture Problem in 1 Compartment model, with Volume and rate constant parameters
;Model Desc: Population Mixture Problem in 1 Compartment model,
; with Volume and rate constant parameters and their inter-subject
; variances modeled from two sub-populations
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# example3 (from ad1tr1m2s)
$INPUT C SET ID JID TIME CONC=DV DOSE=AMT RATE EVID MDV CMT VC1 K101
VC2 K102 SIGZ PROB
$DATA example3.csv IGNORE=C
$SUBROUTINES ADVAN1 TRANS1
; The mixture model uses THETA(5) as the mixture proportion parameter,
; defining the proportion of subjects in sub-population 1 (P(1),
; and in sub-population 2 (P(2)
$MIX
P(1)=THETA(5)
P(2)=1.0-THETA(5)
NSPOP=2
$PK
; The MUs should always be unconditionally defined, that is,
; they should never be defined in IF/THEN blocks
; THETA(1) models the Volume of sub-population 1
MU_1=THETA(1)
; THETA(2) models the clearance of sub-population 1
MU_2=THETA(2)
; THETA(3) models the Volume of sub-population 2
MU_3=THETA(3)
; THETA(4) models the clearance of sub-population 2
MU_4=THETA(4)
VCM=DEXP(MU_1+ETA(1))
K10M=DEXP(MU_2+ETA(2))
VCF=DEXP(MU_3+ETA(3))
K10F=DEXP(MU_4+ETA(4))
Q=1
IF(MIXNUM.EQ.2) Q=0
V=Q*VCM+(1.0-Q)*VCF
K=Q*K10M+(1.0-Q)*K10F
S1=V
$ERROR
Y = F + F*EPS(1)
; Initial THETAs
$THETA
(-1000.0 4.3 1000.0) ;[MU_1]
(-1000.0 -2.9 1000.0) ;[MU_2]
(-1000.0 4.3 1000.0) ;[MU_3]
(-1000.0 -0.67 1000.0) ;[MU_4]
(0.0001 0.667 0.9999) ;[P(1)]
;Initial OMEGA block 1, for sub-population 1
$OMEGA BLOCK(2)
.04 ;[p]
.01 ; [f]
.027; [p]
;Initial OMEGA block 2, for sub-population 2
$OMEGA BLOCK(2)
.05; [p]
.01; [f]
.06; [p]
$SIGMA
0.01 ;[p]
; Prior information setup for OMEGAS only
$PRIOR NWPRI
; Prior OMEGA block 1. Note that because the OMEGA is separated
; into blocks, so their priors should have the same block design.
$OMEGAP BLOCK(2)
0.05 FIX
0.0 0.05
; Prior OMEGA block 2
$OMEGAP BLOCK(2)
0.05 FIX
0.0 0.05
; Degrees of Freedom defined for Priors.
; One for each OMEGA block defining each sub-population
$OMEGAPD (2 FIX) (2 FIX)
$EST METHOD=ITS INTERACTION NITER=20 PRINT=1 NOABORT SIGL=8
FILE=example3.ext CTYPE=3 CITER=10
CALPHA=0.05 NOPRIOR=1
$EST NBURN=500 NITER=500 METHOD=SAEM INTERACTION PRINT=10 SIGL=6
ISAMPLE=2
$EST METHOD=IMP INTERACTION NITER=5 ISAMPLE=1000 PRINT=1 NOABORT
SIGL=6 EONLY=1 MAPITER=0
$EST METHOD=BAYES INTERACTION NBURN=2000 NITER=1000 PRINT=10
FILE=example3.txt SIGL=8 NOPRIOR=0
$EST MAXEVAL=9999 NSIG=3 SIGL=12 PRINT=1 FILE=example3.ext
METHOD=CONDITIONAL INTERACTION NOABORT
NOPRIOR=1
$COV MATRIX=R UNCONDITIONALPopulation Mixture Problem in 1 Compartment model, with rate constant parameter
Population Mixture Problem in 1 Compartment model, with rate constant parameter and its inter-subject variances modeled as coming from two sub-populations.
;Model Desc: Population Mixture Problem in 1 Compartment model,
; with rate constant parameter and its inter-subject variances
; modeled as coming from two sub-populations
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# example4 (from ad1tr1m2t)
$INPUT C SET ID JID TIME CONC=DV DOSE=AMT RATE EVID MDV CMT VC1
K101 VC2 K102 SIGZ PROB
$DATA example4.csv IGNORE=C
$SUBROUTINES ADVAN1 TRANS1
$MIX
P(1)=THETA(4)
P(2)=1.0-THETA(4)
NSPOP=2
$PK
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
K10M=DEXP(MU_1+ETA(1))
V=DEXP(MU_2+ETA(2))
K10F=DEXP(MU_3+ETA(3))
Q=1
IF(MIXNUM.EQ.2) Q=0
K=Q*K10M+(1.0-Q)*K10F
S1=V
$ERROR
Y = F + F*EPS(1)
$THETA ;[P(1)]
-2.9 ; K10M
4.3 ; V
-0.67 ; K10F
(0.0001 0.667 0.9999)
$OMEGA BLOCK(3)
.027 ;[p]
0.01 ;[f]
.04 ;[p]
0.0 ;[f] ; BAND MATRIX: K10M X K10F COVARIANCE IS 0
0.01 ;[f]
0.06 ;[p]
$SIGMA
0.01 ;[p]
; Prior information setup for OMEGAS only
$PRIOR NWPRI
; Prior OMEGA
$OMEGAP BLOCK(3)
0.05 FIX
0.0 0.05
0.0 0.0 0.05
; Degrees of Freedom defined for Priors.
$OMEGAPD (3 FIX)
$EST METHOD=ITS INTERACTION NITER=30 PRINT=5 NOABORT SIGL=6
FILE=example4.ext NOPRIOR=1 CTYPE=3 CITER=10 CALPHA=0.05
$EST METHOD=IMP INTERACTION NITER=20 ISAMPLE=300 PRINT=1
NOABORT SIGL=6 NOPRIOR=1
$EST NBURN=500 NITER=500 METHOD=SAEM INTERACTION PRINT=10 SIGL=6
ISAMPLE=2 NOPRIOR=1 MAPITER=0
$EST METHOD=IMP INTERACTION EONLY=1 NITER=20 ISAMPLE=3000 PRINT=1
NOABORT SIGL=6 NOPRIOR=1
$EST METHOD=BAYES INTERACTION NBURN=2000 NITER=5000 PRINT=10
FILE=example4.txt SIGL=6 NOPRIOR=0
$EST MAXEVAL=9999 NSIG=3 SIGL=10 PRINT=1
METHOD=CONDITIONAL INTERACTION NONINFETA=1
NOABORT FILE=example4.ext NOPRIOR=1
$COV MATRIX=R UNCONDITIONALPopulation Mixture Problem in 1 Compartment model, with rate constant parameter across sub-populations
Population Mixture Problem in 1 Compartment model, with rate constant parameter mean modeled for two sub-populations, but its inter-subject variance is the same in both sub-populations.
;Model Desc: Population Mixture Problem in 1
; Compartment model, with rate constant parameter
; mean modeled for two sub-populations, but its inter-subject
; variance is the same in both sub-populations
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# example5 (from ad1tr1m4t)
$INPUT C SET ID JID TIME CONC=DV DOSE=AMT RATE EVID MDV CMT
VC1 K101 VC2 K102 SIGZ PROB
$DATA example5.csv IGNORE=C
$SUBROUTINES ADVAN1 TRANS1
$MIX
P(1)=THETA(4)
P(2)=1.0-THETA(4)
NSPOP=2
$PK
Q=1
IF(MIXNUM.EQ.2) Q=0
MU_1=THETA(1)
; Note that MU_2 can be modeled as THETA(2) or THETA(3),
; depending on the MIXNUM value.
; Also, we are avoiding IF/THEN blocks.
MU_2=Q*THETA(2)+(1.0-Q)*THETA(3)
V=DEXP(MU_1+ETA(1))
K=DEXP(MU_2+ETA(2))
S1=V
$ERROR
Y = F + F*EPS(1)
$THETA
(-1000.0 4.3 1000.0) ;[MU_1]
(-1000.0 -2.9 1000.0) ;[MU_2-1]
(-1000.0 -0.67 1000.0) ;[MU_2-2]
(0.0001 0.667 0.9999) ;[P(1)]
$OMEGA BLOCK(2)
0.04 ;[p]
0.01 ;[f]
0.04 ;[p]
$SIGMA
0.01 ;[p]
$EST METHOD=ITS INTERACTION NITER=100 PRINT=1 NOABORT SIGL=8
FILE=example5.ext CTYPE=3
$EST METHOD=IMPMAP INTERACTION NITER=20 ISAMPLE=300 PRINT=1 NOABORT SIGL=8
$EST METHOD=IMP INTERACTION NITER=20 MAPITER=0 ISAMPLE=1000 PRINT=1 NOABORT
SIGL=6
$EST NBURN=500 NITER=500 METHOD=SAEM INTERACTION PRINT=10 SIGL=6 ISAMPLE=2
$EST METHOD=IMP INTERACTION NITER=5 ISAMPLE=1000 PRINT=1
NOABORT SIGL=6 EONLY=1
$EST METHOD=BAYES INTERACTION NBURN=2000 NITER=5000 PRINT=10
FILE=example5.txt SIGL=8
$EST MAXEVAL=9999 NSIG=2 SIGL=8 PRINT=10 FILE=example5.ext
METHOD=CONDITIONAL INTERACTION NOABORT
$COV MATRIX=RReceptor Mediated Clearance model with Dynamic Change in Receptors
Receptor Mediated Clearance model with Dynamic Change in Receptors.
;Model Desc: Receptor Mediated Clearance model with Dynamic Change
; in Receptors
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# example6 (from r2compl)
$INPUT C SET ID JID TIME DV=CONC DOSE=AMT RATE EVID MDV CMT
$DATA example6.csv IGNORE=C
; The new numerical integration solver is used, although ADVAN=9
; is also efficient for this problem.
$SUBROUTINES ADVAN13 TRANS1 TOL=4
$MODEL NCOMPARTMENTS=3
$PK
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
MU_5=THETA(5)
MU_6=THETA(6)
MU_7=THETA(7)
MU_8=THETA(8)
VC=EXP(MU_1+ETA(1))
K10=EXP(MU_2+ETA(2))
K12=EXP(MU_3+ETA(3))
K21=EXP(MU_4+ETA(4))
VM=EXP(MU_5+ETA(5))
KMC=EXP(MU_6+ETA(6))
K03=EXP(MU_7+ETA(7))
K30=EXP(MU_8+ETA(8))
S3=VC
S1=VC
KM=KMC*S1
F3=K03/K30
$DES
DADT(1) = -(K10+K12)*A(1) + K21*A(2) - VM*A(1)*A(3)/(A(1)+KM)
DADT(2) = K12*A(1) - K21*A(2)
DADT(3) = -VM*A(1)*A(3)/(A(1)+KM) - K30*A(3) + K03
$ERROR
CALLFL=0
ETYPE=1
IF(CMT.NE.1) ETYPE=0
IPRED=F
Y = F + F*ETYPE*EPS(1) + F*(1.0-ETYPE)*EPS(2)
$THETA
;Initial THETAs
( 4.0 ) ;[MU_1]
( -2.1 ) ;[MU_2]
( 0.7 ) ;[MU_3]
( -0.17 );[MU_4]
( 2.2 ) ;[MU_5]
( 0.14 ) ;[MU_6]
( 3.7 ) ;[MU_7]
( -0.7) ;[MU_8]
;Initial OMEGAs
$OMEGA BLOCK(8)
0.2 ;[p]
-0.0043 ;[f]
0.2 ;[p]
0.0048 ;[f]
-0.0023 ;[f]
0.2 ;[p]
0.0032 ;[f]
0.0059 ;[f]
-0.0014 ;[f]
0.2 ;[p]
0.0029 ;[f]
0.0027 ;[f]
-0.00026 ;[f]
-0.0032 ;[f]
0.2 ;[p]
-0.0025 ;[f]
0.00097 ;[f]
0.0024 ;[f]
0.00197 ;[f]
-0.0080 ;[f]
0.2 ;[p]
0.0031 ;[f]
-0.00571 ;[f]
0.0030 ;[f]
-0.0074 ;[f]
0.0025 ;[f]
0.0034 ;[f]
0.2 ;[p]
0.00973 ;[f]
0.00862 ;[f]
0.0041 ;[f]
0.0046 ;[f]
0.00061 ;[f]
-0.0056 ;[f]
0.0056 ;[f]
0.2 ;[p]
$SIGMA
0.1 ;[p]
0.1 ;[p]
$PRIOR NWPRI
; Omega prior
$OMEGAP BLOCK(8)
0.2 FIX
0.0 0.2
0.0 0.0 0.2
0.0 0.0 0.0 0.2
0.0 0.0 0.0 0.0 0.2
0.0 0.0 0.0 0.0 0.0 0.2
0.0 0.0 0.0 0.0 0.0 0.0 0.2
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.2
; degrees of freedom for OMEGA prior
$OMEGAPD
(8 FIXED) ;[dfo]
; Starting with a short iterative two stage analysis brings the
; results closer so less time needs to be spent during the
; burn-in of the BAYES analysis
$EST METHOD=ITS INTERACTION SIGL=4 NITER=15 PRINT=1
FILE=example6.ext NOABORT NOPRIOR=1
$EST METHOD=BAYES INTERACTION NBURN=4000 SIGL=4 NITER=10000
PRINT=10 CTYPE=3 FILE=example6.txt NOABORT NOPRIOR=0
; By default, ISAMPLE_M* are 2. Since there are many data points
; per subject, setting these to 1 is enough, and it reduces the
; time of the analysis
ISAMPLE_M1=1 ISAMPLE_M2=1 ISAMPLE_M3=1 IACCEPT=0.4
$COV MATRIX=R UNCONDITIONALInter-occasion Variability
;Model Desc: Interoccasion Variability
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB run# example7r
$INPUT C SET ID TIME AMT RATE EVID MDV CMT DV OCC
$ABBR REPLACE ETA(OCC_CL)=ETA(3,4,5)
$DATA example7r.csv IGNORE=C
$SUBROUTINES ADVAN1 TRANS2
$PK
MU_1=THETA(1)
MU_2=THETA(2)
V=DEXP(MU_1+ETA(1))
S1=V
VC=V
CL=DEXP(MU_2+ETA(2))*EXP(ETA(OCC_CL))
$ERROR
IPRED=F
Y = F+F*EPS(1)
;Initial THETAs
$THETA
2.0 ;[MU_1]
2.0 ;[MU_2]
;Initial OMEGAs
$OMEGA BLOCK(2)
.3 ;[p]
-.01 ;[f]
.3 ;[p]
$OMEGA BLOCK(1)
.1 ;[p]
$OMEGA BLOCK(1) SAME(2)
$SIGMA
0.1 ;[p]
$PRIOR NWPRI
; Degrees of freedom for Prior Omega blocks
$OMEGAPD (2.0 FIXED) (1.0 FIXED)
; Prior OMEGAs
$OMEGAP BLOCK(2)
.14 FIX
0.0 .125
$OMEGAP BLOCK(1) .0164 FIX
$OMEGAP BLOCK(1) SAME(2)
$EST METHOD=ITS INTERACTION FILE=example7r.ext NITER=10000
PRINT=5 NOABORT SIGL=8 CTYPE=3 CITER=10
NOPRIOR=1 CALPHA=0.05 NSIG=2
$EST METHOD=SAEM INTERACTION NBURN=30000 NITER=500 SIGL=8
ISAMPLE=2 PRINT=10 SEED=1556678 CTYPE=3
CITER=10 CALPHA=0.05 NOPRIOR=1
$EST METHOD=IMP INTERACTION EONLY=1 NITER=4 ISAMPLE=3000
PRINT=1 SIGL=10 NOPRIOR=1 MAPITER=0
$EST METHOD=BAYES INTERACTION FILE=example7r.txt NBURN=10000
NITER=10000 PRINT=100 CTYPE=3 CITER=10
CALPHA=0.05 NOPRIOR=0
$EST METHOD=COND INTERACTION MAXEVAL=9999 NSIG=3 SIGL=10 PRINT=5
NOABORT NOPRIOR=1
FILE=example7r.ext
$COV MATRIX=R PRINT=E UNCONDITIONALSample History of Individual Values in MCMC Bayesian Analysis
;Model Desc: Two compartment Model, Using ADVAN3, TRANS4
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# Example 8 (from samp5l)
$abbr DECLARE INTEGER FIRST_WRITE
$abbr DECLARE INTEGER FIRST_WRITE2
$INPUT C SET ID JID TIME DV=CONC AMT=DOSE RATE EVID MDV CMT CLX V1X QX V2X
SDIX SDSX
$DATA example8.csv IGNORE=C
$SUBROUTINES ADVAN3 TRANS4
$PK
include nonmem_reserved_general
; Request extra information for Bayesian analysis.
; An extra call will then be made
; for accepted samples
BAYES_EXTRA_REQUEST=1
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
CL=DEXP(MU_1+ETA(1))
V1=DEXP(MU_2+ETA(2))
Q=DEXP(MU_3+ETA(3))
V2=DEXP(MU_4+ETA(4))
S1=V1
; When Bayes_extra=1, then this particular set of individual
; parameters were "accepted"
; So you may record them if you wish
IF(BAYES_EXTRA==1 .AND. ITER_REPORT>=0 .AND. FIRST_WRITE==0 .AND.
NEWIND/=2) THEN
IF(PNM_RUN_MODE==PNM_SINGLE) THEN
" OPEN(UNIT=53,FILE='example8.prr')
; If parallel, need to have different file name for each node
; best if a pnm file with -awnf is used, such as mpiwini8_awnf.pnm,
; and command line option -maxlim=3. This allows all files to be
; written to primary run directory, and not into worker folders.
ELSE
" OPEN(UNIT=53,FILE='example8_'//trim(tfi(pnm_node_number))//'.prr')
ENDIF
IF(PNM_NODE_NUMBER==1) THEN ; Only first file should have header,
; for easy merging later
; When a parallel run is completed you can merge the files with
; (Windows):
; copy example8_*.prr example.prr
; (Linux):
; cat example8_*.prr > example.prr
" WRITE(53,94) 'ITER_REPORT ',' ID ',' CL ',' V1
',' Q ',' V2 ',' OBJI '
" 94 FORMAT(A12,1X,A14,5(1X,A12))
ENDIF
FIRST_WRITE=1
ENDIF
IF(BAYES_EXTRA==1 .AND. ITER_REPORT>=0 .AND. NEWIND/=2 ) THEN
" WRITE(53,95) ITER_REPORT,ID,CL,V1,Q,V2,OBJI(NIREC,1)
" 95 FORMAT(I12,1X,F14.0,5(1X,1PG12.5))
ENDIF
$ERROR
include nonmem_reserved_general
BAYES_EXTRA_REQUEST=1
Y = F + F*EPS(1)
IF(BAYES_EXTRA==1 .AND. ITER_REPORT>=0 .AND. FIRST_WRITE2==0 ) THEN
IF(PNM_RUN_MODE==PNM_SINGLE) THEN
" OPEN(UNIT=54,FILE='example8.iwr')
ELSE
" OPEN(UNIT=54,FILE='example8_'//trim(tfi(pnm_node_number))//'.iwr')
ENDIF
IF(PNM_NODE_NUMBER==1) THEN ; Only one file should have header,
; for easy merging later
; When a parallel run is completed you can merge the files with
; (Windows):
; copy example8_*.iwr example.iwr
; (Linux):
; cat example8_*.iwr > example.iwr
" WRITE(54,96) 'ITER_REPORT ','ID ','TIME ','F '
" 96 FORMAT(A12,1X,A14,2(1X,A12))
ENDIF
FIRST_WRITE2=1
ENDIF
IF(BAYES_EXTRA==1 .AND. ITER_REPORT>=0 ) THEN
" WRITE(54,97) ITER_REPORT,ID,TIME,F
" 97 FORMAT(I12,1X,F14.0,2(1X,1PG12.5))
ENDIF
; Initial values of THETA
$THETA
(2.0) ;[LN(CL)]
(2.0) ;[LN(V1)]
(2.0) ;[LN(Q)]
(2.0) ;[LN(V2)]
;INITIAL values of OMEGA
$OMEGA BLOCK(4)
0.15 ;[P]
0.01 ;[F]
0.15 ;[P]
0.01 ;[F]
0.01 ;[F]
0.15 ;[P]
0.01 ;[F]
0.01 ;[F]
0.01 ;[F]
0.15 ;[P]
;Initial value of SIGMA
$SIGMA
(0.6 ) ;[P]
$PRIOR NWPRI
; Prior information to the THETAS.
$THETAP (2.0 FIX) (2.0 FIX) (2.0 FIX) (2.0 FIX)
$THETAPV BLOCK(4)
10000 FIX
0.00 10000
0.00 0.00 10000
0.00 0.00 0.0 10000
; Prior information to the OMEGAS.
$OMEGAP BLOCK(4)
0.2 FIX
0.0 0.2
0.0 0.0 0.2
0.0 0.0 0.0 0.2
$OMEGAPD (4 FIX)
$EST METHOD=BAYES INTERACTION FILE=example8.ext NBURN=10000 NITER=1000
PRINT=100 NOPRIOR=0 RANMETHOD=P
CTYPE=3 CINTERVAL=100When BAYES_EXTRA=1, this means that the sample was accepted, and you can then store relevant parameters into a file.
Note that the distinct files names are given for each node, so there are no name conflicts during parallelization.
Note the use of the include file nonmem_reserved_general, which for purposes of this example contain the following declarations of reserved variables:
"C ITER_REPORT: Iteration number that is reported to output
"C (can be negative, if during a burn period).
"C BAYES_EXTRA, BAYES_EXTRA_REQUEST, used in example 8
" USE NMBAYES_REAL, ONLY: OBJI
" USE NMBAYES_INT, ONLY: ITER_REPORT,BAYES_EXTRA_REQUEST,BAYES_EXTRA
" USE PNM_CONFIG, ONLY: PNM_NODE_NUMBER
" USE NM_INTERFACE, ONLY: TFI,TFD
" USE PNM_CONFIG, ONLY: PNM_RUN_MODE
" USE PNM_CONFIG, ONLY: PNM_SINGLE
" USE PNM_CONFIG, ONLY: PNM_MANAGER
" USE PNM_CONFIG, ONLY: PNM_WORKER
" USE PNM_CONFIG, ONLY: PNM_NODE_NUMBERA suitable command could be (example8.bat):
nmfe75 example8.ctl example8.res -prdefault -parafile=mpiwini8_awnf.pnm
[nodes]=3 -maxlim=3 -parafprint=100Simulated Annealing For SAEM using Constraint Subroutine
;Model Desc: Two compartment Model, Using ADVAN3, TRANS4
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# Example 9 (from samp5l)
$INPUT C SET ID JID TIME DV=CONC AMT=DOSE RATE EVID MDV CMT CLX V1X QX V2X
SDIX SDSX
$DATA example9.csv IGNORE=C
$SUBROUTINES ADVAN3 TRANS4 OTHER=ANEAL.F90
$PK
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
CL=DEXP(MU_1+ETA(1))
V1=DEXP(MU_2+ETA(2))
Q=DEXP(MU_3+ETA(3))
V2=DEXP(MU_4+ETA(4))
S1=V1
$ERROR
Y = F + F*EPS(1)
; Initial values of THETA
$THETA
(0.001, 2.0) ;[LN(CL)]
(0.001, 2.0) ;[LN(V1)]
(0.001, 2.0) ;[LN(Q)]
(0.001, 2.0) ;[LN(V2)]
;INITIAL values of OMEGA
$OMEGA BLOCK(4)
0.05 ;[P]
0.01 ;[F]
0.05 ;[P]
0.01 ;[F]
0.01 ;[F]
0.05 ;[P]
0.01 ;[F]
0.01 ;[F]
0.01 ;[F]
0.05 ;[P]
;Initial value of SIGMA
$SIGMA
(0.6 ) ;[P]
$EST METHOD=SAEM INTERACTION FILE=example9.ext NBURN=5000 NITER=500 PRINT=10
NOABORT SIGL=6
CTYPE=3 CINTERVAL=100 CITER=10 CALPHA=0.05File Aneal.f90:
SUBROUTINE
CONSTRAINT(THETAS,NTHETAS,SIGMA2,NSIGMAS,OMEGA,NOMEGAS,ITER_NO)
USE SIZES, ONLY: ISIZE,DPSIZE
INCLUDE '..\nm\TOTAL.INC'
INTEGER(KIND=ISIZE) NTHETAS,NSIGMAS,NOMEGAS,ITER_NO
INTEGER I,J,ITER_OLD
DATA ITER_OLD /-1/
REAL(KIND=DPSIZE) ::
OMEGA(MAXOMEG,MAXOMEG),THETAS(MAXPTHETA),SIGMA2(MAXPTHETA)
REAL(KIND=DPSIZE) :: OMEGO(MAXOMEG)
SAVE
!----------------------------------------------------------------------------
-
IF(SAEM_MODE==1 .AND. IMP_MODE==0 .AND. ITS_MODE==0 .AND. ITER_NO<200)
THEN
IF(ITER_NO/=ITER_OLD .OR. ITER_NO==0) THEN
! During burn-in phase of SAEM, and when a new iteration occurs
! (iter_old<>iter_no)
! store the present diagonals of OMEGAs
ITER_OLD=ITER_NO
DO I=1,NOMEGAS
OMEGO(I)=OMEGA(I,I)
ENDDO
ENDIF
IF(ITER_NO /=0) THEN
DO I=1,NOMEGAS
! Use whatever algorithm needed to "slow down" the reduction of Omega
! The expansion of Omega should be less with each iteration.
OMEGA(I,I)=OMEGO(I)*(1.0D+00+10.0D+00/ITER_NO)
ENDDO
ENDIF
ENDIF
RETURN
!
END SUBROUTINE CONSTRAINTOne Compartment First Order Absorption Pharmacokinetics with Categorical Data
$PROB F_FLAG04est2a.ctl
$INPUT C ID DOSE=AMT TIME DV WT TYPE
$DATA example10.csv IGNORE=@
$SUBROUTINES ADVAN2 TRANS2
$PK
CALLFL=1
MU_1=DLOG(THETA(1))
KA=DEXP(MU_1+ETA(1))
MU_2=DLOG(THETA(2))
V=DEXP(MU_2+ETA(2))
MU_3=DLOG(THETA(3))
CL=DEXP(MU_3+ETA(3))
SC=V/1000
$THETA 5.0 10.0 2.0 0.1 0.1
$OMEGA BLOCK (3)
0.5
0.01 0.5
0.01 0.01 0.5
; Because THETA(4) and THETA(5) have no inter-subject variability
; associated with them, the algorithm must use a more computationally
; expensive gradient evaluation for these two parameters
$SIGMA 0.1
$PRIOR NWPRI
; Priors to OMEGAs
$OMEGAP BLOCK (3)
0.09 FIX
0.0 0.09
0.0 0.0 0.09
$OMEGAPD (3 FIX)
$ERROR
EXPP=THETA(4)+F*THETA(5)
IF (TYPE.EQ.0) THEN
; PK Data
F_FLAG=0
Y=F+F*ERR(1) ; a prediction
ELSE
; Categorical data
F_FLAG=1
; Use protected exponent PEXP, to avoid numerical overflow
A=PEXP(EXPP)
B=1+A
Y=DV*A/B+(1-DV)/B ; a likelihood
ENDIF
$EST METHOD=ITS INTER LAP NITER=1000 PRINT=5 SIGL=6 NSIG=2
NOABORT NOPRIOR=1 CTYPE=3 CITER=10 CALPHA=0.05
FILE=example10.ext
; Because of categorical data, which can make conditional density highly
; non-normal, select a t-distribution with 4 degrees of freedom for
; importance sampling proposal density
$EST METHOD=IMP INTER LAP NITER=1000 PRINT=1 ISAMPLE=300 DF=4
IACCEPT=1.0
$EST METHOD=IMP EONLY=1 NITER=5 ISAMPLE=1000 PRINT=1 DF=4
IACCEPT=1.0 MAPITER=0
$EST METHOD=SAEM EONLY=0 INTER LAP NBURN=2000 NITER=1000 PRINT=50
DF=0 IACCEPT=0.4
$EST METHOD=IMP EONLY=1 NITER=5 ISAMPLE=1000 PRINT=1 DF=4
IACCEPT=1.0 MAPITER=0
$EST METHOD=BAYES NBURN=3000 NSAMPLE=3000 PRINT=100
FILE=example10.txt DF=0 IACCEPT=0.4 NOPRIOR=0
$EST METHOD=COND LAP INTER MAXEVAL=9999 PRINT=1 FILE=example10.ext
NOPRIOR=1 NOHABORT
$COV UNCONDITIONAL PRINT=E MATRIX=R SIGL=10
$TABLE ID DOSE WT TIME TYPE DV A NOPRINT FILE=example10.tab