129 Examples
Generated from nmhelp plain-text source files (*.exa) of the NONMEM 7.6.0 distribution. These mirror the content presented by nmhelp and complement the formatted Part VIII Help Guide.
1 EXAMPLES
2 ABBREVIATED FUNCTION EXAMPLES
Here is an example of the use and construction of a function of both time and the first-order elimination rate constant of a simple PK one-compartment model, to compute the amount in the system at a given time. (See abbreviated_function).
$PRED
...
VECTRA(1)=THETA(1)*EXP(ETA(1))
VECTRA(2)=TIME
...
A=DOSE*FUNCA(VECTRA)
...
With the $ABBREVIATED REPLACE feature, the function name can be
changed to one that is perhaps easier to understand.
$ABBR REPLACE ONE_COMPARTMENT_LINEAR=FUNCA
...
$PRED
...
VECTRA(1)=THETA(1)*EXP(ETA(1))
VECTRA(2)=TIME
...
A=DOSE*ONE_COMPARTMENT_LINEAR(VECTRA)
...
Either way, one needs to supply the function.
FUNCTION FUNCA(X,X1,X2)
! Implements the One Compartment Linear model
USE SIZES, ONLY: DPSIZE
REAL(KIND=DPSIZE), INTENT(IN) :: X
REAL(KIND=DPSIZE), INTENT(IN OUT) :: X1,X2
REAL(KIND=DPSIZE) :: FUNCA
DIMENSION :: X(9),X1(9),X2(9,9)
REAL(KIND=DPSIZE) :: EXPT
! THE FUNCTION ITSELF
EXPT=EXP(-X(1)*X(2))
FUNCA=EXPT
! 1ST. PARTIALS
X1(1)=-EXPT*X(2)
X1(2)=-EXPT*X(1)
! 2ND. PARTIALS
X2(1,1)=EXPT*X(2)*X(2)
X2(1,2)=EXPT*X(1)*X(2)
X2(2,1)=EXPT*X(1)*X(2)
X2(2,2)=EXPT*X(1)*X(1)
RETURN
END
Suppose the above code is in file ONE_COMPARTMENT_LINEAR.f90.
Then the control stream should also include
$SUBROUTINES OTHER=ONE_COMPARTMENT_LINEAR.f90
Multiple functions can be used with the $ABBREVIATED REPLACE feature.
Suppose two functions are needed, called (for example) RED and GREEN.
The control stream might contain the following fragment:
$SUBR OTHER=red.f90 OTHER=green.f90
...
$ABBR REPLACE RED=FUNCA
$ABBR REPLACE GREEN=FUNCB
...
$PRED
...
VECTRA(1)=THETA(1)*EXP(ETA(1))
VECTRA(2)=TIME
A=DOSE*RED(VECTRA)
B=DOSE*GREEN(VECTRA)
...
The file red.f90 contains code for FUNCA, and the file green.f90 con-
tains code for FUNCB.
The $ABBREVIATED REPLACE feature could also be used to provide more
meaningful names for reserved vectors. For example, suppose the argu-
ments of the functions are different vectors.
$ABBR REPLACE RED=FUNCA, REDARG=VECTRA
$ABBR REPLACE GREEN=FUNCB, GREENARG=VECTRB
...
$PRED
...
REDARG(1)=...
REDARG(2)=...
GREENARG(1)=...
GREENARG(2)=...
A=DOSE*RED(REDARG)
B=DOSE*GREEN(GREENARG)
...REFERENCES: Guide Introduction_7
3 BAYES EXAMPLE 1
This is example1.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;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.ETAREFERENCES: Guide Introduction_7
4 BAYES EXAMPLE 10
This is example10.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
$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.tabREFERENCES: Guide Introduction_7
5 BAYES EXAMPLE 10L
This is example10l.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
$PROB F_FLAG04est2a.ctl
$INPUT C ID DOSE=AMT TIME DV WT TYPE
$DATA example10l.csv IGNORE=@
$SUBROUTINES ADVAN2 TRANS2
$PK
CALLFL=1
MU_1=THETA(1)
KA=DEXP(MU_1+ETA(1))
MU_2=THETA(2)
V=DEXP(MU_2+ETA(2))
MU_3=THETA(3)
CL=DEXP(MU_3+ETA(3))
SC=V/1000
$THETA 1.6 2.3 0.7 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=example10l.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=example10l.txt DF=0 IACCEPT=0.4 NOPRIOR=0
$EST METHOD=COND LAP INTER MAXEVAL=9999 PRINT=1 FILE=example10l.ext
NOPRIOR=1
$COV UNCONDITIONAL PRINT=E MATRIX=R SIGL=10
$TABLE ID DOSE WT TIME TYPE DV A NOPRINT FILE=example10l.tabREFERENCES: Guide Introduction_7
6 BAYES EXAMPLE 1B
This is example1B.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;Model Desc: Two compartment Model, Using ADVAN3, TRANS4
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# Example 1B (from samp5l)
$INPUT C SET ID JID TIME DV=CONC AMT=DOSE RATE EVID MDV CMT CLX
V1X QX V2X SDIX SDSX
$DATA example1b.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]
;NTHETA=number of Thetas to be estimated
;NETA=number of Etas to be estimated
; (and to be described by NETAxBETA OMEGA matrix)
;NTHP=number of thetas which have a prior
;NETP=number of Omegas with prior
;Prior information is important for MCMC Bayesian analysis,
; not necessary for maximization methods
$PRIOR NWPRI NTHETA=4, NETA=4, NTHP=4, NETP=4
; Prior information of THETAS
$THETAP (2.0 FIX) (2.0 FIX) (2.0 FIX) (2.0 FIX)
; Variance to prior information of THETAS.
; Because variances are very large, this
; means that prior information to the THETAS is highly uninformative.
$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
; 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)
; 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
; example1b.ext
$EST METHOD=ITS MAPITER=0 INTERACTION FILE=example1b.ext NITER=500
PRINT=5 CINTERVAL=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
; example1b.ext.
$EST METHOD=SAEM INTERACTION NBURN=3000 NITER=500
PRINT=250 CINTERVAL=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 MAPITER=0 NITER=5
ISAMPLE=3000 PRINT=1 SIGL=8 NOPRIOR=1
$EST METHOD=DIRECT INTERACTION EONLY=1 NITER=3 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=example1b.txt NBURN=10000
NITER=10000 PRINT=1000 CINTERVAL=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 example1b.ext with the EM methods is fine.
$EST METHOD=COND INTERACTION MAXEVAL=9999 NSIG=3 SIGL=10
PRINT=5 NOABORT NOPRIOR=1
FILE=example1b.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=example1b.TAB NOPRINT
$TABLE ID CL V1 Q V2 FIRSTONLY NOAPPEND NOPRINT FILE=example1b.PAR
FORMAT=,1PE15.8
$TABLE ID ETA1 ETA2 ETA3 ETA4 FIRSTONLY NOAPPEND NOPRINT
FILE=example1b.ETA FORMAT=,1PE15.8REFERENCES: Guide Introduction_7
7 BAYES EXAMPLE 2
This is example2.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;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 UNCONDITIONALREFERENCES: Guide Introduction_7
8 BAYES EXAMPLE 3
This is example3.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;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-popluation
$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 UNCONDITIONALREFERENCES: Guide Introduction_7
9 BAYES EXAMPLE 4
This is example4.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;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)
V=DEXP(MU_1+ETA(1))
K10M=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
(-1000.0 4.3 1000.0) ;[MU_1]
(-1000.0 -2.9 1000.0) ;[MU_2]
(-1000.0 -0.67 1000.0) ;[MU_3]
(0.0001 0.667 0.9999) ;[P(1)]
$OMEGA BLOCK(3)
.04 ;[p]
0.01 ;[f]
.027 ;[p]
0.01 ;[f]
0.001 ;[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=12 PRINT=1
METHOD=CONDITIONAL INTERACTION
NOABORT FILE=example4.ext NOPRIOR=1
$COV MATRIX=R UNCONDITIONAL SIGL=10REFERENCES: Guide Introduction_7
10 BAYES EXAMPLE 5
This is example5.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;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 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 MAPITER=0
$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=RREFERENCES: Guide Introduction_7
11 BAYES EXAMPLE 6
This is example6.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;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 UNCONDITIONALREFERENCES: Guide Introduction_7
12 BAYES EXAMPLE 6S
This is example6s.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory. Note the old style syntax for priors.
;Model Desc: Receptor Mediated Clearance model with Dynamic Change
; in Receptors
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# example6s (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
$PRIOR NWPRI NTHETA=8, NETA=8, NTHP=0, NETP=8, NPEXP=1
$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
DEL=1.0E-10
DADT(1) = -(K10+K12)*A(1) + K21*A(2) - VM*A(1)*A(3)/(A(1)+KM+DEL)
DADT(2) = K12*A(1) - K21*A(2)
DADT(3) = -VM*A(1)*A(3)/(A(1)+KM+DEL) - 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]
; degrees of freedom for OMEGA prior
(8 FIXED) ;[dfo]
;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]
; Omega prior
$OMEGA 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
$SIGMA BLOCK(1) 0.1 ;[p]
;$SIGMA BLOCK(1) 0.1 ;[p]
$SIGMA BLOCK(1) SAME
; 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=30 PRINT=1
FILE=example6s.ext NOABORT NOPRIOR=1 CTYPE=3
$EST METHOD=IMP INTERACTION SIGL=4 NITER=30 PRINT=1
FILE=example6s.ext NOABORT NOPRIOR=1 ISAMPLE=100
$EST METHOD=BAYES INTERACTION NBURN=4000 SIGL=4 NITER=1000
PRINT=10 CTYPE=3 FILE=example6s.txt NOABORT NOPRIOR=0
PSAMPLE_M1=1 PSAMPLE_M2=1 PSAMPLE_M3=1
; 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 UNCONDITIONALREFERENCES: Guide Introduction_7
13 BAYES EXAMPLE 6SB
This is example6sb.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory. Note the old style syntax for priors.
;Model Desc: Receptor Mediated Clearance model with Dynamic Change
; in Receptors
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# example6sb (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
$PRIOR NWPRI NTHETA=8, NETA=8, NTHP=0, NETP=8, NPEXP=1
$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(1)
$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]
; degrees of freedom for OMEGA prior
(8 FIXED) ;[dfo]
;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]
; Omega prior
$OMEGA 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
$SIGMA BLOCK(1) 0.1 ;[p]
; 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=30 PRINT=1
FILE=example6sb.ext NOABORT NOPRIOR=1 CTYPE=3
$EST METHOD=IMP INTERACTION SIGL=4 NITER=30 PRINT=1
FILE=example6sb.ext NOABORT NOPRIOR=1 ISAMPLE=100
$EST METHOD=BAYES INTERACTION NBURN=4000 SIGL=4 NITER=1000 PRINT=10
CTYPE=3 FILE=example6sb.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 UNCONDITIONALREFERENCES: Guide Introduction_7
14 BAYES EXAMPLE 7R
This is example7r.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;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 UNCONDITIONALREFERENCES: Guide Introduction_7
15 BAYES EXAMPLE 8
This is example8.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;Model Desc: Two compartment Model, Using ADVAN3, TRANS4
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# Example 8 (from samp5l)
$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. TIME==0.0) THEN
" WRITE(51,98) ITER_REPORT,ID,CL,V1,Q,V2
" 98 FORMAT(I12,1X,F14.0,4(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 ) THEN
" WRITE(52,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)x4
$THETAPV BLOCK(4) FIX VALUES(10000.0,0.0)
; 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 CTYPE=3 CINTERVAL=100REFERENCES: Guide Introduction_7
16 BAYES EXAMPLE 9
This is example9.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;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.05REFERENCES: Guide Introduction_7
17 BIND EXAMPLE
EXAMPLES OF THE USE OF $BIND Case 1. Suppose a fragment of the dose records looks like this: TIME AMT PREP X 0 100 4 . 10 100 1 90 20 0 0 80 Suppose $INPUT and $PK include this code:
$INPUT ID TIME AMT DV PREP X
$PK
F1=1
IF (PREP.EQ.4) F1=THETA(3)
ALAG1=THETA(5)
CL=THETA(1)*X*EXP(ETA(1))
It is intended that theta(3) be the ratio of bioavailabilities of
preparation 4 to preparation 1. However, calls to PK are made only at
event times 0, 10 and 20. Suppose the value of theta(5) (lag time) is
5. After 5 time units, the first dose enters the system. The value
of PREP (dose preparation type) is taken from the event record at time
10, which is 1. F1 is computed incorrectly. If the estimation step
is run, the gradient with theta(3) will be 0, indicating the problem.
Case 2.
If lag time is not modeled (no ALAG1), but the doses are multiple
(ADDL>0, II>0), a more subtle error can occur. Suppose a fragment of
the data records looks like this:
TIME AMT PREP X II ADDL
0 100 4 . 5 1
10 100 1 90 5 1
20 0 0 80 . .
Suppose $INPUT and $PK include this code:
$INPUT ID TIME AMT DV PREP X II ADDL
$PK
F1=1
IF (PREP.EQ.4) F1=THETA(3)
CL=THETA(1)*X*EXP(ETA(1))
F1 is computed correctly for the initiating dose records at times 0
and 10, and for the additional dose at time 15, but incorrectly for
the additional dose at time 5. Theta(3) applies only to the first
dose having PRED=4. (The gradient for theta(3) is not 0.) The same
problem could have occurred with the dose having PREP=1, except that
by an accident of coding, PREP=0 (on the record having TIME=20) is
treated the same as PREP=1, so that the additional dose at time 15 has
F1=1 as it should.
On the other hand, the values of X are presumably recorded so that 90
is valid for the time 0 through 10 and 80 is valid for the time 10
through 20. Therefore, the model for CL is correct in both examples.
In either of these two examples, the model can be corrected by the
insertion of $BIND and CALLFL=-2:
$INPUT ID TIME AMT DV PREP X II ADDL ; II & ADDL only in case 2
$BIND - - - - DOSE NEXT
$PK
CALLFL=-2
Because of CALLFL=-2 in the abbreviated code, additional calls to PK
are made at the nonevent dose times (5 in case 1; 5 and 15 in case 2).
Because of the $BIND record, PREP has the value at these calls from
the initiating dose record (4 at time 5, 1 at time 15). F1 is com-
puted appropriately for the dose preparation. X still has values from
the next event record, so the computation of CL is unchanged, e.g., it
is the same from time=0 to time=5 as it is from time=5 to time=10.
(See $BIND).REFERENCES: Guide IV, section V.C.2 , V.C.5; Guide VI, section III.B.2 , V.K
18 BIVARIATE FUNCTION
This is a fully-worked out example of the use of the $ABBREVIATED FUNCTION option to declare a user-defined function. In this example, the function BIVARIATE is used to compute the integral of the bivariate normal distribution of two correlated data points, DV1 and DV2, with correlation RHO. This is bivariate.ctl from the NONMEM 7 distribution medium. This file, the data file bivariate.csv, and the function’s source code bivariate.f90 can be found in the examples directory.
$PROB BIVARIATE EXAMPLE
; THESE DECLARATIONS ALLOW ANY FUNCTION TO HAVE
; ALTERNATIVE DIMENSIONS FOR THEIR ARRAYS
; BUT, USER DEFINED DIMENSIONS ARE PASSED AS THE
; LAST ARGUMENT TO FUNC, SUCH AS:
; BV=BIVARIATE(VBI(1),FNC001_1(1,1),FNC001_2(1,1,1),5)
$ABBR FUNCTION BIVARIATE(VBI,5)
$INPUT SIM ID DOSE DV TIME
$DATA bivariate.csv IGNORE=
$SUBROUTINES OTHER=bivariate.f90
$PRED
B1=THETA(1)
B2=THETA(2)
B3=THETA(3)
K =LOG(2)/EXP(THETA(4))
ED50=EXP(THETA(5))
U =(1-EXP(-K*TIME ))
MU_1=B1+B3*DOSE/(DOSE+ED50)
MU_2=B2
MXB=MU_1+ETA(1)
MXU=MU_2+ETA(2)
MX =MXB + MXU*U ;***Current model prediction***;
PHIMX=PHI(MX)
IF(NEWIND.NE.2) THEN
TIMEP=0
MXP=0
DVP=0
PHIMXP=0.5
ENDIF
RHOB=(2/(1+EXP(-THETA(8)))-1)
IF(RHOB>0.0) RHO=RHOB**(TIME-TIMEP)
IF(RHOB==0.0) RHO=0.0
IF(RHOB<0.0) RHO=-(-RHOB)**(TIME-TIMEP)
PC =(1-PHIMX) *(1-DV ) + PHIMX*DV
IF(PC.LE.0.0) EXIT
V=SQRT(1+OMEGA(1,1)+OMEGA(2,2)*U**2)
POPP = (B1+B2*U +B3*DOSE/(DOSE+ED50))/V ;*Population mean prediction*;
IF (TIME.EQ.1) THEN
JP=PC
PCP=1.0
ELSE
;***Pass information to bivariate normal***;
VBI(1)=RHO
VBI(2)=MX
VBI(3)=MXP
VBI(4)=1 ;*0 = Upper tail as in Drezner & Wesolowsky; 1 = Bottom tail*;
VBI(5)=1 ;*0 = 3 pt approximation; 1 = 5 point approximation*;
BV=BIVARIATE(VBI)
JP=((DV-1)*(DVP-1)+(DV-1)*(1-2*DVP)*PHIMXP+(DVP-1) &
*(1-2*DV)*PHIMX+(1-2*DV)*(1-2*DVP)*BV)
ENDIF
IF(JP.LE.0.0) EXIT
LOGL=LOG(JP/PCP)
Y = -2*LOGL
MXP=MX
PCP=PC
DVP=DV
TIMEP=TIME
PHIMXP=PHIMX
$THETA
-1.7 ; 1 B1
1.2 ; 2 B2
2.9 ; 3 B3
1.4 ; 4 LOG(B4)
1.2 ; 5 LOG(B5)
(0.0 FIXED) ; 6 LOG SQRT VAR(ETA1)
(0.0 FIXED) ; 7 LOG SQRT VAR(ETA2)
2.2 ; 8 RHO parameter
$OMEGA DIAGONAL(2)
0.8 ; V1
0.8 ; V2
;$EST METHOD=IMP LAPLACE -2LL PRINT=1 NITER=300 ISAMPLE=300
;SIGL=6 CTYPE=3 NOHABORT
$EST MAX=0 PRINT=1 METHOD=1 LAPLACE -2LL SIGL=10 NOHABORT
$COV COMPRESS MATRIX=R PRINT=E UNCONDITIONAL
Here is a fragment of the data file bivariate.csv:
SIM,ID,DOSE,DV,TIME
1,1,0,0,1
1,1,0,0,2
1,1,0,0,4
1,1,0,0,8
1,1,0,0,16
1,1,0,0,24
1,1,0,0,36
This is the file bivariate.f90. (The extension .f90 in the file name
is not actually required because the contents of the file is copied
to FSUBS and compiled as part of FSUBS. It is not compiled indepen-
dently.)
! Integral of Bivariate normal distribution of two correlated data points,
! DV1, DV2, correlation RHO.
! BV1, BV2 and GAUSS FROM Drezner and Wesolowsky,
! J. Statistical Computat. Simul. 35, pp. 101-107, 1990.
! BV2 is more accurate for extreme values of rho
! Partial derivatives worked out by Robert Bauer.
!
! How to use in NONMEM 7.4
! VBI is a suggested vector name. Could be any name.
! But the function name must match what is in bivariate.f90
! $ABBR FUNCTION BIVARIATE(VBI,5)
! $SUBROUTINES ... OTHER=bivariate.f90
! $PK
! ...
! DV1 AND DV2 SHOULD BE TWO SEPARATE DATA ITEMS ON EACH RECORD.
! VBI(1)=RHO
! VBI(2)=DV1 ! H
! VBI(3)=DV2 !K
! VBI(4)=DTYPE !=0 FOR H TO INF, K TO INF, OR 1 FOR -INF TO H, -INF TO K
! VBI(5)=BVTYPE ! BVTYPE=0 USES SIMPLE INTEGRATOR, BVTYPE=1 USES MORE
! ACCURATE ONE FOR LARGE VALUES OF RHO
! F_FLAG=1
! IT IS SAFER TO AVOID BV BEING EXACTLY 0, JUST MAKE SURE
! IT IS REALLY SMALL, TO AVOID LOG() ERRORS
! BV=BIVARIATE(VBI)+1.0D-30
! IF YOU USE -2LL FORMAT (DEFAULT):
! Y=-2.0D+00*LOG(BV)
FUNCTION BIVARIATE(X,X1,X2,NDIM)
! NDIM SHOULD EQUAL 5
! IMPLICIT REAL*8 (A-H,O-Z)
IMPLICIT NONE
INTEGER NDIM
REAL*8 X(NDIM),X1(NDIM),X2(NDIM,NDIM)
REAL*8 R,H,K,ITYPE,PI,PI2,SPI2,R1,RSQR,HRK,KRH,PHRK,PKRH,XH,XK
REAL*8 XHRK,XKRH,XC,BV,BV2,BIVARIATE,GAUSS,BV1
EXTERNAL GAUSS,BV1,BV2
INTEGER BVTYPE
! unload the arguments
R=X(1)
H=X(2)
K=X(3)
ITYPE=X(4)
BVTYPE=X(5)
! ITYPE=0 INTGRATE H TO INF, K TO INF
! ITYPE=1 INTGRATE -INF TO H -INF TO K
! BVTYPE=0 USES THE SIMPLER BV1
! BVTYPE=1 USES THE MORE ACCURATE BV2
! SET UP SOME USEFUL CALCUATIONS FOR FIRST
! AND SECOND DERIVATIVES OF THE ARGUMENTS
! DERIVATIVES NOT NEEDED IF YOU USE IMP, SAEM, BAYES, OR IF
! YOU USE LAPLACE OR ITS WITH OPTMAP=1 ETADER=3 (NONMEM 7.3), OR
! DERIVATIVES NEEDED IF YOU USE STANDARD LAPLACE
PI=3.141592653589793238D+00
PI2=2.0D+00*pi
spi2=sqrt(pi2)
R1=1.0D+00-R*R
RSQR=SQRT(ABS(R1))
IF(RSQR.LE.0.0D+00) RSQR=1.0D-100
IF(R1.LE.0.0D+00) R1=1.0D-100
HRK=(H-R*K)/RSQR
KRH=(K-R*H)/RSQR
PHRK=1.0D+00-GAUSS(HRK)
PKRH=1.0D+00-GAUSS(KRH)
XH=EXP(-H*H/2.0D+00)/SPI2
XK=EXP(-K*K/2.0D+00)/SPI2
XHRK=EXP(-HRK*HRK/2.0D+00)/SPI2/RSQR
XKRH=EXP(-KRH*KRH/2.0D+00)/SPI2/RSQR
XC=EXP(-(H*H-2.0D+00*R*H*K+K*K)/2.0D+00/R1)/RSQR/PI2
! parital F WRT RHO, from Drezner and Wesolowsky, equation (4)
X1(1)=XC
! parital F WRT H
X1(2)=XH*PKRH+(ITYPE-1.0D+00)*XH
! parital F WRT K
X1(3)=XK*PHRK+(ITYPE-1.0D+00)*XK
! 2ND parital F WRT RHO,RHO
X2(1,1)=R/R1*XC+XC*(H*K*(1.0D+00-R)*(1.0D+00-R)-R*(H-K)*(H-K))/R1/R1
! 2ND parital F WRT RHO,H
X2(1,2)=XC*(R*K-H)/R1
! 2ND parital F WRT RHO,K
X2(1,3)=XC*(R*H-K)/R1
X2(2,1)=X2(1,2)
X2(3,1)=X2(1,3)
! 2ND parital F WRT H,K
X2(2,3)=XC
X2(3,2)=X2(2,3)
! 2ND parital F WRT H,H
X2(2,2)=-X1(2)*H-XC*R
! 2ND parital F WRT K,K
X2(3,3)=-X1(3)*K-XC*R
! write(50,*) 'A ',x1(1),x1(2),x1(3)
! write(50,*) 'B ',x2(1,1),x2(1,2),x2(1,3),x2(2,1), &
! x2(2,2),x2(2,3),x2(3,1),x2(3,2),x2(3,3)
! bV2 is more accurate for exterme rho values.
IF(BVTYPE==0) THEN
BV=BV1(H,K,R)
ELSE
BV=BV2(H,K,R)
ENDIF
IF(ITYPE==1.0D+00) BV=BV-GAUSS(H)-GAUSS(K)+1.0D+00
BIVARIATE=BV
RETURN
END
FUNCTION BV1(H1,H2,R)
IMPLICIT NONE
INTEGER I
! IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 H1,H2,R,H12,BV,BV1,RR,RR2,H3,GAUSS
EXTERNAL GAUSS
REAL*8 X(5),W(5),DVAL
DATA X/.04691008,.23076534,.5,.76923466,.95308992/
DATA W/.018854042,.038088059,.0452707394,.038088059,.018854042/
H12 =(H1*H1+ H2*H2)/2.0D+00
H3 = H1*H2
BV = 0.0D+00
DO 1 I= 1,5
RR = R*X(I)
RR2= 1.0D+00-RR*RR
IF(RR2.LE.0.0D+00) RR2=1.0D-100
DVAL=(RR*H3 - H12)/RR2
IF(DVAL.GT.300.0D+00) DVAL=300.0D+00
BV = BV + W(I)*EXP(DVAL)/SQRT(RR2)
1 CONTINUE
BV = BV*R + GAUSS(H1)*GAUSS(H2)
BV1=BV
RETURN
END
FUNCTION BV2(H1,HK,R)
IMPLICIT NONE
! IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 X(5),W(5)
REAL*8 H1,HK,R,H2,H12,BV,BV2,R2,R3,H3,H7,H6,H5,AA,AB,R1,RR,GAUSS,RR2
REAL*8 DVAL,DVAL2
EXTERNAL GAUSS
INTEGER I
DATA X/.04691008,.23076534,.5,.76923466,.95308992/
DATA W/.018854042,.038088059,.0452707394,.038088059,.018854042/
H2 = HK
H12 = (H1*H1+ H2*H2)/2.0D+00
BV = 0.0D+00
IF(ABS(R).LT.0.7D+00)GO TO 4
R2= 1.0D+00-R*R
R3 = SQRT(R2)
IF(R.LT.0.0D+00)H2 = -H2
H3=H1*H2
H7 = EXP( -H3/2.0D+00)
IF(R2.EQ.0.0D+00) GO TO 3
H6 = ABS(H1 - H2)
H5 = H6*H6/2.0D+00
H6 = H6/R3
AA = 0.5D+00 - H3/8.0D+00
AB = 3.0D+00 - 2.0D+00*AA*H5
BV = .13298076D+00*H6*AB*GAUSS(H6) &
-EXP(-H5/R2)*(AB + AA*R2)*.053051647D+00
DO 2 I= 1,5
R1= R3*X(I)
RR = R1*R1
R2 = SQRT(DABS(1.0D+00 - RR))
IF(R2.EQ.0.0D+00) R2=1.0D-100
IF(RR.EQ.0.0D+00) RR=1.0D-100
DVAL=-H5/RR
IF(DVAL.GT.300.0D+00) DVAL=300.0D+00
DVAL2=-H3/(1.0D+00 + R2)
IF(DVAL2.GT.300.0D+00) DVAL2=300.0D+00
BV = BV - W(I)*EXP(DVAL)*(EXP(DVAL2)/R2/H7 - 1.0D+00 - AA*RR)
2 CONTINUE
3 IF(R.GT.0.0D+00)BV = BV*R3*H7 + GAUSS(MAX(H1,H2))
IF(R.LT.0.0D+00)BV = MAX(0.0D+00,GAUSS(H1) - GAUSS(H2)) - BV*R3*H7
BV2=BV
RETURN
4 H3=H1*H2
DO 1 I= 1,5
R1= R*X(I)
RR2= 1.0D+00-R1*R1
IF(RR2.EQ.0.0D+00) RR2=1.0D-100
1 BV = BV + W(I)*EXP((R1*H3 -H12)/RR2)/SQRT(RR2)
BV2 =GAUSS(H1)*GAUSS(H2) + R*BV
RETURN
END
FUNCTION GAUSS(Z)
IMPLICIT NONE
! IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 Z,GAUSS,X,G
REAL*8 A(4)
INTEGER I
DATA A/ -.72657601,.71070688,-.142248368,.127414796/
X= 1.0D+00/(1.0D+00 + .23164189D+00*ABS(Z))
G = .53070271D+00
DO 1 I= 1,4
1 G=G*X+A(I)
GAUSS = G*X*EXP( - Z*Z/2.0D+00)
IF(Z.LT.0.0D+00)GAUSS = 1.0D+00 -GAUSS
RETURN
END
! USE PROGRAM HEADER FOR STAND-ALONE EXECUTABLE TESTING
! PROGRAM BIVTEST
SUBROUTINE BIVTEST
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 A(9),A1(9),A2(9,9)
REAL*8 AH(9),AH1(9),AH2(9,9)
REAL*8 AJ(9),AJ1(9),AJ2(9,9)
REAL*8 B(9),B1(9),B2(9,9)
REAL*8 H,K
HDEL=1.0D-04
1 CONTINUE
WRITE(*,*) 'ENTER RHO,H,K,INTEGRAND-TYPE,BVTYPE'
READ(*,*) RHO,H,K,DTYPE,BVTYPE
A(1)=RHO
A(2)=H
A(3)=K
A(4)=DTYPE
A(5)=BVTYPE
BV=FUNCA(A,A1,A2)
WRITE(*,*) 'VALUE ',BV
WRITE(*,*) GAUSS(-1.96D+00),GAUSS(1.96D+00)
AH(1:5)=A(1:5)
AH(1)=RHO-HDEL
AJ(1:5)=A(1:5)
AJ(1)=RHO+HDEL
BVH=FUNCA(AH,AH1,AH2)
BVJ=FUNCA(AJ,AJ1,AJ2)
B1(1)=(BVJ-BVH)/2.0D+00/HDEL
B2(1,1)=(AJ1(1)-AH1(1))/2.0D+00/HDEL
B2(1,2)=(AJ1(2)-AH1(2))/2.0D+00/HDEL
B2(1,3)=(AJ1(3)-AH1(3))/2.0D+00/HDEL
WRITE(*,*)
WRITE(*,*) 'ANALOG RHO',A1(1),A2(1,1),A2(1,2),A2(1,3)
WRITE(*,*) 'NUMER RHO',B1(1),B2(1,1),B2(1,2),B2(1,3)
AH(1:5)=A(1:5)
AH(2)=H-HDEL
AJ(1:5)=A(1:5)
AJ(2)=H+HDEL
BVH=FUNCA(AH,AH1,AH2)
BVJ=FUNCA(AJ,AJ1,AJ2)
B1(2)=(BVJ-BVH)/2.0D+00/HDEL
B2(2,1)=(AJ1(1)-AH1(1))/2.0D+00/HDEL
B2(2,2)=(AJ1(2)-AH1(2))/2.0D+00/HDEL
B2(2,3)=(AJ1(3)-AH1(3))/2.0D+00/HDEL
WRITE(*,*)
WRITE(*,*) 'ANALOG H',A1(2),A2(2,1),A2(2,2),A2(2,3)
WRITE(*,*) 'NUMER H',B1(2),B2(2,1),B2(2,2),B2(2,3)
AH(1:5)=A(1:5)
AH(3)=K-HDEL
AJ(1:5)=A(1:5)
AJ(3)=K+HDEL
BVH=FUNCA(AH,AH1,AH2)
BVJ=FUNCA(AJ,AJ1,AJ2)
B1(3)=(BVJ-BVH)/2.0D+00/HDEL
B2(3,1)=(AJ1(1)-AH1(1))/2.0D+00/HDEL
B2(3,2)=(AJ1(2)-AH1(2))/2.0D+00/HDEL
B2(3,3)=(AJ1(3)-AH1(3))/2.0D+00/HDEL
WRITE(*,*)
WRITE(*,*) 'ANALOG K',A1(3),A2(3,1),A2(3,2),A2(3,3)
WRITE(*,*) 'NUMER K',B1(3),B2(3,1),B2(3,2),B2(3,3)
GO TO 1
ENDREFERENCES: Guide Introduction_7
19 CCONTR EXAMPLE
This example illustates the use of the CONTR and CCONTR stubs with NONMEM 7. This example was developed by Bill Frame.
$PROBLEM
$DATA NMDATA7.CSV
$INPUT ID DV
$SUB CONTR=CONTR.txt CCONTR=CCONTR.txt
$PRED
W=THETA(2) ;SD
CL=THETA(1)*EXP(ETA(1)) ;CL/F
PRE=1/CL ; SS ASSUMING INPUT RATE = 1
LAM=THETA(3)
Y=(PRE**LAM-1)/LAM+EPS(1)*W
RES1=(DV-PRE)/W
$THETA
(0,.1) ;CL/F
(0,2) ;SD ADDITIVE
(0,.251) ;BOX COX LAMBDA PARAMETER
$OMEGA
.003
$SIGMA
1 FIX
$EST MAXEVALS=9999 METH=1 PRINT=1
$COV PRINT=E
Files CONTR.txt and CCONTR.txt are as follows:
SUBROUTINE CONTR (ICALL,CNT,IER1,IER2)
DOUBLE PRECISION CNT
CALL NCONTR (CNT,IER1,IER2,L2R)
RETURN
END
SUBROUTINE CCONTR(I,CNT,P1,P2,IER1,IER2)
USE SIZES, ONLY: ISIZE,DPSIZE
USE ROCM_REAL, ONLY: THETA=>THETAC,Y=>DV_ITM2
USE NM_INTERFACE,ONLY: CELS
IMPLICIT NONE
INTEGER(KIND=ISIZE), INTENT(IN OUT) :: I,IER1,IER2
REAL(KIND=DPSIZE), INTENT(IN OUT) :: CNT,P1(:),P2(:,:)
REAL(KIND=DPSIZE) :: ONE,TWO,W
DATA ONE,TWO/1.00D+00,2.00D+00/
SAVE
IF (I.LE.1) RETURN
W=Y(1)
Y(1)=(Y(1)**THETA(3)-ONE)/THETA(3)
CALL CELS (CNT,P1,P2,IER1,IER2)
Y(1)=W
CNT=CNT-TWO*(THETA(3)-ONE)*LOG(Y(1))
RETURN
END
The first individual's data follows.
id,dv
1,223.28
1,41.59
1,87.6
1,72.11
1,216.32
1,118.1
1,232.29REFERENCES: Guide Introduction_7
20 EXAMPLES USING MTIME TO MODEL PERIODIC DISCONTINUITIES IN $DES
Some systems defined with differential equations can be discontinuous with respect to time. Discontinuities are typically introduced in the system by suddenly changing the value of one or more model variables at specific points of time. Such changes can be periodic due to e.g. a circadian rhythm. It is important to define discontinuous variables in
$PK using model event time (MTIME) variables rather than updating them
in $DES using IF-ELSE-ENDIF tests of T or using the discontinuous INT
function. In fact, when integrating from time ta to time tb, rou-
tines such as DVERK (ADVAN6) may go slightly beyond time tb (i.e., it
may happen that $DES is called at T>tb). The values of DADT should be
continuous. If any element of DADT changes at T=b, this should happen
at the next integration interval, integrating from time b to time c.
Two examples are provided for the implementation of periodic disconti-
nuities using MTIME variables. The first (step_circexa.ctl) shows
how to model the daily reset of a step function and illustrates how
NONMEM update the MTIME and closely related MNOW variables at and in
between event records. The second (idr_circexa.ctl) applies this step
function in an indirect response model that describes a truncated
sinusoidal baseline reponse. It also shows how values of T that
range from 0 to infinity can be transformed into repeated values of
0-24h. All files may be found in the examples directory. These exam-
ples were suggested by Sebastien Bihorel, Luann Phillips and Jill
Fiedler-Kelly.
The data set used in both examples is circadian.csv. This data set has
4 subjects with sparse data (one observation every 24 hours from time
0 to time 168) and 4 subjects with finely-spaced data (every 0.5 hours
from time 0 to time 48). It is with the second group of subjects that
scatters of the step function, variable T24, and indirect response vs.
time can best be seen.
There are no dose event records (i.e. AMT=0 for all records) because
the step function and indirect response are generated endogenously.
If this model is incorporated in a larger model, dose records can be
used for other compartments.
step_circexa.ctl
This example illustrates how to turn on a step function, FLAG, at
a certain time of the day, SHIFT, turn it off after a certain
duration, DUR, and repeat this every 24 hours.
^
1 -| --------- --------- ---------
| | | | | |
| | | | | |
0 -|--- -------- -------- -------- FLAG
---------------------------------------------------> TIME
| | |
SHIFT-24 0 24 48
<- SHIFT -><- DUR ->
The MTIME(1) and MTIME(2) variables respectively define the times
at which FLAG is set to 1 and 0, such as MTIME(2)-MTIME(1)=DUR.
After MTIME(2) is reached, both variables are incremented by 24,
and the FLAG update process perpetuates itself every 24 hours.
Compartment 1 is a dummy compartment intended to monitor the pro-
per update of the FLAG variable: the amount A1 starts at 0 and is
incremented by DUR every day.
Because the FLAG variable is not necessarily updated at the time
of event records, both $TABLE and WRITE statements within the ab-
breviated code are used in step_circexa.ctl to report the values
of the variables of interest. WRITE statements are used because
lines in a table file are produced only for event records, but
not for nonevent times such as Model times (MTIMES) at which the
variables computed in $PK are of special interest here.
Besides MTIME(1) and MTIME(2), the following variables are of
interest in step_circexa.ctl:
* MPAST(1) and MPAST(2): variables automatically updated based on
MTIME(1) and MTIME(2), and used to set FLAG in $DES. MPAST(i)=0
until the call to PK subsequent to the one for which MNOW=i. At
that call and until MTIME(i) is redefined, MPAST(i)=1.
* INTMTIME ("interval MTIME") is the value of MTIME(1) for the
entire integration interval during which FLAG=1, ie from MTIME(1)
up to and including the endpoint at MTIME(2). It is not specifi-
cally used in step_circexa.ctl, but will be in the next example.
INTMTIME is computed and displayed here so that it can be dis-
cussed with the other variables.
* TSTATE is a reserved variable giving the time at which the
current state vector (compartment amounts) was computed. It gives
the time to which the system was most recently advanced.
The FLAG variable can be computed in either $PK or $DES, but the
latter is preferred.
Background details:
Suppose that the values of TIME on the event records are t1 t2 t3
etc. and that PK is called with every event record and NONEVENT
(ADDITIONAL AND LAGGED) DOSE TIMES AND AT MODEL TIMES.
Suppose there are k MTIME variables with values
t1 <= MTIME(1) <= MTIME(2) <= ... <= MTIME(k) <= t2
The interval [t1,t2] is integrated by smaller intervals
[t1, MTIME(1)], [MTIME(1), MTIME(2)] ... [MTIME(k), t2].
Calls to DES during any integration interval [ta, tb] will have
ta <= T, but it may happen that T > tb. For calls to PK and
ERROR, the sequence of calls for the record with TIME=t2 is:
* for each i from 1 to k: call PK with the record at MTIME(i)
(MNOW=i, TSTATE=MTIME(i-1)) (if i=1, TSTATE=t1)
* call PK with the record at t2
(MNOW=0; TSTATE=MTIME(k))
* call ERROR with the record at t2
(MNOW is set to 9 to identify the WRITE in $ERROR; TSTATE=t2)
The same calls to PK occur even if MTIME(k) happens to corres-
pond to t2. In such a case, the call with MNOW=k precedes the
call with MNOW=0, and the values in the table file are from the
second call.
The lines of the file step_circexa.txt that corresponding to the
second event record (TIME=24) for ID 1 (set with SHIFT=0.9088 and
DUR=20.8385), are:
MNOW TIME MT1 MT2 MP1 MP2 TSTATE INTMTIME
(1) 1.0000 24.0000 0.9088 21.7473 0.0000 0.0000 0.0000 0.9088
(2) 2.0000 24.0000 24.9088 45.7473 1.0000 0.0000 0.9088 0.9088
(3) 0.0000 24.0000 24.9088 45.7473 0.0000 0.0000 21.7473 24.9088
(4) 9.0000 24.0000 24.9088 45.7473 0.0000 0.0000 24.0000 24.9088
where MTi and MPi stand for MTIME(i) and MPAST(i).
At line 1, FLAG=0. This corresponds to the advance from time 0 to
MTIME(1): no update of MTIME variables is performed.
At line 2, FLAG=1. This corresponds to the advance from MTIME(1)
to MTIME(2), at which point the MTIME(1) and MTIME(2) variables
are updated, but INTMTIME retains its value.
At line 3, FLAG=0. This corresponds to the advance from MTIME(2)
to time 24, i.e. the call to PK with MNOW=0, in which INMTIME is
updated to the new value of MTIME(1) defined at the previous call
to PK.
At line 4, FLAG=0. This corresponds to the call to ERROR.
The corresponding line of the table file step_circexa_debug.tab
is:
TIME MT1 MT2 STS INTMTIME
24.0000 24.9088 45.7473 21.7473 24.9088
Note that TSTATE values are:
0 (time of 1st. event record)
0.9088 after advance to MTIME(1)=0.9088
21.7473 after advance to MTIME(2)=21.7473
24.0000 after advance to t2=24.
There will be as many lines with MNOW=0.0 in step_circexa.txt, as
there are lines in step_circexa_debug.tab. These lines will be
consistent except for the values of MNOW.
INMTIME retains the value of MTIME(1)=0.9088 that pertains to the
entire integration interval, up to and including the end point at
MTIME(2). It is set in $PK before the MTIME values are changed
when MNOW=2. The values in the table are those set in $PK when
MNOW=0 (the final call to PK with this data record) by which time
INTMTIME and the MTIME(i) have their new values for the next set
of advances.
In this particular example, with simulated etas, it sometimes
happens that MTIME(1)<t1. From a biological point of view, the
step function (and the response) likely started well before the
data collection. The beginning of the data collection may occur
with FLAG=1 (in the next example, this is during the oscillation
of the response) or with FLAG=0 (during the flat part of the
response). When MTIME(1) < t1, MTIME(1) is ignored during the
set of calls to PK described above.
idr_circexa.ctl
This example extends the concepts introduced in step_circexa.ctl
and illustrates how to apply the daily reset of the step function
to model a process whose baseline is characterized by a sinusoid
function during part of the day and a flat line for the rest of
the day. It also show how to create a 24h clock time variable T24
in $DES by transforming the T variable ranging from 0 to infinity
into repeating intervals of continuous time between 0 and 24.
^
Rmin+Amp -| --- --- ---
| | | | | |
| | | | | |
| | | | | |
Rmin -|-- -------- -------- ---- Response R
-------------------------------------------------> TIME
| | |
SHIFT-24 0 24 48
<- SHIFT -><- DUR ->
The baseline response depicted above can be described by:
R(t) = | Rmin + Amp*sin(z(t)) for t in [SHIFT+n*24,SHIFT+DUR+n*24] (Eq.1)
| Rmin for any other value of t
where
Rmin is the minimum response
Amp is the amplitude of the sinusoid function of time
z(t) is a circadian function of time which scales the time
intervals [SHIFT+n*24, SHIFT+DUR+n*24] to [0, PI], that is
an interval of x such that sin(x)>=0:
z(t) = PI * (t-(SHIFT+n*24))/DUR
n is (-1,0,1,2,...) if (SHIFT+DUR)>24
is ( 0,1,2,...) if (SHIFT+DUR)<=24
The case when SHIFT+DUR>24 is when the flag function is already 1
at time 0, that is it started before the first event. Therefore,
an offset of 24 is needed.
z(t) is the function that drives the circadian rhythm of R and is
implemented using the INTMTIME variable described in the previous
example and that represents SHIFT+n*24. MTIME(1) and MTIME(2) are
defined as in step_circexa.ctl to periodically update INTMTIME
and the FLAG variable used later in the $DES block.
The baseline function, R(t), can be fitted using Eq.1 in absence
of drug effect. However, some biological processes which follow
such circadian rhythms can be influenced by drugs (eg, metabolic
effects of corticoids). If the effect of the drug is not direct,
one can express the biological response and the drug effect in
terms of an indirect response model (IDR). Assuming that the drug
stimulates the formation of the biological response, one can
parameterize the system using the following differential equa-
tion:
dR/dt = KIN(t)*(1+S(Cp)) - KOUT*R (Eq.2)
where
KOUT is the elimination rate of the biological response and
is stationary
KIN(t) is the formation rate of the biological response and is
time-varying dependent on the step function FLAG:
FLAG = MPAST(1)-MPAST(2)
S(Cp) is a stimulation function of the drug concentration.
In absence of drug, there is no stimulation (i.e., S(Cp) = 0) and
the baseline response R(t) described in Eq.1 can be reparameter-
ized as follows:
dR/dt = KIN(t) - KOUT*R (Eq.3)
The model is defined such as, when FLAG=0, KIN(t) is constant and
equal to KOUT*Rmin, and such as, when FLAG=1, KIN(t) has a sinu-
soidal shape, increasing from its minimum value to its maximum,
and decreasing back to its minimum within the duration DUR. One
can obtain the explicit expression of KIN(t) within the time when
FLAG=1, by determining the first derivative of R(t), as defined
in the first part of Eq.1 and solving for KIN(t) using Eq.3. This
expression of KIN(t) is used in $DES.
The $DES block also includes an independent line of code meant to
define a continuous 24h clock time. This code is supported by the
definition and circadian reset of the MTIME(3) variable in $PK.
More specifically, this code is dependent on INTMTIME3, the value
of MTIME(3) for the 24-hour intervals.
Note: Original idea for this example stems from Andreas Krause's
post on the NMUsers email list dated from 5/26/2011, subject
"Coding INTEGER Function in NONMEM". The original code was modi-
fied to use model event time variables.
(See mtime).REFERENCES: None
21 CONTROL3 EXAMPLE
This is file CONTROL3 from the NONMEM distribution medium.
$PROBLEM THEOPHYLLINE SINGLE SUBJECT DATA
$INPUT DOSE=AMT TIME CP=DV
$DATA DATA3
$SUBROUTINES ADVAN2
$PK
CALLFL=1
KA=THETA(1)
K=THETA(2)
SC=THETA(3)
$ERROR
Y=F+ERR(1)
$THETA (0,1.7) (0,.102) (0,29)
$ESTIMATION MAXEVAL=240 PRINT=2
$COVR
$TABLE TIME
$SCAT CP VS TIME
$SCAT PRED VS TIME
$SCAT RES VS TIME
$SCAT PRED VS CP UNIT
The data is file DATA3 from the NONMEM distribution media.
320 .0 .
. .27 1.71
. .52 7.91
. 1. 8.31
. 1.92 8.33
. 3.5 6.85
. 5.02 6.08
. 7.03 5.4
. 9. 4.55
. 12. 3.01
. 24.3 .90REFERENCES: Guide III, section VI.1
22 CONTROL4 EXAMPLE
This is file CONTROL4 from NONMEM distribution media
$PROB THEOPHYLLINE POPULATION DATA
$INPUT ID DOSE TIME CP=DV WT
$DATA THEO
$PRED
;THETA(1)=MEAN ABSORPTION RATE CONSTANT (1/HR)
;THETA(2)=MEAN ELIMINATION RATE CONSTANT (1/HR)
;THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
;DOSE=WT-ADJUSTED DOSE (MG/KG)
;DS=NON-WT-ADJUSTED DOSE (MG)
IF (DOSE.NE.0) THEN
DS=DOSE*WT
W=WT
ENDIF
KA=THETA(1)+ETA(1)
KE=THETA(2)+ETA(2)
CL=THETA(3)*W+ETA(3)
D=EXP(-KE*TIME)-EXP(-KA*TIME)
E=CL*(KA-KE)
F=DS*KE*KA/E*D
Y=F+EPS(1)
$THETA (.1,3,5) (.008,.08,.5) (.004,.04,.9)
$OMEGA BLOCK(3) 6 .005 .0002 .3 .006 .4
$SIGMA .4
$EST MAXEVAL=450 PRINT=5
$COV
$TABLE ID DOSE WT TIME
$SCAT (RES WRES) VS TIME BY ID
The data is file THEO from the NONMEM distribution media. The first
individual's data follows.
1 4.02 0. .74 79.6
1 . 0.25 2.84 .
1 . 0.57 6.57 .
1 . 1.12 10.5 .
1 . 2.02 9.66 .
1 . 3.82 8.58 .
1 . 5.1 8.36 .
1 . 7.03 7.47 .
1 . 9.05 6.89 .
1 . 12.12 5.94 .
1 . 24.37 3.28 .REFERENCES: Guide III, section VI.3; Guide IV, Appendix II
23 CONTROL5 EXAMPLE
This is file CONTROL5 from NONMEM distribution media
$PROB THEOPHYLLINE POPULATION DATA
$INPUT ID DOSE=AMT TIME CP=DV WT
$DATA THEOPP
$SUBROUTINES ADVAN2
$PK
;THETA(1)=MEAN ABSORPTION RATE CONSTANT (1/HR)
;THETA(2)=MEAN ELIMINATION RATE CONSTANT (1/HR)
;THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
;SCALING PARAMETER=VOLUME/WT SINCE DOSE IS WEIGHT-ADJUSTED
CALLFL=1
KA=THETA(1)+ETA(1)
K=THETA(2)+ETA(2)
CL=THETA(3)*WT+ETA(3)
SC=CL/K/WT
$THETA (.1,3,5) (.008,.08,.5) (.004,.04,.9)
$OMEGA BLOCK(3) 6 .005 .0002 .3 .006 .4
$ERROR
Y=F+EPS(1)
$SIGMA .4
$EST MAXEVAL=450 PRINT=5
$COV
$TABLE ID DOSE WT TIME
$SCAT (RES WRES) VS TIME BY ID
The data is file THEOPP from the NONMEM distribution media. The first
individual's data follows.
1 4.02 0. . 79.6
1 . 0. .74 .
1 . 0.25 2.84 .
1 . 0.57 6.57 .
1 . 1.12 10.5 .
1 . 2.02 9.66 .
1 . 3.82 8.58 .
1 . 5.1 8.36 .
1 . 7.03 7.47 .
1 . 9.05 6.89 .
1 . 12.12 5.94 .
1 . 24.37 3.28 .24 CONTROL6 EXAMPLE
This is file CONTROL6 from NONMEM distribution media. It is discussed in Guide IV, Appendix VIII.
$PROB THEOPHYLLINE POPULATION DATA
$INPUT ID DOSE=AMT TIME CP=DV WT
$DATA THEOPP
$SUBROUTINES ADVAN7
$MODEL COMP=(DEPOT,INITIALOFF,DEFDOSE) COMP=(CENTRAL,DEFOBS,NOOFF)
$PK
;THETA(1)=MEAN ABSORPTION RATE CONSTANT (1/HR)
;THETA(2)=MEAN ELIMINATION RATE CONSTANT (1/HR)
;THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
;SCALING PARAMETER=VOLUME/WT SINCE DOSE IS WEIGHT-ADJUSTED
CALLFL=1
K12=THETA(1)+ETA(1)
K20=THETA(2)+ETA(2)
CL=THETA(3)*WT+ETA(3)
S2=CL/K20/WT
$THETA (.1,3,5) (.008,.08,.5) (.004,.04,.9)
$OMEGA BLOCK(3) 6 .005 .0002 .3 .006 .4
$ERROR
Y=F+EPS(1)
$SIGMA .4
$EST MAXEVAL=450 PRINT=5
$COV
$TABLE ID DOSE WT TIME
$SCAT (RES WRES) VS TIME BY ID
The data is file THEOPP from the NONMEM distribution media. The first
individual's data follows.
1 4.02 0. . 79.6
1 . 0. .74 .
1 . 0.25 2.84 .
1 . 0.57 6.57 .
1 . 1.12 10.5 .
1 . 2.02 9.66 .
1 . 3.82 8.58 .
1 . 5.1 8.36 .
1 . 7.03 7.47 .
1 . 9.05 6.89 .
1 . 12.12 5.94 .
1 . 24.37 3.28 .25 CONTROL7 EXAMPLE
This is file CONTROL7 from NONMEM distribution media It is discussed in Guide IV, Appendix VIII.
$PROB THEOPHYLLINE POPULATION DATA
$INPUT ID DOSE=AMT TIME CP=DV WT
$DATA THEOPP
$SUBROUTINES ADVAN6 TOL=5
$MODEL COMP=(DEPOT,INITIALOFF,DEFDOSE) COMP=(CENTRAL,DEFOBS,NOOFF)
$PK
;THETA(1)=MEAN ABSORPTION RATE CONSTANT (1/HR)
;THETA(2)=MEAN ELIMINATION RATE CONSTANT (1/HR)
;THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
;SCALING PARAMETER=VOLUME/WT SINCE DOSE IS WEIGHT-ADJUSTED
CALLFL=1
KA=THETA(1)+ETA(1)
KE=THETA(2)+ETA(2)
CL=THETA(3)*WT+ETA(3)
S2=CL/KE/WT
$THETA (.1,3,5) (.008,.08,.5) (.004,.04,.9)
$OMEGA BLOCK(3) 6 .005 .0002 .3 .006 .4
$DES
DADT(1)=-KA*A(1)
DADT(2)= KA*A(1)-KE*A(2)
$ERROR
Y=F+EPS(1)
$SIGMA .4
$EST MAXEVAL=450 PRINT=5
$COV
$TABLE ID DOSE WT TIME
$SCAT (RES WRES) VS TIME BY ID
The data is file THEOPP from the NONMEM distribution media. The first
individual's data follows.
1 4.02 0. . 79.6
1 . 0. .74 .
1 . 0.25 2.84 .
1 . 0.57 6.57 .
1 . 1.12 10.5 .
1 . 2.02 9.66 .
1 . 3.82 8.58 .
1 . 5.1 8.36 .
1 . 7.03 7.47 .
1 . 9.05 6.89 .
1 . 12.12 5.94 .
1 . 24.37 3.28 .REFERENCES: Guide III, section VI.3; Guide IV, Appendix VI , VIII
26 DATA_AVERAGE EXAMPLE
A template record is a data record that contains a non-zero value for the raw-data (RAW_) data item. Only the DV values with those observation records (i.e., records with MDV=0) matching the template data record will be included in the average. This average replaces the DV item in the template record, and it may be displayed in a table or scatterplot. An observation record matches a template if all its data items are identical to all those on the template, except that the values of NONMEM data items do not affect the match. (Any nonNONMEM data item type may in fact be excluded from the match (See omit)). Example 1 Here is a fragment of a data set (which may or may not actually be for PREDPP): #ID TIME AMT WT AGE EVID MDV DV RAW_ GEND 1 5 0 70 50 0 1 2 1 1 template record 4 5 0 70 50 0 0 4 0 1 observation record These two records match. The NONMEM items are the ID, MDV, DV, and RAW_ items. The DV item with the observation record will be included in the average. The DV item in the template record is ignored for the purpose of computing the average (MDV=1); the displayed DV item with this record will be the average value, and the displayed RES item will be the difference between this average and the PRED item. The average is taken over all records matching the TIME, AMT, WT, AGE, EVID, and GEND items of the template record. If the omit record
$OMIT GEND
appears, then the two records
#ID TIME AMT WT AGE EVID MDV DV RAW_ GEND
1 5 0 70 50 0 1 2 1 1 template record
4 5 0 70 50 0 0 4 0 2 observation record
also match. In this case the average is taken over all records match-
ing the TIME, AMT, WT, AGE, and EVID items of the template record.
Example 2
Here is an example of categorical dose-response with four different
bolus doses per individual (separated by washouts), and data averages
are formed for each of the four doses.
#ID DOSE DV MDV RAW_
100 5 . 1 1 set DV item to avr. for DOSE=5
100 10 . 1 1 set DV item to avr. for DOSE=10
100 20 . 1 1 set DV item to avr. for DOSE=20
100 40 . 1 1 set DV item to avr. for DOSE=40
1 5 1 0 0
1 10 1 0 0
1 20 2 0 0
1 40 2 0 0
2 5 1 0 0
2 10 1 0 0
2 20 1 0 0
2 40 1 0 0
3 5 1 0 0
3 10 2 0 0
3 20 2 0 0
3 40 3 0 0
.
.
Example 3
An average of a PRED-defined item may also be displayed.
$ABBREV COMRES=1 COMSAV=1
...
$PRED
COM(1)=LOG(DV)
...
$TABLE COM(1)=LDV
For template records, the value of LDV displayed in the table is the
average of LOG DV over all observation records matching the template.
(See data average block, raw).REFERENCES: none
27 INTERPOLATION IN $DES
This is an example of code that can be used in $DES to interpolate a covariate such as WT between values that are recorded on the data records. This example is new to NONMEM 7.4, although the code can be used with previous versions. If D_WT is used in the model for volume of distribution, then this would allow a time-varying volume of distribution to be integrated. All the code in $PK to compute OLDTIME and OLDWT and SLOPE, and the code for D_WT in $DES, could be copied to user’s control stream. Additional code (for integrating D_WT in $DES and the geometric solution in $ERROR) is for testing and would not be part of the user’s control stream. The code cannot be used when any WT values are missing. $PK sees only the previous and current record. If WT is missing from the current record, $PK has no way of knowing what WT will be on the next record. The code in NONMEM 7.3 help file INFN_INTERPOLATION EXAMPLE 2 (infn2.exa) can be used to fill in missing WT values in the event records. The $INFN block makes two passes through the data set. After the first pass, it knows what the WT values are on all the event records, so that it can interpolate to fill in a missing WT. The finedata utility (fine1.exa) can also be used to fill in the missing values, and performs interpolation if necessary. (See infn2.exa). (See fine1.exa). Here is the control file:
$PROB INTERPOLATE WT IN $DES
; this example shows how to interpolate WT in $DES.
; it is assumed that WT is recorded on every data record.
; As a test, the value of D_WT in $DES is integrated to
; obtain AUC of WT VS. T
; This is also calculated geometrically in $ERROR.
$INPUT ID TIME WT DV
$DATA desinterp.dat
$SUBROUTINES ADVAN6 TOL=5
$MODEL
COMP=(AUC_WT DEFOBS)
$PK
; initialize OLDTIME and OLDWT
IF (NEWIND.LE.1) THEN
OLDTIME=TIME
OLDWT=WT
ENDIF
; calculate the slope for $DES
DELTA_TIME=TIME-OLDTIME
DELTA_WT=WT-OLDWT
IF (DELTA_TIME>0) THEN
SLOPE=DELTA_WT/DELTA_TIME
ELSE
SLOPE=0.
ENDIF
; save wt and time for next $PK record
OLDTIME=TIME
OLDWT=WT
$DES
D_WT=OLDWT+SLOPE*(T-OLDTIME) ; D_WT is the value of WT at time T
DADT(1)=D_WT ; compute AUC of D_WT as a test
$ERROR
Y=F+ETA(1)+EPS(1)
; Compute geometric solution as a test.
; Does not use compartment amounts.
; Use only the values of WT and TIME on event records.
; Suppose WT vs T looks like this:
;
;
; WT
;
; | w3 w4
; | w2 w5
; | w1
;
; --------------------------------> TIME
; t1 t2 t3 t4 t5
;
; at t2, the contribution to the sum is
; the rectangle w1 x (t2-t1)
; plus the triangular piece
; (w2-w1)/(t2-t1) / 2
;
; w2
; /|
; /
; /__|
; w1
; |
; |
; ------------
; t1 t2
IF (NEWIND.LE.1) THEN
PREV_WT=WT ; Initialize WT from previous data record
SUM=0
ELSE
SUM=SUM+PREV_WT*DELTA_TIME+DELTA_WT*DELTA_TIME/2
ENDIF
PREV_WT=WT ; save WT from previous data record
$OMEGA 1
$SIGMA 1
$TABLE ID TIME WT PRED=AUC_WT SUM FILE=desinterp.tbl NOPRINT NOAPPEND
; The following two values should always be equal:
; PRED (which is the AUC of WT obtained by integrating WT)
; SUM (which is the geometric solution) computed in $ERROR
Here is the data file for the first subject. Note that WT sometimes is
constant and sometimes decreases:
1 0. 10 0
1 1. 20 0
1 2. 35 0
1 2. 35 0
1 4. 45 0
1 5. 40 0
Here is the table file:
TABLE NO. 1
ID TIME WT AUC_WT SUM
1.0000E+00 0.0000E+00 1.0000E+01 0.0000E+00 0.0000E+00
1.0000E+00 1.0000E+00 2.0000E+01 1.5000E+01 1.5000E+01
1.0000E+00 2.0000E+00 3.5000E+01 4.2500E+01 4.2500E+01
1.0000E+00 2.0000E+00 3.5000E+01 4.2500E+01 4.2500E+01
1.0000E+00 4.0000E+00 4.5000E+01 1.2250E+02 1.2250E+02
1.0000E+00 5.0000E+00 4.0000E+01 1.6500E+02 1.6500E+0228 ENTERHEPATIC CIRCULATION EXAMPLES
The help item “Model Time examples” gives a fragment of code for modelling EHC (Enterohepatic Circulation, which is also called Enterohepatic Recycling). It illustrates the use of MTIME parameters to model instantaneous changes in differential equations. This help item, “Enterhepatic circulation examples”, describes two fully-worked out control streams in the examples directory. They demonstrate how to generalize the fragment to multiple sequential doses. Both examples use the same data. (hillss.dat and mtimess.dat are identical). There is a steady-state bolus dose at time 0. The inter-dose interval II is 12, and there are enough additional doses (ADDL=100) to continue the dosing pattern throughout the data set. There are “other” records every 4 units till time 140 to allow compartment amounts to be displayed and there is a final observation record at time 144. mtimess.ctl This example incorporates the fragment into a complete control stream. MTIME parameters are used to turn on and off the EHC terms in the differential equations. The variable FLAG is 1 between times MTIME(1) and MTIME(2) after each dose event and turns on the EHC terms. After time MTIME(2) is reached, a new set of MTIME’s is defined which affect the next dosing interval. MTIME parameters are not dose-related parameters and have no effect on steady-state dose events. Even if PK computes MTIME(i)< II, this produces future changes in the system, and does not apply retroactively to the preceding implied doses. (See Guide VI, Section V.F.4, Note 4). A steady-state dose record should not be used. Instead, the SS dose record is described as a transient dose with SS=DROP on the
$INPUT record.
hillss.ctl
This control stream does not use MTIME. Instead, a smooth step
model using Hill terms in a sigmoid emax model is used. The $DES
code has to compute all the necessary variables. Flag1 and flag2
are continuous variables that change from 0 to 1 at the times
corresponding to the MTIME's. The FLAG variable is similar to
FLAG in mtimess.ctl. The changes to the differential equations
are not instantaneous, but they are continuous. If the exponent
in the Hill term is made larger, the predictions approach those
of the MTIME model. However, very large values of the exponent
can lead to numerical difficulties in PREDPP. Smaller values of
the exponent may be more realistic physiologically. A Steady-
State dose event record is used with this model.
Note that SS dose records should only be used when the kinetics
implemented in the model coincides with the II (interdose inter-
val) of the SS record. Just as MTIME's cannot affect the differ-
ential equations retroactively, changes to the differential equa-
tions that happen in the future cannot affect the Steady-state
calculations. For example, with II=12, then the kinetics should
not be different in the interval of time 0 to 12 vs. time 12 to
24 or time 24 to 36, etc. In each interval, changes occur at
theta(8) and theta(8)+theta(9) after the start of the interval.
$DES computes the first change time for each interval using the
INT function:
mt1=II*INT(T/II)+theta(8)
NM-TRAN gives a warning about the use of the INT function in
$DES:
(WARNING 68) THE INT, MOD, MIN, OR MAX FUNCTION IS BEING USED OUT-
SIDE OF A
SIMULATION BLOCK. IF THE FUNCTION VALUE AFFECTS THE VALUE OF THE
OBJECTIVE FUNCTION, THEN AN ERROR WILL PROBABLY OCCUR.
This warning may be disregarded. Discontinuties occur at the
ends of the integration intervals, but the kinetics are unaf-
fected. For example, when T<12, the value of 12*INT(T/12) is 0.
The values of FLAG1 and FLAG2 are initially 0 and FLAG is 0. As
T approaches the end point T=12, the values of FLAG1 and FLAG2
both become 1 and FLAG is 0. At the end point when T=12, the
value of 12*INT(T/12) is 12 and both FLAG1 and FLAG2 are 0. The
discontinuity in mt1 and mt2 does not affect the FLAG variable
because FLAG1 and FLAG2 are both 0 or both 1 in the neighborhood
of the discontinuity, and the kinetics are continuous.
In these examples, EHC is driven by the dose events. The EHC changes
can also driven by the clock. Suppose every 12 hours, a new EHC cycle
begins. (The value 12 is chosen so that the two versions will give
the same predictions.) E.g., instead in mtimess.ctl, instead of
MTIME(1) = MTIME(1)+II
compute
MTIME(1) = MTIME(1)+12
In hillss.ctl, instead of
mt1=inter*INT(T/inter)+theta(8)
compute
mt1=12*INT(T/12)+theta(8)
There is no difference in the results.REFERENCES: Guide VI, section III.F.9; Guide VI, section V.F; Guide IV, section V.C.5
29 EXOGENOUS SUPPLEMENTATION EXAMPLE
DISCUSSION:
In this example, an oral “drug” is given exogenously, and it also exists as an endogenous substance. There is an unknown dosing history prior to the observation period (i.e., prior to time zero). This example illustrates how three sources of drug can be modeled: preexisting endogenous drug, pre-existing drug from an unknown prior dosing history, and drug from known doses. Any combination of the three could be modeled without the others.
The rate of endogenous drug production is assumed to be constant, with no feedback control of production. Thus endogenous drug is at steadystate, and, with linear kinetics, its effect is simply to add a constant increment to exogenous drug in the sampled compartment (the increment is modeled as theta(7)).
For the drug with unknown dosing history, it is assumed that the subject is at steady state with respect to this drug. This part of total drug is modeled by a steady state infusion dose into the depot compartment, ending at time 0, and having an unknown rate. The result of the SS dose is to introduce drug into all compartments of the system (not just the central compartment) because it is distributed throughout the system and is subject to elimination from the system. The unknown rate is modeled as theta(5). NONMEM will adjust theta(5) to best fit not only the “baseline” observation at time 0, but also the later observations.
Note that if samples are not taken sufficiently long after the time of the last dose ( > 4 half-lives), then theta(7) and theta(5) may not be separately identifiable. Note that the value of theta(7) may be determined by the residual concentration after all exogenous drug has disappeared.
A combined additive and ccv error model is used. Theta(8) is the ratio of the C.V. of the proportional component to the standard deviation of the additive component.
Any ADVAN/TRANS combination could be used. Population data could also be modeled in this manner, with eta variables in the $PK block.
$PROBLEM Example of pre-existing drug.
$INPUT ID TIME DV AMT SS II RATE
$DATA DATA1
$SUBROUTINES ADVAN4 TRANS5$PK
AOB=THETA(1)
ALPHA=THETA(2)
BETA=THETA(3)
KA=THETA(4)
R1=THETA(5)
S2=THETA(6)$ERROR
FP=THETA(7)+F ; adds endogenous component
W=(1+THETA(8)*THETA(8)*FP*FP)**.5
Y=FP+W*ERR(1)Note that, if there are other doses into the depot compartment with modeled rates, it is necessary to assign a value to R1 conditionally. E.g.,
IF (TIME.EQ.0) THEN
R1=THETA(5) ; rate for SS infusion record at time 0ELSE
R1=....; rate of other kind of dose
ENDIFNote also that the combined additive and ccv error model can also be modeled using two random variables:
Y = F*(1+ERR(1)) + ERR(2)A fragment of the data follows. Record 1 specifies the SS infusion for the pre-existing drug, which ends at time 0. Record 2 gives the baseline observation. Record 3 specifies an oral bolus dose. Record 4 gives an observation.
1 0 . 0 1 0 -1
1 0 62.2 . . . .
1 0.01 . 95 . . .
1 0.50 235.93 . . . .REFERENCES: Guide V, section 8; Guide VI, section III.F.5
30 EXPECTATION EXAMPLE
These are examples of the use of the Marginal (MRG_) data item. Example 1 Suppose that the probability that a particular subject experiences a pain relief score of 2 is computed. Suppose also one wants to compute the (posterior population) expectation of the probability with each of 4 different bolus doses, not all of which are among those used to obtain observations. A fragment of the control file follows.
$INPUT ID DOSE MDV MRG_ ...
$PRED
....
Y = likelihood of observation given ETA
$EST METH=COND LAPLACE LIKELIHOOD
$TABLE DOSE
A fragment of the data follows, with comments following ";".
#ID DOSE DV MDV MRG_
100 5 2 1 1 ; PRED item is set to expectation for DOSE=5
100 10 2 1 1 ; PRED item is set to expectation for DOSE=10
100 20 2 1 1 ; PRED item is set to expectation for DOSE=20
100 40 2 1 1 ; PRED item is set to expectation for DOSE=40
1 3 1 0 0
1 10 1 0 0
1 25 2 0 0
1 30 2 0 0
2 3 1 0 0
2 10 1 0 0
2 25 1 0 0
2 30 1 0 0
3 3 1 0 0
3 10 2 0 0
3 25 2 0 0
3 30 3 0 0
... etc ...
Example 2
This example produces a plot of four residuals, formed by the differ-
ences between the raw-data-averages and their (posterior population)
expectations, versus the doses used to obtain the data.
$INPUT ID DOSE MDV MRG_ RAW_
$PRED
...
IF (ICALL.EQ.5) THEN
Y = expectation of observation given ETA
ELSE
Y = likelihood of observation given ETA
ENDIF
$EST METHOD=COND LAPLACE LIKELIHOOD
$SCAT RES VS DOSE
#ID DOSE DV MDV MRG_ RAW_
100 5 . 1 1 1
100 10 . 1 1 1
100 20 . 1 1 1
100 40 . 1 1 1
1 5 1 0 0 0
1 10 1 0 0 0
1 20 2 0 0 0
1 40 2 0 0 0
... etc ...
Example 3
This example produces a plot of four residuals, formed by the differ-
ences between the proportion of subjects in the data set with pain
relief score 2 and the (posterior population) expectation of the prob-
ability that a subject experiences a pain score of 2, versus the doses
used to obtain the data.
$INPUT ID DOSE MDV MRG_ RAW_
$PRED
...
Y = likelihood of observation given ETA
IF (ICALL.EQ.6) THEN
DVR=DV
DV=0
IF (DVR.EQ.2) DV=1
ENDIF
$EST METHOD=COND LAPLACE LIKELIHOOD
$SCAT RES VS DOSE
#ID DOSE DV MDV MRG_ RAW_
100 5 2 1 1 1
100 10 2 1 1 1
100 20 2 1 1 1
100 40 2 1 1 1
1 5 1 0 0 0
1 10 1 0 0 0
1 20 2 0 0 0
1 40 2 0 0 0
2 5 1 0 0 0
2 10 1 0 0 0
2 20 1 0 0 0
2 40 3 0 0 0
3 5 1 0 0 0
3 10 2 0 0 0
3 20 2 0 0 0
3 40 3 0 0 0
... etc ...
(See mrg, expectation block, data average block, raw).REFERENCES: none
31 FINEDATA EXAMPLE
This example is new to NONMEM 7. It is an example of the use of the finedata utility program to fill in missing values of an independent variable (covariate). Other examples (infn1, infn2) show the use of the INFN routines for interpolation. In data file THEOPP used for CONTROL5, where only one value of WT is present, the INFN routine is used to fill in missing values of WT. Presumably these examples would be part of a NONMEM run involving Estimation. In this example, finedata is used to fill in the missing values of WT and create a new data set, THEOPPfull. This data set can be used in subsequent NONMEM runs for estimation.
USAGE:
finedata theofine.ctl
Contents of theofine.ctl are as follows:
$PROB THEOPHYLLINE POPULATION DATA
$INPUT ID DOSE=AMT TIME CP=DV WT
$DATA THEOPP$FINEDATA AXIS=TIME FILE=THEOPPfull WT=LIN MISSING=. NEVAL=-1
(See control5). (See infn1 example). (See infn2 example).
REFERENCES: Guide Introduction_7
32 FOCE_PARALLEL EXAMPLE
The example control stream file foce_parallel.ctl is in the ..directory.
; Used for comparing single versus parallel computing for FOCE method.
;$SIZES LVR=30
;$SIZES LTH=15
;$SIZES LIM1=100
$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
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]
;$EST METHOD=CHAIN FILE=foce_parallel.chn ISAMPLE=3 NSAMPLE=0
$EST METHOD=1 INTERACTION NSIG=3 PRINT=1
NOABORT FORMAT=s1PE23.16 MAXEVAL=9999 ; PARAFILE=parallel_file.pnm
$COV MATRIX=S PRINT=E UNCONDITIONAL ; PARAFILE=OFFREFERENCES: Guide Introduction_7
33 IGNORE_ACCEPT EXAMPLE
The options IGNORE=(list) and ACCEPT=(list) of the $DATA record allow for simple selection criteria for records that are to be excluded or included, respectively, in the NONMEM data set. A list has the format: cond,cond,cond,… where each cond is a test of a data item: “label operator value”. The commas are implied .OR. operators. The .AND. and .NOT. operators are not permitted. Suppose an .AND. is needed. E.g., IGNORE=(GEN.EQ.1.AND.AGE.GT.60) If the expressions are simple, this may be done by negating both sides: .NOT.IGNORE =.NOT.(GEN.EQ.1.AND.AGE.GT.60) ACCEPT=.NOT.GEN.EQ.1.OR..NOT.AGE.GT.60 Thus, ACCEPT=(GEN.NE.1,AGE.LE.60) But logical operators inside parenthesis are not permitted, so that the following cannot be coded:
$DATA ACCEPT=((A==1.OR.A==2).AND.B<100)
Even if negated, this is:
IGNORE=(.NOT.(A==1.OR.A==2)).OR..NOT.B<100)
which is:
IGNORE=((A.NE.1.AND.A.NE.2).OR.B>=100)
There is still an .AND. inside parentheses.
There is a workaround that takes advantage of the fact that A takes
numeric values and there are no values between 1 and 2.
(A.NE.1.AND.A.NE.2)
can be implemented with an "or":
(A.LT.1.OR.A.GT.2).
The following are equivalent:
IGNORE=(B.GE.100) IGNORE=(A.LT.1) IGNORE=(A.GT.2)
IGNORE=(B.GE.100,A.LT.1,A.GT.2)
This workaround cannot handle continuous values of A or more compli-
cated conditions.
With NONMEM through 7.3, two separate NONMEM problems can be used, in
two separate runs, for a complicated condition. In the first problem,
abbreviated code is used to append a new variable to the data with a
value indicating ignore/accept; in the second run the variable is used
in the $DATA statement to choose from the data set.
Suppose ACC is the variable that is to have values 0/1 for
ignore/accept. The desired code is:
ACC=0
IF ((A == 1.OR.A == 2).AND.B<100) ACC=1
The above can be coded using the FORTRAN language. But NM-TRAN
abbreviated code is a subset of FORTRAN. Logical operators .NOT.,
.AND., .OR. may not be used within parentheses. There are two work-
arounds.
One is to clear the parentheses:
IF (A==1.AND.B<100.OR.A==2.AND.B<100) ACC=1
(this is always possible, no matter how complicted the conditional
expression).
The second is to use several statements. There are may ways to do
this. In the following, the .AND. of multiple conditions is false if
any of the conditions is false.
ACC=1
IF (A.NE.1.AND.A.NE.2) ACC=0
IF (B.GE.100) ACC=0
In prob1.ctl, the first problem appends ACC to the (new) data set. In
prob2.ctl, the $DATA record uses ACC to select records.
All required data items for NONMEM (and PREDPP, if used) must be
present in the original data file testab.dat. The control stream for
the first run:
$PROB test of accept Run #1
; Example of implementation of
; $DATA ACCEPT=((A == 1.OR.A == 2).AND.B<100)
$INPUT NO A B ID DV MDV
$DATA testab.dat IGNORE C
$PRED
; Implements:
; ACC=1 IF ((A == 1.OR.A == 2).AND.B < 100)
ACC=1
IF (A.NE.1.AND.A.NE.2) ACC=0
IF (B.GE.100) ACC=0
; The model for this run is unimportant. Keep it simple.
Y=THETA(1)+ETA(1)+EPS(1)
$THETA 1
$OMEGA 1
$SIGMA .4
$TABLE NO A B ID DV MDV ACC
NOPRINT NOAPPEND NOHEADER FILE=run1.tab
The control stream for the second run:
$PROB test of accept Run #2
$INPUT NO A B ID DV MDV ACC
; Must do a numeric test for ACC==1, because it appears in
; table file run1.tab as 1.0000E+00
$DATA run1.tab ACCEPT=(ACC.EQN.1)
$PRED
; the model in this problem is the model
; used for simulation, analysis, etc.
Y=THETA(1)+ETA(1)+EPS(1)
$THETA 1
$OMEGA 1
$SIGMA .4
$TABLE NO A B ID DV MDV
NOPRINT NOAPPEND NOHEADER FILE=run2.tab
With NONMEM 7.4, it is possible to combine both problems into one NON-
MEM run using the $TABLE record option "EXCLUDE_BY list". If any
variable in list is non-0, the record is excluded from the table file
probe.tab. The $DATA ACCEPT/IGNORE option is not needed. Because
probe.tab does not exist when NM-TRAN processes the second problem in
the control stream, $DATA .. NOOPEN is needed and a format specifica-
tion is needed.
$PROB test of NONMEM 7.4 exclude_by feature
$INPUT NO A B ID DV MDV
$DATA testab.dat IGNORE C
$PRED
; sets ACC=0 if (a==1.or.a==2).and.b<100)
; uses the clear the parenthesis approach.
ACC=1
IF (A==1.AND.B<100.OR.A==2.AND.B<100) ACC=0
; The following is the model for both problems.
; It may be more complicated than this.
Y=THETA(1)+ETA(1)+EPS(1)
$THETA 1
$OMEGA 1
$SIGMA .4
; probe.tab has only those records that have ACC=0.
$TABLE NO A B ID DV MDV EXCLUDE_BY ACC
NOPRINT NOAPPEND NOHEADER FILE=probe.tab
$PROB test of exclude_by #2
$INPUT NO A B ID DV MDV
$DATA probe.tab (10F12.0) NOOPEN
$THETA 1
$OMEGA 1
$SIGMA .4
; This problem should have $ESTIMATION and other tasks
$TABLE NO A B ID DV NOAPPEND NOHEADER NOPRINT FILE=probe2.tab
Here is a fragment of the data testab.dat
C NO A B ID DV MDV
1 1 100 1 0 0
2 2 100 1 0 0
3 3 100 1 0 0
4 1 101 1 0 0
5 2 101 1 0 0
6 3 101 1 0 0
7 1 99 1 0 0
8 2 99 1 0 0
9 3 99 1 0 0
(See $data, $table).REFERENCES: none
34 INFN_INTERPOLATION EXAMPLE 1
This example is adapted from Guide VI (PREDPP), Figure 37. It was revised for NONMEM 7.3. Earlier versions were untested and did not work correctly. This version allows both single-subject and population data. It requires a “missing independent variable” data item. In this example, the INFN subroutine is used to fill in (interpolate) missing values of WT in the THEOPPMV data set. THEOPPMV is a version of the data set THEOPP used in CONTROL5, (See control5). If WT is recorded on every record, then the $PK block can be simpler. The $PK block can be called with every record (no need for CALLFL=1) and no special code is needed to save the value of WT for subsequent records. Example infn2 shows how the INFN routine may be implemented in abbreviated code. Example fine1 shows how the finedata utility program can be used to fill in the values without the use of the INFN routine.
USAGE:
$PROB THEOPHYLLINE POPULATION DATA
$INPUT ID DOSE=AMT TIME CP=DV WT MISS
$DATA THEOPPMV
$SUBROUTINES ADVAN2 INFN=infnsub
The file infnsub contains the following:
! THIS EXAMPLE IS ADAPTED FROM GUIDE VI (PREDPP), FIGURE 37.
! FOR HELP FILE infn1.exa
! MODIFIED BY AJB 9/2012
! INFN ROUTINE FOR COMPUTING LINEARLY INTERPOLATED VALUES
! OF AN INDEPENDENT VARIABLE V. ILLUSTRATES USE OF ROUTINE PASS.
! USE FOR SINGLE-SUBJECT OR POPULATION DATA
! WORKS WITH 1 OR MORE NON-MISSING VALUES.
! A SINGLE SUCH VALUE MAY BE ANYWHERE IN DATA SET.
! DATREC(UI)=TIME DATA ITEM
! DATREC(VI)=INDEPENDENT VARIABLE DATA ITEM
! DATREC(MI)=MISSING INDEPENDENT VARIABLE DATA ITEM
! =0 INDEP VAR NOT MISSING
! >0 IF THIS DATA RECORD HAS MISSING INDEP VAR AND:
! =1 THIS DATA RECORD PRECEDES FIRST NON-MISSING VALUE
! =3 THIS DATA RECORD FOLLOWS LAST NON-MISSING VALUE
! =2 THIS DATA RECORD IS BETWEEN NON-MISSING VALUES
!
SUBROUTINE INFN (ICALL,THETA,DATREC,INDXS,NEWIN)
USE SIZES, ONLY: ISIZE,DPSIZE,NO,MAXIDS
! NOTE THAT NEWIND IS NOT OBTAINED FROM THE SUBROUTINE ARGUMENT
USE NMPRD_INT, ONLY: NEWIND=>NWIND
INTEGER(KIND=ISIZE), INTENT(IN) :: ICALL,INDXS(*),NEWIN
REAL(KIND=DPSIZE), INTENT(IN) :: THETA(*)
REAL(KIND=DPSIZE), INTENT(IN OUT) :: DATREC(*)
REAL(KIND=DPSIZE) :: U,V
DIMENSION U(NO,MAXIDS),V(NO,MAXIDS),DEPVAR(MAXIDS),LASTI(MAXIDS)
INTEGER UI,VI,MI,IS,LASTI
! IMPORTANT: CHANGE UI, VI, MI ACCORDING TO LAYOUT OF DATA SET
UI=3
VI=5
MI=6
!
IF (ICALL.NE.1) RETURN
I=0
IS=0
! INITIALIZE PASS
MODE=0
CALL PASS (MODE)
MODE=2
! PASS THROUGH DATA
5 CALL PASS (MODE)
IF (MODE.EQ.0) GO TO 10
IF (NEWIND<2) THEN
I=0
IS=IS+1
ENDIF
! IF INDEP VAR IS PRESENT, STORE TIME AND VALUE
IF (DATREC(MI).EQ.0.) THEN
I=I+1
U(IS,I)=DATREC(UI)
V(IS,I)=DATREC(VI)
LASTI(IS)=I
ENDIF
IF (I == 1) THEN ! SAVE VALUE IN CASE ONLY RECORDED ONCE
DEPVAR(IS)=V(IS,1)
ENDIF
GO TO 5
! INITIALIZE PASS A SECOND TIME
10 I=0
IS=0
MODE=0
CALL PASS (MODE)
MODE=2
! PASS THROUGH DATA A SECOND TIME
15 CALL PASS (MODE)
IF (MODE.EQ.0) RETURN
IF (NEWIND<2) THEN
I=0
IS=IS+1
ENDIF
! IF INDEP VAR IS MISSING AND ONLY RECORDED ONCE, COPY IT
IF (DATREC(MI).NE.0. .AND. LASTI(IS) == 1 ) THEN
DATREC(VI)=DEPVAR(IS)
GO TO 15
ENDIF
! IF INDEP VAR IS MISSING, STORE INTERPOLATED VALUE
IF (DATREC(MI).EQ.0.) THEN
I=I+1
ELSE
IF (DATREC(MI).EQ.1.) THEN ! EXTRAPOL. FROM FIRST 2 VALUES
K=1
L=2
ELSEIF (DATREC(MI).EQ.2.) THEN ! INTERPOL. FROM BEFORE AND AFTER
K=I
L=I+1
ELSEIF (DATREC(MI).EQ.3.) THEN ! EXTRAPOL. FROM LAST 2 VALUES
K=I-1
L=I
ENDIF
B=(V(IS,K)-V(IS,L))/(U(IS,K)-U(IS,L))
DATREC(VI)=V(IS,K)+B*(DATREC(UI)-U(IS,K))
ENDIF
GO TO 15
END
Data file THEOPPMV is identical to THEOPP, but has an extra data item
MISS (Missing Dependent Variable). The data for the first subject is:
1 4.02 0. . 79.6 0
1 . 0. .74 . 3
1 . 0.25 2.84 . 3
1 . 0.57 6.57 . 3
1 . 1.12 10.5 . 3
1 . 2.02 9.66 . 3
1 . 3.82 8.58 . 3
1 . 5.1 8.36 . 3
1 . 7.03 7.47 . 3
1 . 9.05 6.89 . 3
1 . 12.12 5.94 . 3
1 . 24.37 3.28 . 3
(See infn2 example).
(See fine1 example).REFERENCES: Guide VI, section VI.A , Figure 37
35 INFN_INTERPOLATION EXAMPLE 2
This example is new to NONMEM 7. It is an improved version of infn1.exa. In this example, instead of using a stand-alone INFN subroutine, $INFN abbreviated code is used to fill in missing values of WT in the THEOPP data set. (See control5). If WT is recorded on every record, then the $PK block can be simpler. The $PK block can be called with every record (no need for CALLFL=1) and no special code is needed to save the value of WT for subsequent records. No “missing independent variable” data item is needed, so the THEOPP data file can be used as-is. The $ABBR abbreviated code is located in file infnabbr for convenient re-use. This is also an example of the
$ABBR REPLACE feature, which allows easy replacement of character
strings in the included code, and is also an example of declared sub-
scripted variables.
Example fine1 shows how the finedata utility program can be used to
fill in the values without the use of the INFN routine.USAGE:
$PROB THEOPHYLLINE POPULATION DATA
$INPUT ID DOSE=AMT TIME CP=DV WT
$DATA THEOPP
$SUBROUTINES ADVAN2
$ABBR REPLACE INTVBL=WT
$ABBR REPLACE NULLVAL=0.0
INCLUDE 'infnabbr'
The file infnabbr contains the following:
; INFNABBR: ABBREVIATED CODE FOR HELP infn2.exa
; THIS INFN BLOCK INTERPOLATES AND EXTRAPOLATES VALUES OF A SINGLE
; DATA ITEM CALLED INTVBL.
; BY DEFAULT, THE AXIS VARIABLE IS CALLED TIME.
; THE CHARACTER STRING NULLVAL DENOTES A NULL VALUE.
;
; IF THERE IS ONLY ONE NON-NULL VALUE, THIS VALUE IS COPIED TO
; ALL NULL VALUES. OTHERWISE:
; LINEAR INTERPOLATION IS PERFORMED FOR NULL VALUES THAT
; LIE BETWEEN TWO NON-NULL VALUES.
; LINEAR EXTRAPOLATION IS PERFORMED FOR VALUES BEFORE THE
; FIRST NON-NULL VALUE
; LINEAR EXTRAPOLATION IS PERFORMED FOR VALUES AFTER THE
; LAST NON-NULL VALUE
; INFNABBR CAN BE USED WITH BOTH SINGLE
; SUBJECT AND POPULATION DATA SET.
; WRITTEN BY AJB 9/2012
;
; TYPICAL USAGE:
; $ABBR REPLACE INTVBL=WT
; $ABBR REPLACE NULLVAL=0.0
; $ABBR REPLACE TIME=T (ONLY IF T IS THE AXIS VARIABLE RATHER THAN TIME)
; INCLUDE 'INFNINTERP'
;
;
$ABBR DECLARE U(MAXIDS,NO), V(MAXIDS,NO)
$ABBR DECLARE INTEGER IS
$ABBR DECLARE INTEGER I, INTEGER J, INTEGER L
$ABBR DECLARE IVAL(NO)
$INFN
; $INFN FOR COMPUTING INTERPOLATED VALUES OF AN
; INDEPENDENT VARIABLE NAMED INTVBL
; E.g., INTVBL MIGHT BE WEIGHT
; ASSUME THAT INTVBL = NULLVAL WHEN MISSING.
; IS=SUBJECT #
; IVAL(IS)= # OF NON-NULL VALUES OF THE INDEP VAR. FOR SUBJECT IS
; VI(IS,I) = ITH. NON-NULL VALUE OF THE INDEPENDENT VARIABLE
; UI(IS,I) = TIME OF THE ITH. NON-NULL VALUE OF THE INDEPENDENT VARIABLE
IS=0
IF (ICALL.EQ.1) THEN
; FIRST PASS. SAVE VALUES OF TIME & INDEP VAR. WHEN INDEP VAR. IS NON-NULL
DO WHILE (DATA)
IF (NEWIND < 2) THEN ; INITIALIZE NEW INDIVIDUAL
IS=IS+1
IVAL(IS)=0
ENDIF
IF (INTVBL /= NULLVAL) THEN ; SAVE NON-NULL VALUE
IVAL(IS)=IVAL(IS)+1
I=IVAL(IS)
U(IS,I)=TIME
V(IS,I)=INTVBL
ENDIF
ENDDO
; SECOND PASS
IS=0
DO WHILE (DATA)
IF (NEWIND < 2) THEN ; INITIALIZE NEW INDIVIDUAL
IS=IS+1
I=0
ENDIF
; IF INDEP VAR IS MISSING AND WAS ONLY NON-NULL ONCE, COPY IT
IF (INTVBL == NULLVAL .AND. IVAL(IS) == 1) THEN
INTVBL=V(IS,1)
ELSE
IF (INTVBL /= NULLVAL) THEN ; COUNT ANOTHER NON-NULL VALUE
I=I+1
ELSE ; CURRENT RECORD HAS NULL VALUE
IF (I==0) THEN ; EXTRAPOLATE FROM FIRST 2 VALUES
J=1
L=2
ELSEIF (I>0 .AND. I<IVAL(IS)) THEN ; INTERPOL. BEFORE AND AFTER
J=I
L=I+1
ELSEIF (I==IVAL(IS)) THEN ; EXTRAPOL. FROM LAST 2 VALUES
J=I-1
L=I
ENDIF
SLOPE=(V(IS,J)-V(IS,L))/(U(IS,J)-U(IS,L))
INTVBL=V(IS,J)+SLOPE*(TIME-U(IS,J))
ENDIF
ENDIF
ENDDO
ENDIF
(See control5).
(See infn1 example).
(See fine1 example).REFERENCES: Guide VI, section VI.A , Figure 37; Guide Introduction_7
36 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 (52,*) THETAFR
WRITE (53,*) SETHET
WRITE (53,*) SETHETR
WRITE (54,*) OMEGA(BLOCK)
WRITE (55,*) SEOMEG(BLOCK)
WRITE (56,*) SIGMA(BLOCK)
WRITE (57,*) SESIGM(BLOCK)
WRITE (58,*) IERE,IERC
ENDIFREFERENCES: Guide VI, section VI.A
37 INTEROCCASION_VARIABILITY EXAMPLE
This example is based on work described in Karlsson and Sheiner, JPB 21(6):735-750 (1994). The model accounts for random variability in a subject’s parameters between study occasions. It illustrates a data file in which EVID data items are present. In the fragment of data below, EVID has values 1 (dose event), 0 (observation event), and 3 (reset event). The reset event is present at the start of the subject’s second study occasion. It causes PREDPP to re-initialize the kinetic system as if for a new individual, and permits time to be reset to 0. Note that the subject has the same dosage regimen prior to both study occasions (a bolus dose of AMT 1000 every 12 hours leading to steady state). However, it would have been possible to specify a different dosing regimen prior to the second occasion. In the abbreviated code, EXIT statements are used to constrain TVCL and TVV (typical value of clearance and volume) to positive values. Because the models for TVCL and TVV include both thetas and covariates, it is difficult to bound the thetas so as to insure that TVCL and TVV are non-negative for all subjects. The NOABORT option is present on the $ESTIMATION record. Hence NONMEM will perform PRED error recovery, and will avoid values of theta (subsequent to the initial values) that result in non-positive TVCL or TVV for any subject. Explanation of user data items: AGE age RACE race (1 or 2) HT height HCTZ hydrochlorothiazide (0 or 1) OCC occasion (0 or 1)
;THIS IS A NMTRAN CONTROL STREAM For THE PRAZOSIN PK ANALYSIS
$PROB prazosin data
$INPUT ID AGE RACE HT HCTZ AMT TIME DV EVID SS II OCC
$DATA praz21
$SUBROUTINE ADVAN2 TRANS2
$PK
OCC2=1-OCC
R=0
IF (RACE.EQ.2) R=1
TMP=(HT-160)*THETA(1)
IF (TMP.LE.0.) TMP=0.0
TMP2=(AGE-60)*THETA(2)
IF (AGE.LE.60.) TMP2=0.0
TVCL=THETA(3)+TMP-TMP2
TVCL=TVCL+THETA(4)*R+THETA(5)*HCTZ
IF (TVCL.LE.0.) EXIT 1 100
TVV=THETA(6)+HCTZ*THETA(7)
IF (TVV.LE.0.) EXIT 1 200
TVKA=THETA(8)
CL=TVCL*EXP(ETA(3)*OCC+ETA(5)*OCC2+ETA(1))
V =TVV *EXP(ETA(4)*OCC+ETA(6)*OCC2+ETA(2))
KA=TVKA*EXP(ETA(8)*OCC+ETA(9)*OCC2+ETA(7))
S2=V
$THETA (0,0.4) (0.,0.3) 20. 4 -4.
$THETA 80 -25
$THETA (0,1.)
$OMEGA BLOCK(2) 0.2 0.1 0.2
$OMEGA BLOCK(2) 0.2 0.1 0.2
$OMEGA BLOCK(2) SAME
$OMEGA 1.2
$OMEGA BLOCK(1) .5
$OMEGA BLOCK(1) SAME
$ERROR
Y=F*(1+EPS(1))
$SIGMA 0.1
$EST NOABORT SIG=3 MAX=3000 PRINT=10 POSTHOC
$COVARIANCE
$TABLE ID OCC ETA1 ETA3 ETA5 ETA2 ETA4 ETA6 NOPRINT FILE=praz32.tab
Data for one subject follow.
21 57 1 180 1 1000 0 0.62 1 1 12 0
21 57 1 180 1 0.0 1 12.06 0 0 0 0
21 57 1 180 1 0.0 3 6.81 0 0 0 0
21 57 1 180 1 0.0 4 4.89 0 0 0 0
21 57 1 180 1 0.0 5 4.04 0 0 0 0
21 57 1 180 1 0.0 6 2.82 0 0 0 0
21 57 1 180 1 0.0 7 2.72 0 0 0 0
21 57 1 180 1 0.0 8 1.78 0 0 0 0
21 57 1 180 1 0.0 10 1.07 0 0 0 0
21 57 1 180 1 0.0 12 0.75 0 0 0 0
21 57 1 180 1 0.0 0 2.83 3 0 0 0
21 57 1 180 1 1000 0 2.83 1 1 12 1
21 57 1 180 1 0.0 1 6.27 0 0 0 1
21 57 1 180 1 0.0 2 12.72 0 0 0 1
21 57 1 180 1 0.0 3 9.99 0 0 0 1
21 57 1 180 1 0.0 4 6.90 0 0 0 1
21 57 1 180 1 0.0 5 6.59 0 0 0 1
21 57 1 180 1 0.0 6 4.07 0 0 0 1
21 57 1 180 1 0.0 7 3.11 0 0 0 1
21 57 1 180 1 0.0 8 2.79 0 0 0 1
21 57 1 180 1 0.0 10 1.72 0 0 0 1
21 57 1 180 1 0.0 12 1.17 0 0 0 1
(See $estimation, abbreviated code, evid).REFERENCES: Guide VI, section V.B , VII.C.2; Guide IV, section III.B.10 , IV.G; Guide V, section 12.4.15
38 INITIAL STEADY STATE EXAMPLE
This is an example of the use of the Initial Steady State (I_SS) reserved variable in $PK. It was provided by Nick Holford, Dept Pharmacology & Clinical Pharmacology, University of Auckland. Glucose and insulin are both endogenous substances, whose rates of formation are given by terms RGLU and RINS in the differential equations. In the absence of dose events, both glucose and insulin are assumed to be at steady state levels. Exogenous doses of both glucose and insulin are also present in the first subject’s data. It is possible to intialize the steady state condition with a dose record having AMT=0, RATE=0, SS=1, but this is artificial. The initial condition that the system is at steady state is a feature of the model, not of the data. In this example, the $PK block contains the statement I_SS=1. This causes the system to be initialized appropriately. The statement I_SS=1 could also have been placed among the
$MODEL statements. However the initialization is performed, the
result is the same: non-zero concentrations GLU and INS are computed
at TIME=0.
The SS data item need not be present in the data set for the I_SS com-
putation. It is present in this example in order to demonstrate that
SS may also be computed using a dose record (see the third and fourth
subjects' data).
Note that SIGMA and OMEGA are fixed to 0 for purposes of illustration,
so that the simulated values are exactly equal to those that arise
from the model without the random variability due to ETA and EPS.
NONMEM issues warning messages ("INITIAL ESTIMATE OF WITHIN INDIVIDUAL
VARIANCE IS ZER0"). These may be ignored.
This example also illustrates the use of A(n) on the right in the
$ERROR block to obtain the compartment amounts that are used for the
computation of concentrations GLU and INS.
$PROB nicksiss6.ctl I_SS=1 in $PK ADVAN6
$DATA Nicksiss.dat
$INPUT ID TIME CMT AMT RATE SS DV
$SIM (200070927) ONLYSIM NSUB=1
;Silber HE, Jauslin PM, Frey N, Gieschke R, Simonsson US, Karlsson MO.
;An integrated model for glucose and insulin regulation in healthy
;volunteers and type 2 diabetic patients following intravenous glucose
;provocations. J Clin Pharmacol. 2007 Sep;47(9):1159-71.
$THETA
40 ; POP_RGLU MMOL/H/70KG - glucose input rate
40 ; POP_VGLU L/70KG - volume of distribution for glucose
5 ; POP_CLGLU L/H/70KG - glucose clearance
1 ; POP_EMXGLU - EMAX for effect of glucose on insulin input
10 ; POP_C50GLU MMOL/L - EC50 for effect of glucose on Insulin input
5 ; POP_HILGLU - Hill parameter for effect of glucose on Insulin input
5000 ; POP_RINS PMOL/H/70KG - insulin input rate
5 ; POP_VINS L/70KG - volume of distribution for insulin
70 ; POP_CLINS L/H/70KG - insulin clearance
2 ; POP_EMXINS - EMAX for effect of insulin on glucose clearance
50 ; POP_C50INS PMOL/L - EC50 for effect of insulin on glucose clearance
2 ; POP_HILINS - Hill parameter for effect of insulin on glucose clearance
$OMEGA
0 FIX ;PPV_RGLU
0 FIX ;PPV_RINS
$SIGMA
0 FIXED ;0.1 ;G_EXP_RUV
0 FIXED ;1 ;G_ADD_RUV MMOL/L
0 FIXED ;0.1 ;I_EXP_RUV
0 FIXED ;1 ;I_ADD_RUV PMOL/L
$SUBR ADVAN6 TOL=3
$MODEL
COMP (GLUCOSE)
COMP (INSULIN)
$PK
I_SS=1
; GLUCOSE
RGLU=THETA(1)*EXP(ETA(1))
VGLU=THETA(2)
CLGLU=THETA(3)
EMXGLU=THETA(4)
C50GLU=THETA(5)
HILGLU=THETA(6)
;INSULIN
RINS=THETA(7)*EXP(ETA(2))
VINS=THETA(8)
CLINS=THETA(9)
EMXINS=THETA(10)
C50INS=THETA(11)
HILINS=THETA(12)
S1=VGLU
S2=VINS
$DES
DGLU=A(1)/VGLU
DINS=A(2)/VINS
DGLUH=DGLU**HILGLU
DGEFF=EMXGLU*DGLUH/(C50GLU**HILGLU+DGLUH) ; effect of glucose on insulin
DINSH=DINS**HILINS
DIEFF=EMXINS*DINSH/(C50INS**HILINS+DINSH) ; effect of insulin on glucose
DADT(1)=RGLU - CLGLU*(1+DIEFF)*DGLU ; glucose
DADT(2)=RINS*(1+DGEFF)- CLINS*DINS ; insulin
$ERROR
GLU=A(1)/VGLU ; glucose concentration
INS=A(2)/VINS ; insulin concentration
IF (CMT.EQ.1)THEN
Y=GLU*(1+ERR(1))+ERR(2)
ENDIF
IF (CMT.EQ.2)THEN
Y=INS*(1+ERR(3))+ERR(4)
ENDIF
$TABLE ID TIME CMT AMT RATE SS GLU INS
ONEHEADER NOPRINT FILE=nicksiss6.tab
The data (Nicksiss.dat) follows.
#ID,TIME,CMT,AMT,RATE,SS,DV
1 0 1 0 0 0 .
1 0 2 0 0 0 .
1 .1 1 25 0 0 .
1 1.0 2 10 0 0 .
1 3 1 0 0 0 .
1 3 2 0 0 0 .
1 100 1 . . . .
1 100 2 . . . .
2 0 1 0 0 0 .
2 0 2 0 0 0 .
2 100 1 . . . .
2 100 2 . . . .
3 0 1 0 0 1 .
3 100 1 . . . .
3 100 2 . . . .
4 0 1 0 0 1 .
4 100 1 . . . .
4 100 2 . . . .
(See Initial Steady State,I_SS,ISSMOD)
(See i_ss, initial_condition).REFERENCES: none
39 LOGISTIC REGRESSION EXAMPLE
Here is a typical NM-TRAN control stream for dichotomous response population data. Such data is referred to as “odd-type” data. Note that the eta is a population eta but there are no epsilons.
$PROB DICHOTOMOUS RESPONSE
$DATA data
$INPUT ID DOSE DV
$PRED
LOGIT=THETA(1)+THETA(2)*DOSE**THETA(3)+ETA(1)
A=EXP(LOGIT)
P=A/(1+A)
IF (DV.EQ.1) THEN
Y=P
ELSE
Y=1-P
ENDIF
$THETA .1 (25,100) (1,3,5)
$OMEGA 1
$ESTIMATION MET=COND LAPLACE LIKE MAX=500
$SCAT ETA(1) VS DOSE
Here is a typical NM-TRAN control stream when, one wants to first sim-
ulate the data.
$PROB DICHOTOMOUS RESPONSE
$DATA data
$INPUT ID DOSE DV
$PRED
LOGIT=THETA(1)+THETA(2)*DOSE**THETA(3)+ETA(1)
A=EXP(LOGIT)
P=A/(1+A)
IF (ICALL.EQ.4) THEN
CALL RANDOM (2,R)
IF (R.LE.P) DV=1
IF (R.GT.P) DV=0
ENDIF
IF (DV.EQ.1) Y=P
IF (DV.EQ.0) Y=1-P
$THETA .1 (25,100) (1,3,5)
$OMEGA 1
$SIMULATION (7755399) (45211 UNIFORM)
$ESTIMATION MET=COND LAPLACE LIKE MAX=500
$SCAT ETA(1) VS DOSEREFERENCES: None
40 MICHAELIS_MENTEN EXAMPLE
The kinetic system in this example is composed of depot, peripheral, and central compartments, and includes Michaelis-Menten Elimination from the central compartment. Either ADVAN6 or ADVAN8 could be used. As with any ADVAN, one or more doses could be input to any of the compartments. Steady state doses are also possible.
$PROBLEM EXAMPLE OF M-M ELIMINATION WITH DEPOT AND PERIPH
$INPUT ID TIME AMT DV
$DATA datafile
$SUBROUTINES ADVAN6 TRANS1 TOL=4
$MODEL COMP=(DEPOT,DEFDOS),COMP=(CENTRAL,DEFOBS) COMP=(PERIPH)
$PK
VM = THETA(1)*EXP(ETA(1))
KM = THETA(2)*EXP(ETA(2))
S2 = THETA(3)*EXP(ETA(3))
K12 = THETA(4)*EXP(ETA(4))
K23 = THETA(5)*EXP(ETA(5))
K32 = THETA(6)*EXP(ETA(6))
$ERROR
Y = F + ERR(1)
$DES
C2 = A(2)/S2
DADT(1) = -K12*A(1)
DADT(2) = K12*A(1) -K23*A(2) +K32*A(3) -C2*VM/(KM+C2)
DADT(3) = K23*A(2) -K32*A(3)
; requires $THETA, $OMEGA, $SIGMA etc.
If the units of KM is amount, rather than concentration, the $DES
block need not compute C2. The change is as follows:
DADT(2) = K12*A(1) -K23*A(2) +K32*A(3) -A(2)*VM/(KM+A(2))
(ADVAN10 assumes the units of KM is amount.)REFERENCES: Guide VI, section VII.C.6 , VII.C.10
41 MIXTURE_MODEL EXAMPLE
In this example of a mixture model, it is assumed that some unknown fraction of the population has one set of typical values of CL and Vd, and that the remaining fraction has another set of typical values. Reserved variables of interest: MIXNUM NONMEM sets MIXNUM to 1 or 2 according to whether PK is to compute parameter values for the 1st or 2nd subpopulations, respectively. MIXEST For each individual, NONMEM computes an estimate of the number of the subpopulation of which each individual is a member, and stores this estimate in the integer variable MIXEST. EST In the $PK block, MIXEST is stored in the PK-defined variable EST, which enables MIXEST to be displayed in tables and scatterplots.
$PROB PHENOBARB POPULATION DATA MIXTURE MODEL
$DATA DATA2 (6F7.0)
$INPUT ID TIME AMT WT APGR CP=DV
$SUBROUTINES ADVAN1 TRANS2
$PK
CALLFL=1
EST=MIXEST
IF (MIXNUM.EQ.2) THEN
CL=THETA(2)*THETA(1)*EXP(ETA(3))
V=THETA(4)*THETA(3)*EXP(ETA(4))
ELSE
CL=THETA(1)*EXP(ETA(1))
V=THETA(3)*EXP(ETA(2))
ENDIF
S1=V
$ERROR
Y=F*EXP(EPS(1))
$MIX
P(1)=THETA(5)
P(2)=1.-THETA(5)
NSPOP=2
$THETA (0,.0047) (0,1) (0,.99) (0,1) (0,.5,1)
$OMEGA BLOCK(2) .05 .01 .03
$OMEGA BLOCK(2) SAME
$SIGMA .02
$ESTM MAXEVAL=500 PRINT=5
$TABLE ID EST FIRSTONLY NOAPPEND
(See displayed_values, $mix, mixnum mixest).REFERENCES: Guide VI, section III.L.2 , Figures 6, 23
42 MODEL TIME EXAMPLES
These are examples of the use of model times in PREDPP. Enterohepatic Recycling This fragment of abbreviated code may be used to model EHC. The transfer of drug from compt. 4 to 1 is controlled by FLAG, which is 1 between the times specified by MTIME(1) and MTIME(2), and is 0 otherwise.
$PK
MTIME(1)=THETA(8)
MTIME(2)=MTIME(1)+THETA(9)
....
$DES
FLAG=MPAST(1)-MPAST(2)
DADT(1)=-KA*A(1)+K41*A(4)*FLAG
DADT(4)=K1G*A(2)-K41*A(4)*FLAG
....
Flexible Input Rate
This fragment of abbreviated code may be used to model a flexible
modeled infusion rate R1.
The rate is 400*EXP(ETA(1)) from time 0 to 1.5
The rate is 300*EXP(ETA(2)) from time 1.5 to 2.5
The rate is 200*EXP(ETA(3)) from time 2.5 till end of infusion.
In the three assignment statments for R1, exactly one of the
right-side expressions is non-zero at each call to PK.
$PK
MTIME(1)=1.5
MTIME(2)=2.5
R1= 400*EXP(ETA(1))*(1-MPAST(1))
R1=R1+300*EXP(ETA(2))*(MPAST(1)-MPAST(2))
R1=R1+200*EXP(ETA(3))*MPAST(2)
Changing a Model Time
This fragment of abbreviated code shows how a model time parame-
ter can be changed. Suppose there are events at times 0 and 10
but one wants to advance in increments of 1 with stops at times
1, 2, 3, .... , 9.
IF (TIME.EQ.0) TEMP=0
TEMP=TEMP+1
MTIME(1)=TEMP
MTDIFF=1
Changing Absorption Rate
This fragment of abbreviated code shows how the Absorption Rate
KA can be changed.
KA is XKA1 from time 0 to THETA(5).
KA is XKA2 after time THETA(5).
$PK
MTIME(1)=THETA(5) ; change point for KA
XKA1=THETA(3)*EXP(ETA(3))
XKA2=THETA(4)*EXP(ETA(4))
KA=XKA1*(1-MPAST(1))+XKA2*(MPAST(1))
(See mtime).
See also "Cirdadian example: Examples Using MTIME to Model Periodic
Discontinuities in $DES".REFERENCES: None
43 MOVING AREA UNDER THE CURVE EXAMPLE
The following set of examples was developed in response to a question to the NONMEM Users about modeling a moving value of the area under the curve. The original question (from Pavel Belo) asked: “Efficacy is frequently considered a function of AUC. A disadvantage of this model of efficacy is that the effect is irreversable because AUC of concentration can only increase; it cannot decrease. In many cases, a more meaningful model is one where AUC is calculated from time t -a to t (kind of”moving average”), where t is time in the system of differential equations (variable T in NONMEM).” DUPLICATE SYSTEM METHOD Bob Bauer proposed a data set DELAYDATA and control file for computing a moving value of AUC. Bob wrote “In the following simple absorption model example developed by me and Alison Boeckmann for illustration purposes, compartments 1, 2, and 3 are the”real time” depot, central, auc, and compartments 4,5,6 are the “delayed time” depot, central, auc. So, the base model (non-time delay) system (compartments 1,2,3) is replicated (compartments 4,5,6) for the time delay portion. In addition, the data set duplicates the dose information of compartment 1 into compartment 4, and setting ALAG4 to a non-zero value in the control stream file provides a lag time to any doses inputted into compartment 4 (so this would take care of multiple dose problems as well). This allows for assessment and availability of AUC(t) and AUCT(t-time-delay) at any time t.” DELAYDATA and aucdelay1.ctl The data for one individual (DELAYDATA) and the control file (aucdelay1.ctl) from the email are provided in the examples directory of NONMEM 7.4 and higher. Note that aucdelay1.ctl and all the examples discussed below may be used with data from more than one individual and different numbers of event records per individual. In aucdelay1.ctl, the model for time delay is TDY=THETA(1)*EXP(ETA(1)) and the variable DAUC (“difference in area under the curve”) gives the difference in the area under the curve DAUC=A(3)-A(6) ; AUC(T)-AUC(T-TDY) One significant modification to the control file aucdelay1.ctl was made. The $ERROR model was changed from Y=F+EPS(1) to Y=F+DAUC+EPS(1) This demonstrates that DAUC can modify Y (although the relationship would not usually be this simple). SINGLE SYSTEM METHOD In aucdelay1.ctl the value of DAUC is available continuously at every time value T during the integration, but in fact it is not used in the differential equations themselves. DAUC is used only in the $ERROR block. In this case, an alternate method is possible. Only the “real time” single system of differential equations is integrated, and the needed values of AUC are captured “on the fly” in the $PK block. This should result in faster run times because fewer differential equations need to be numerically integreated. The model event time variable (MTIME) is used instead of absorption lag (ALAG) and is modelled with the same time delay variable TDY. For each event record, MTIME is used to interrupt the integration at the time of interest TIME-TDY. Variables called SAVEA are used to capture the value of A(3) at model event times for use in the $ERROR block. The general algorithm is rather complicated, so a simple algorithm will be described first. SINGLE SYSTEM METHOD: SIMPLE MODEL A simple model is possible when TDY is always smaller than the time between any two event records. aucdelay1S.ctl This is identical to aucdelay1.ctl except that THETA(1) has initial estimate .5 and is constrained to be between 0 and 1, and TDY is constrained to be less than 1, which is the change in time between any two event records. aucdelay2S.ctl This implements the simple model. Only the original (single) set of compartments 1-3 are defined. aucdelay2S.ctl uses the same data file DELAYDATA, but ignores records with CMT=4 (doses into the delayed dose compartment). For each event record, it is necessary to have available NEXTT, the TIME of the next event record. This data item is not recorded in DELAYDATA. Instead, NEXTT is listed in the $INPUT record as an extra data item beyond the data items (columns) listed in the data file, and its entries and are created (“transgenerated”) during each run by code in the $INFN block. The control stream has the same constraints on THETA(1) and TDY as aucdelay1S.ctl. Only one model event time variable, MTIME(1) is needed. Model event time MTIME(1) is modelled as NEXTT-TDY and changes with every event record. SAVEA is a variable that saves the value of the AUC compartment A(3) when $PK is evaluated at the model event time. At such times, MTIME(1) is reset for the next record. The important lines of code are relatively simple, and use TSTATE (“State Time”; the time at which the state-vector A was last computed.)
$PK
TDY=THETA(1)*EXP(ETA(1))
MTIME(1)=NEXTT-TDY
IF (TSTATE==TIME-TDY) SAVEA=A(3)
...
$ERROR
DAUC=A(3)-SAVEA
...
The two control files generate the same predictions.
SINGLE SYSTEM METHOD: GENERAL MODEL
The model is more complicated if TDY may be greater than the time
between event records, because each value of compartment A(3) must be
saved until needed. For the ith. event record, with TIME(i), the
value of SAVEAi at TIME(i)-TDY is needed, and is obtained by using a
value of MTIME(i) equal to TIME(i)-TDY.
aucdelay2.ctl and aucdelay3.ctl
The general model is implemented in aucdelay2.ctl and aucdelay3.ctl,
These versions integrate a single system of differential equations.
File aucdelay3.ctl is a short-hand version of aucdelay2.ctl. After
processing by the nmtemplate and doexpand utilites, it gives a control
stream that is essentially identical to aucdelay2.ctl. It may be eas-
ier to look at the more compact code in aucdelay3.ctl.
The maximum number of records in any individual record is set using
nmtemplate. DELAYDATA contains 16 records. Commands to run aucde-
lay3.ctl are:
nmtemplate aucdelay3.ctl temp1 maxrecs=16
doexpand temp1 temp2
nmfe75 temp2 aucdelay3.res
(The most up-to-date doexpand and nmtemplate utilities are provided in
https://nonmem.iconplc.com/utilities/)
Like aucdelay2S.ctl, aucdelay3.ctl uses DELAYDATA, and ignores records
with CMT=4. An $INFN block is used to save all values of TIME in an
array TIMES(i). There is no need for an explicit NEXTT variable
because all values of TIME are available in this array. Note that the
TIMES array contains all TIME values in the data set, which typically
contains more than one subject's data. In the abbreviated code,
reserved variable IRECIDX gives the position in TIMES of the start of
the current individual record. (LIREC gives the number of records in
the current individual record and is needed if different subjects have
different numbers of data records so that maxrecs is too large for
some subjects.)
At the start of an individual record, a value of MTIME(i) is computed
for every data record as TIMES(KREC)-TDY. (KREC uses reserved vari-
able IRECIDX, the first record in the data set, and hence in TIMES,
for the current individual record.) Some values of MTIME(i) may be
<=0; such values are ignored by PREDPP.
When PK is called at a model event time MTIME(i) the value of A(3) is
stored as SAVEAi. There may be several calls to PK for a given value
of TSTATE and MTIME(i) but SAVEAi must be assigned on only the first
of these calls because only at this call will A(3) have derivaties
with respect to ETA (if MTIME depends on ETA). This is why MTIME(i)=0
is necessary after the value of A(3) (and its ETA derivatives, if any)
has been saved.
REMARKS
1. RECURSIVE CODE
Random variables are saved for future use by using with recursive
code, such as in aucdelay2.ctl:
IF (NEWIND.LE.1.AND.1<=LIREC) THEN
MTIME(1)=TIMES(KREC)-TDY ; SET NEW VALUE OF MTIME (TDY MAY DEPEND ON ETA)
SAVEA1=0 ; ERASE OLD VALUE FROM PREVIOUS PASS
ELSE
MTIME(1)=MTIME(1)
SAVEA1=SAVEA1
ENDIF
It is instructive to note that there is one warning that may be
ignored:
(WARNING 3) THERE MAY BE AN ERROR IN THE ABBREVIATED CODE. THE
FOLLOWING
ONE OR MORE RANDOM VARIABLES ARE DEFINED WITH "IF" STATEMENTS
THAT DO NOT
PROVIDE DEFINITIONS FOR BOTH THE "THEN" AND "ELSE" CASES. IF ALL
CONDITIONS FAIL, THE VALUES OF THESE VARIABLES WILL BE ZERO.
SA
The code in $ERROR is such that SA is always given a value and SA
never defaults to zero. Recursive code for SA could be used to
avoid this warning if desired:
IF ([I]<=LIREC.AND.TIME==TIMES(KREC)) THEN
SA=SAVEA[I]
ELSE
SA=SA
ENDIF
2. SAVEA variables
The variables SAVEAi save the value of compartment 3 (AUC) at
TIME(i)-TDY. They must be assigned individually because there is
no true array feature in NMTRAN for random variables.
3. SIMULATION
The control streams can be used to simulate values of DV. There
will be a difference between simulated values from aucdelay1.ctl
vs. aucdelay2.ctl (2 and 3 are the same). The extra record in
aucdelay1.ctl's FDATA causes an extra call to SIMEPS (even though
it is a dose record), so the simulated etas and eps are not iden-
tical. To have identical simulated values, then do not use
$DATA DELAYDATA IGNORE= IGNORE=(CMT==4)
Instead, use a version of DELAYDATA in which the second record
has EVID=2 ("other-type event record"), so that it is physically
present but does nothing. This adds another record to the data
set. Instead of
nmtemplate aucdelay3.ctl temp1 maxrecs=16
the value maxrecs=17 should be used.
All files may be found in the examples directory:
DELAYDATA
aucdelay1S.ctl
aucdelay2S.ctl
aucdelay1.ctl
aucdelay2.ctl
aucdelay3.ctlREFERENCES: Guide Introduction_7
44 MULTIPLE_DOSE EXAMPLE
This example is based on work described in Fattinger K. et al., Netilmicin in the neonate - population pharmacokinetic analysis and dosing recommendations, Clinical Pharmacology and Therapeutics, 1991, Jul, V50; N1:55-65. It illustrates a data file with multiple doses of different kinds. A FORTRAN format statement is needed to describe the data file because each data record spans two FORTRAN records and also because blanks are sometimes used for null data items. Explanation of data items: TIME hours KG bodyweight(kg) LENG length(cm) SEX 0 for female and 1 for male GAGE gestational age at birth AGE actual age (in days) AGE1 age at start of therapy IMIV im or iv dose (1.0 is iv and 0.0 is im) CREA creatinine value TAD time after last dose (h) All doses are given as zero order infusions. The RATE data item is -2, indicating that duration D1 is modeled in the $PK block. There are two types of doses: IV (intravenous) and IM (intramuscular). The duration of IV doses is 2 minutes. The duration of IM doses is unknown and is given by theta(7). The Estimation Step produces a Model Specification File, which may be used in in a later run; (See $msfi, model specification file).
$PROBLEM NETILMICIN
$INPUT ID MDV TIME AMT RATE CP=DV KG LENG SEX
GAGE AGE AGE1 IMIV CREA TAD
$DATA data (6F7.0/7X,9F7.0)
$SUBROUTINE ADVAN3,TRANS4
$PK
TVCL1=THETA(1)*KG
GAGEN=GAGE/39.0
PAN=(AGE+1.0)/7.0
IF(GAGE.LT.39.0) THEN
TVCL=TVCL1*(GAGEN**THETA(3))*(PAN**THETA(2))
ELSE
TVCL=TVCL1*(PAN**THETA(2))
ENDIF
CL=TVCL*EXP(ETA(1))
TVV1=THETA(4)*KG
V1=TVV1*EXP(ETA(2))
Q=THETA(5)
V2=THETA(6)*KG
S1=V1
IF(IMIV.GE.1.0) THEN
D1=.0333333 ; 2./60.
ELSE
D1=THETA(7)
ENDIF
$ERROR
Y=F*(1+EPS(1))
$THETA (0,0.07) (0,0.1) (0,1.1) (0,0.4) (0,0.07) (0,0.8) (0,.033)
$OMEGA 0.02 0.02
$SIGMA 0.02
$ESTIMATION PRINT=5 MSFO=msfo1
$COVARIANCE
$TABLE MDV ID KG LENG SEX GAGE AGE TAD BY MDV ID NOPRINT FILE=NTABLE.DAT
$SCATTER PRED VS CP UNIT
Data for one subject follow. This subject received both IM and IV
doses.
6.0 1.0 0.0 13.0 -2.0
6.0 2.590 47.0 0.0 36.0 0.0 0.0 1.0 56.0 0.0
6.0 1.0 24.0 10.4 -2.0
6.0 2.590 47.0 0.0 36.0 1.0 0.0 1.0 56.0 0.0
6.0 0.0 26.0 0.0 0.0 6.95
6.0 2.590 47.0 0.0 36.0 1.0 0.0 1.0 56.0 2.0 0.0
6.0 0.0 32.0 0.0 0.0 3.45
6.0 2.590 47.0 0.0 36.0 1.0 0.0 1.0 56.0 8.0 0.0
6.0 0.0 40.0 0.0 0.0 1.8
6.0 2.590 47.0 0.0 36.0 2.0 0.0 1.0 56.0 16.0 0.0
6.0 1.0 48.0 10.4 -2.0
6.0 2.590 47.0 0.0 36.0 2.0 0.0 0.0 56.0 0.0
6.0 1.0 72.0 10.4 -2.0
6.0 2.590 47.0 0.0 36.0 3.0 0.0 0.0 56.0 0.0
6.0 1.0 96.0 10.4 -2.0
6.0 2.590 47.0 0.0 36.0 4.0 0.0 0.0 56.0 0.0
6.0 1.0 120.0 10.4 -2.0
6.0 2.590 47.0 0.0 36.0 5.0 0.0 0.0 56.0 0.0
6.0 0.0 122.0 0.0 0.0 5.9
6.0 2.590 47.0 0.0 36.0 5.0 0.0 0.0 56.0 2.0 0.0
6.0 0.0 128.0 0.0 0.0 2.1
6.0 2.590 47.0 0.0 36.0 5.0 0.0 0.0 56.0 8.0 0.0
6.0 0.0 136.0 0.0 0.0 0.75
6.0 2.590 47.0 0.0 36.0 6.0 0.0 0.0 56.0 16.0 0.0
(See duration, $data).REFERENCES: Guide VI, section V.D; Guide V, section 6.8 , 6.10
45 NONPARAMETRIC EXAMPLE
DISCUSSION:
Three examples are given of the NONPARAMETRIC step in NONMEM. In examples 1 and 2, it is assumed that there only two etas in the problem. With more than two etas, insert the appropriate additional code.
Example 1
This example illustrates how the cumulative density can be plotted versus the conditional etas.
The following code is inserted in the $PRED or $PK block.
DN1=CDEN_(1)
DN2=CDEN_(2)The following control records are used. Note that with the $NONPARAMETRIC record, the MARGINALS option is the default and need not be coded explicitly.
$ESTIMATION POSTHOC PRINT=1
$NONPARAMETRIC MARGINALS
$SCAT DN1 VS ETA1
$SCAT DN2 VS ETA2Example 2
This example illustrates how both the conditional etas and the conditional nonparametric etas (CNPE) can be obtained, and how they can be compared. (See comact).
The following code is inserted in the $PRED or $PK block:
;CONDITIONAL PARAMETRIC ETAS
IF (COMACT.EQ.2) ET1A=ETA(1)
IF (COMACT.EQ.2) ET2A=ETA(2)
;CONDITIONAL NONPARAMETRIC ETAS (CNPE)
IF (COMACT.EQ.3) ET1=ETA(1)
IF (COMACT.EQ.3) ET2=ETA(2)The following control statements may be used. The first two are required; the $SCATTER and $TABLE statements are examples of what can be used.
$ESTIMATION POSTHOC
$NONP ETAS; Examine distribution of nonparametric etas $SCAT ET1 VS ID $SCAT ET2 VS ID ; Compare parametric and nonparametric etas $SCAT ET1A VS ET1 $SCAT ET2A VS ET2 $TABLE ID ET1A ET2A ET1 ET2 NOPRINT
Example 3
This example shows the specific statements needed to simulate from a nonparametric distribution using a Model Specification file generated from a previous analysis using $NONPARAMETRIC.
; Model specification file from a previous analysis using $NONPARAMETRIC $MSFI nonpfit.msf
; Second seed will be used for simulation using the nonparametric density ; Must include TRUE=FINAL, otherwise initial instead of final parameter ; estimates will be used.
$SIMULATION (3432192) (6348371 NONPARAMETRIC) TRUE=FINALREFERENCES: None
46 NWPRI EXAMPLES
These examples illustrate the use of the NONMEM utility routine NWPRI by the NM-TRAN $PRIOR record, and also by an equivalent user-supplied routine PRIOR. For a discussion of PRIOR: (See prior). For a discussion of the NWPRI routine and in particular, a discussion of the
$THETA, $OMEGA, $SIGMA records: (See nwpri).
Example 1
Use of PRIOR and NWPRI and a control stream for single-subject data:
This example obtains parameter estimates from single subject data,
using a two-compartment PK model, and it incorporates a prior of mul-
tivariate normal form for all of the THETA vector.
$PROB SINGLE SUBJECT DATA WITH PRIOR ON THETA
$INPUT ID DOSE TIME DV WT
$DATA data1
$PRIOR NWPRI NTHETA=3,NETA=1,NTHP=3,NPEXP=1
$PRED
;
; THETA(1)=MEAN ABSORPTION RATE CONSTANT
; THETA(2)=MEAN ELIM. RATE CONSTANT (1/HR)
; THETA(3)=CLEARANCE (LITERS/HR)
; DOSE=WEIGHT-ADJUSTED DOSE (MG/KG)
;
IF (NEWIND.EQ.0) AMT=DOSE*WT
T1=THETA(1)
T2=THETA(2)
T3=THETA(3)
A=EXP(-T2*TIME)
B=EXP(-T1*TIME)
C=T1-T2
D=A-B
E=T3*C
Y=AMT*T1*T2/E*D+ETA(1)
$THETA (.4,1.7,7) (.025,.102,.4) (.3,3,30)
;mode for prior with THETA:
$THETA 2.77 FIX .0781 FIX 2.63 FIX
$OMEGA .388
;var-cov for prior with THETA:
$OMEGA BLOCK (3) 5.55 .00524 .00024 -.128 .00911 .515 FIX
$EST
Instead of a $PRIOR record, the following may be used:
$SUBROUTINE PRIOR=prior
The prior routine is as follows:
SUBROUTINE PRIOR (ICALL,CNT)
USE SIZES, ONLY: DPSIZE,ISIZE
REAL(KIND=DPSIZE) :: CNT
REAL(KIND=DPSIZE) :: PLEV
INTEGER(KIND=ISIZE) :: ICALL
NTHETA=3
NETA=1
NTHP=3
NPEXP=1
PLEV=0.
ITYP=0
NSAM=0
ISS=0
NEPS=0
NETP=0
NEPP=0
CALL NWPRI(NTHETA,NETA,NEPS,NTHP,NETP,NEPP,NPEXP,ITYP,PLEV, &
NSAM,ISS,CNT)
RETURN
END
Example 2
First example of use of PRIOR and NWPRI and a control stream for popu-
lation data:
This example obtains parameter estimates from population data, using a
two-compartment PK model, and it incorporates a prior of multivariate
normal/inverse Wishart form for all of the THETA vector and OMEGA
matrix.
$PROB POPULATION DATA WITH PRIOR ON THETA AND OMEGA
$INPUT ID DOSE TIME DV WT
$DATA data2
$PRIOR NWPRI NTHETA=3 NETA=3 NTHP=3 NETP=3 NPEXP=1
$PRED
;
; THETA(1)=MEAN ABSORPTION RATE CONSTANT - MEAN ELIM. RATE CONSTANT (l/HR)
; THETA(2)=MEAN ELIM. RATE CONSTANT (1/HR)
; THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
; DOSE=WEIGHT-ADJUSTED DOSE (MG/KG)
;
IF (NEWIND.NE.2) THEN
AMT=DOSE*WT
W=WT
ENDIF
T0=THETA(1)*EXP(ETA(1))
T2=THETA(2)*EXP(ETA(2))
T1=T2+T0
T3=THETA(3)*W*EXP(ETA(3))
A=EXP(-T2*TIME)
B=EXP(-T1*TIME)
C=T1-T2
D=A-B
E=T3*C
Y=AMT*T1*T2/E*D+EPS(1)
$THETA (0,4,5) (0,.09,.5) (.004,.01,.9)
;mode for prior on THETA:
$THETA 3 FIX .08 FIX .04 FIX
;df for prior on OMEGA:
$THETA 12 FIX
$OMEGA BLOCK (3) .7 .04 .05 .02 .06 .08
;var-cov for prior on THETA:
$OMEGA BLOCK (3) .494 .00207 .0000847 .000692 .0000471 .0000292 FIX
;mode for prior on OMEGA:
$OMEGA BLOCK (3) .7 .04 .05 .02 .06 .08 FIX
$SIGMA .4
$EST
Instead of a $PRIOR record, the following may be used:
$SUBROUTINE PRIOR=prior
The prior routine is as follows:
SUBROUTINE PRIOR (ICALL,CNT)
USE SIZES, ONLY: DPSIZE,ISIZE
REAL(KIND=DPSIZE) :: CNT
REAL(KIND=DPSIZE) :: PLEV
INTEGER(KIND=ISIZE) :: ICALL
NTHETA=3
NETA=3
NTHP=3
NETP=3
NPEXP=1
PLEV=0.
ITYP=0
NSAM=0
ISS=0
NEPS=0
NEPP=0
CALL NWPRI(NTHETA,NETA,NEPS,NTHP,NETP,NEPP,NPEXP,ITYP,PLEV, &
NSAM,ISS,CNT)
RETURN
END
Example 3
Second example of use of PRIOR and NWPRI and a control stream for pop-
ulation data:
This example simulates the THETA and OMEGA parameters of the same two-
compartment population PK model, as well as the population data from
this model. During problem finalization, the entire post-simulation
theta vector and omega array in common NMPR16 are written to external
files.
$PROB POPULATION DATA WITH PRIOR ON THETA AND OMEGA
$INPUT ID DOSE TIME DV WT
$DATA data2
$PRIOR NWPRI
NTHETA=3 NETA=3 NTHP=3 NETP=3 NPEXP=1 ITYP=0 PLEV=.9999
$PRED
;
; THETA(1)=MEAN ABSORPTION RATE CONSTANT - MEAN ELIM. RATE CONSTANT (l/HR)
; THETA(2)=MEAN ELIM. RATE CONSTANT (1/HR)
; THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
; DOSE=WEIGHT-ADJUSTED DOSE (MG/KG)
;
IF (NEWIND.NE.2) THEN
AMT=DOSE*WT
W=WT
ENDIF
T0=THETA(1)*EXP(ETA(1))
T2=THETA(2)*EXP(ETA(2))
T1=T2+T0
T3=THETA(3)*W*EXP(ETA(3))
A=EXP(-T2*TIME)
B=EXP(-T1*TIME)
C=T1-T2
D=A-B
E=T3*C
Y=AMT*T1*T2/E*D+EPS(1)
IF (ICALL.EQ.3) THEN
WRITE (97,*) THSIMP,THSIMPR
WRITE (98,*) OMSIMP(BLOCK)
ENDIF
$THETA (0,4,5) (0,.09,.5) (.004,.01,.9)
;mode for prior on THETA:
$THETA 3 FIX .08 FIX .04 FIX
;df for prior on OMEGA:
$THETA 12 FIX
$OMEGA BLOCK (3) .7 .04 .05 .02 .06 .08
;var-cov for prior on THETA:
$OMEGA BLOCK (3) .494 .00207 .0000847 .000692 .0000471 .0000292 FIX
;mode for prior on OMEGA:
$OMEGA BLOCK (3) .7 .04 .05 .02 .06 .08 FIX
$SIGMA .4
$SIM (547676) ONLY
Suppose a $TABLE record is added to the above control stream. To
obtain values of DV computed with the simlulated values of eta,
include the following code:
IF (ICALL.EQ.4) THEN
T0=THSIMP(1)*EXP(ETA(1))
T2=THSIMP(2)*EXP(ETA(2))
T3=THSIMP(3)*W*EXP(ETA(3))
ELSE
T0=THETA(1)*EXP(ETA(1))
T2=THETA(2)*EXP(ETA(2))
T3=THETA(3)*W*EXP(ETA(3))
ENDIF
Instead of a $PRIOR record, the following may be used:
$SUBROUTINE PRIOR=prior
The prior routine is as follows:
SUBROUTINE PRIOR (ICALL,CNT)
USE SIZES, ONLY: DPSIZE,ISIZE
REAL(KIND=DPSIZE) :: CNT
REAL(KIND=DPSIZE) :: PLEV
INTEGER(KIND=ISIZE) :: ICALL
PLEV=.9999
NTHETA=3
NETA=3
NTHP=3
NETP=3
NPEXP=1
NEPS=0
NEPP=0
ITYP=0
NSAM=0
ISS=0
CALL NWPRI(NTHETA,NETA,NEPS,NTHP,NETP,NEPP,NPEXP,ITYP,PLEV, &
NSAM,ISS,CNT)
RETURN
ENDREFERENCES: none
47 PHENOBARB EXAMPLE
This example illustrates simple kinetics.
$PROBLEM PHENOBARB SIMPLE MODEL
$INPUT ID TIME AMT WGT APGR DV
$DATA PHENO
$SUBR ADVAN1 TRANS2
$PK
TVCL=THETA(1)
TVV=THETA(2)
CL=TVCL*EXP(ETA(1))
V=TVV*EXP(ETA(2))
S1=V
$THETA (0,.0105) (0,1.05)
$OMEGA .25 .25
$ERROR
W=F
Y=F+W*EPS(1)
IPRED=F ; individual-specific prediction
IRES=DV-IPRED ; individual-specific residual
IWRES=IRES/W ; individual-specific weighted residual
$SIGMA .04
$ESTIMATION SIGDIGITS=4 POSTHOC
$COVARIANCE
$TABLE ID WGT APGR TIME AMT CL V ETA1 ETA2 IPRED IWRES
NOPRINT ONEHEADER FILE=TABLE1
$SCATTER DV VS PRED UNIT
$SCATTER (RES,WRES) VS PRED
$SCATTER (CL,V) VS (WGT,APGR)
$SCATTER IWRES VS IPRED BY MDV ORDZERO
Data for the first individual follow.
1 0. 25.0 1.4 7 .
1 2.0 . 1.4 7 17.3
1 12.5 3.5 1.4 7 .
1 24.5 3.5 1.4 7 .
1 37.0 3.5 1.4 7 .
1 48.0 3.5 1.4 7 .
1 60.5 3.5 1.4 7 .
1 72.5 3.5 1.4 7 .
1 85.3 3.5 1.4 7 .
1 96.5 3.5 1.4 7 .
1 108.5 3.5 1.4 7 .
1 112.5 . 1.4 7 31.0REFERENCES: Guide III, section VI.2; Guide V, section 2.3
48 PK_PD_SEQUENTIAL_1 EXAMPLE
This is an example for analysis of direct PD using previously estimated individuals’ PK parameters. (See PK_PD_sequential_2 example).
$PROBLEM EMAX MODEL APPLIED TO EFFECT SITE CONCENTRATION
; NOTE: INDIVIDUAL K, V ESTIMATED FROM PREVIOUS PK FIT
; AND RECORDED IN DATA
;
; THE DATA FILE CONTAINS *ONLY* EFFECT OBSERVATIONS,
; CMT = 2 NEEDED TO TELL PREDPP TO SET F = CE
$DATA data
$INPUT ID TIME DV AMT=DOSE CMT KK V1
; |---|
; KK AND V1 ARE ESTIMATES OF K AND V
$SUBROUTINES ADVAN3
$PK
K = KK
K12= .001*K ; TRIVIAL LOSS TO EFFECT COMPT.
K21= THETA(1)*EXP(ETA(1))
EMAX=THETA(2)*EXP(ETA(2))
C50= THETA(3)*EXP(ETA(3))
S2= V1*K12/K21 ; SO THAT CESS = CPSS
$ERROR
Y =EMAX*F/(C50+F)*EXP(ERR(1))
; ADD APPROPRIATE $THETA, $OMEGA, $SIGMA, $ESTIMATION, ETC.
If, moreover, an absorption lag (ALAG1) parameter is modeled in $PK
with parameters estimated from the PD data only, then this introduces
a lag in the dose - Ce relationship. Since this lag was not present
in the dose - Cp relationship, it may be interpreted as a lag in the
Cp - Ce relationship.REFERENCES: Guide IV, section V.C.6; Guide VI, section IV.B.2 , VII.C.3
49 PK_PD_SEQUENTIAL_2 EXAMPLE
This is an example for analysis of direct PD using previously estimated population PK parameters. (See PK_PD_sequential_1 example).
$PROBLEM EMAX MODEL APPLIED TO EFFECT SITE CONCENTRATION
; THE DATA FILE CONTAINS *ONLY* EFFECT OBSERVATIONS,
; POP PK PARAMS ESTIMATED IN PREVIOUS RUN & FIXED HERE.
; CMT = 2 NEEDED TO TELL PREDPP TO SET F = CE
;
$DATA data
$INPUT ID TIME DV AMT=DOSE CMT
$SUBROUTINES ADVAN3
$PK
K = THETA(1)*EXP(ETA(1))
K12= .001*K ; TRIVIAL LOSS TO EFFECT COMPT.
K21= THETA(2)*EXP(ETA(2))
V1= THETA(3)*EXP(ETA(3))
EMAX=THETA(4)*EXP(ETA(4))
C50= THETA(5)*EXP(ETA(5))
S2= V1*K12/K21 ; SO THAT CESS = CPSS
$ERROR
Y= EMAX*F/(C50+F)*EXP(ERR(1))
$THETA
.124 FIXED ; K - FIXED TO PREV EST
.1 ; KEO - ESTIMATED
28.6 FIXED ; V1 - FIXED TO PREV EST
100 ; EMAX- ESTIMATED
50 ; C50 - ESTIMATED
$OMEGA
.113 FIXED ; OMEGA(K) - FIXED TO PREV EST
.04
.627 FIXED ; OMEGA(V1)- FIXED TO PREV EST
.04
.04
$SIGMA .16
; ADD $ESTIMATION, ETC, AS DESIRED.
If, moreover, an absorption lag (ALAG1) parameter is modeled in $PK
with parameters estimated from the PD data only, then this introduces
a lag in the dose - Ce relationship. Since this lag was not present
in the dose - Cp relationship, it may be interpreted as a lag in the
Cp - Ce relationship.REFERENCES: Guide IV, section V.C.6; Guide VI, section IV.B.2 , VII.C.3
50 PK_PD_SIMULTANEOUS_1 EXAMPLE
This is an example for simultaneous analysis of PK and direct PD. Kinetics are one compartment. (See PK_PD_simultaneous_2, PK_PD_simultaneous_3).
$PROBLEM SIMULT PK/PD - ADVAN3
; THE DATA FILE CONTAINS BOTH CP AND EFFECT OBSERVATIONS.
; WHEN DV IS A CP OBSERVATION, CMT = 1 (OR 0)
; WHEN DV IS AN EFFECT OBSERVATION, CMT =2
$DATA data
$INPUT ID TIME DV AMT=DOSE CMT
$SUBROUTINES ADVAN3
$PK
K = THETA(1)*EXP(ETA(1))
K12= .001*K ; TRIVIAL LOSS TO EFFECT CMPT
K21= THETA(2)*EXP(ETA(2)) ; KEO
S1= THETA(3)*EXP(ETA(3)) ; V1 FOR DRUG
EMAX=THETA(4)*EXP(ETA(4))
C50= THETA(5)*EXP(ETA(5))
S2= S1*K12/K21 ; SO THAT CESS = CPSS
$ERROR
IF (CMT.EQ.2) THEN
Y = EMAX*F/(C50+F)*EXP(ERR(2)) ; DV = EFFECT IF CMT=2
ELSE
Y = F*EXP(ERR(1)) ; DV = CP IF CMT=1 (OR 0)
ENDIF
; ADD APPROPRIATE $THETA, $OMEGA, $SIGMA, $ESTIMATION, ETC.REFERENCES: Guide IV, section V.C.6; Guide VI, section IV.B.2 , VII.C.3
51 PK_PD_SIMULTANEOUS_2 EXAMPLE
This is an example for simultaneous analysis of PK and PD. Kinetics are one-compartment. (See PK_PD_simultaneous_1, PK_PD_simultaneous_3).
$PROBLEM SIMULTANEOUS PK/PD - ADVAN6
; NOTE: EXAMPLE USES 1-COMPT PK;
; MORE COMPLEX COMPARTMENTAL MODEL CAN BE USED
; CMT = 1 (OR 0) FOR DV = CP; CMT = 2 FOR DV = EFFECT
$DATA data
$INPUT ID TIME DV AMT=DOSE CMT
$SUBROUTINE ADVAN6 TOL=5
$MODEL
COMP=(CENTRAL,DEFDOSE,DEFOBS)
COMP=EFFECT
$PK
K10= THETA(1)*EXP(ETA(1))
V1= THETA(2)*EXP(ETA(2))
S1= V1
K20= THETA(3)*EXP(ETA(3))
C50= THETA(4)*EXP(ETA(4))
EMAX=THETA(5)*EXP(ETA(5))
$DES
DADT(1)=-K10*A(1) ;NOTE NO LOSS TO A(2)
DADT(2)= K20*(A(1)/V1-A(2)) ;A(2) IS CE
$ERROR
IF (CMT.EQ.2) THEN
Y = EMAX*F/(C50+F)*EXP(ERR(2)) ; DV = EFFECT IF CMT=2
ELSE
Y = F*EXP(ERR(1)) ; DV = CP IF CMT=1 (OR 0)
ENDIF
; ADD APPROPRIATE $THETA, $OMEGA, $SIGMA, $ESTIMATION, ETC.REFERENCES: Guide IV, section V.C.6 , V.C.7; Guide VI, section IV.B.2 , VII.C.6
52 PK_PD_SIMULTANEOUS_3 EXAMPLE
This is an example for simultaneous analysis of PK and direct PD. Kinetics are two-compartment. Either ADVAN5 or ADVAN7 could be used. (See PK_PD_simultaneous_1, PK_PD_simultaneous_2).
$PROBLEM SIMULTANEOUS PK/PD - ADVAN7
; NOTE: MODEL FOR CP CAN BE MORE COMPLEX BY ADDING CMPTS
;
; THE DATA FILE CONTAINS BOTH CP AND EFFECT OBSERVATIONS.
; WHEN DV IS A CP OBSERVATION, CMT = 1 (OR 0),
; WHEN DV IS AN EFFECT OBSERVATION, CMT =3.
$DATA data
$INPUT ID TIME DV AMT=DOSE CMT
$SUBROUTINES ADVAN7
$MODEL
COMP=(CENTRAL,DEFDOSE,DEFOBS)
COMP=PERIPH
COMP=EFFECT
$PK
K10= THETA(1)*EXP(ETA(1))
K12= THETA(2)*EXP(ETA(2))
K13= .001*K10 ; TRIVIAL LOSS TO EFFECT COMPT
K21= THETA(3)*EXP(ETA(3))
S1= THETA(4)*EXP(ETA(4)) ; V1 FOR DRUG
K30= THETA(5)*EXP(ETA(5)) ; KEO
EMAX=THETA(6)*EXP(ETA(6))
C50= THETA(7)*EXP(ETA(7))
S3=S1*K13/K30 ; SO THAT CESS = CPSS
$ERROR
IF (CMT.EQ.3) THEN
Y =EMAX*F/(C50+F)*EXP(ERR(2)) ; CMT = 3; DV = EFFECT
ELSE
Y = F*EXP(ERR(1)) ; CMT = 1 (OR 0); DV = CP
ENDIF
; ADD APPROPRIATE $THETA, $OMEGA, $SIGMA, $ESTIMATION, ETC.REFERENCES: Guide IV, section V.C.6; Guide VI, section IV.B.2 , VII.C.7
53 PK_PD_SIMULTANEOUS_4 EXAMPLE
This is an example for simultaneous analysis of PK and indirect PD, where the PREDPP PK prediction is modified by $ERROR code.
$PROB WARFARIN
$INPUT ID TIME DV AMT CMT
; DV = DRUG CONC OR EFFECT (PROTHROMBIN TIME)
$DATA data
$SUBROUTINE ADVAN6 TOL=5
$MODEL
COMP=(CENTRAL,DEFDOSE,DEFOBS)
; CAN ADD ADD'L COMPARTMENTS IF PK MULTICOMPARTMENTAL
COMP=(CLOTFAC)
$PK
K10= THETA(1)*EXP(ETA(1)) ;LOSS RATE DRUG
S1= THETA(2)*EXP(ETA(2)) ;VD DRUG
SYNTH=THETA(3)*EXP(ETA(3)) ;SYNTH RATE OF CLOTING FACTOR
LOSS= SYNTH
C50= THETA(4)*EXP(ETA(4))
$DES
DADT(1)=-K10*A(1)
EFF = C50/(C50+A(1)/S1)
DADT(2)= SYNTH*EFF - LOSS*A(2) ;CF LEVEL
$ERROR
PT = THETA(5)/(THETA(6)+F) ; HYPERBOLIC MODEL FOR PRO-TIME
IF (CMT.EQ.2) THEN
Y = PT*EXP(EPS(2)) ; OBS IS PRO-TIME
ELSE
Y = F*EXP(EPS(1)) ; OBS IS DRUG CONCENTRATION
ENDIF
; ADD APPROPRIATE $THETA, $OMEGA, $SIGMA, $ESTIMATION, ETC.
The system should be initialized with a bolus dose of AMT=1 into com-
partment 2 (CLOTFAC) at time zero to establish the initial condition
A(1)=0, A(2)=1 at time zero. This is appropriate because, in the
absence of drug, the system is assumed to be at steady-state
(DADT(2)=0).REFERENCES: Guide IV, section V.C.6 , V.C.7; Guide VI, section IV.B.2 , VII.C.6
54 PLASMA URINE EXAMPLE
Here are several examples of simultaneous modelling of plasma and urine data. The examples are adapted from a model for the analysis of Dextromethorphan data provided by Khaled Mohammed Abduljalil. The first involves data from the parent drug. The second involves both parent and metabolite data. In the control streams of these examples plasma and urine data are simulated. The same control streams could also be used for analysis of the data. MU_ modelling is used for the simulation to maintain consistency with subsequent analysis runs but MU_ modelling has no effect on simulation. If MU_ variables are renamed to avoid the character “_” (e.g., rename MU_1 as MU1), the same control streams may be run with NONMEM VI. EXAMPLE 1 sim_parent_a2.ctl The data dextroparent.dat contains data from only the parent drug. The parent drug is elimated directly to urine, and also via conversion to metabolite, but there are no measurements of metabolite in plasma or urine. The fraction F0 of parent drug in the urine is computed as in Guide V, Chapter 4, Section 3.1 equation (4.6) and Section 3.4.1 equation (4.7). (See Guide V, Chapter 4.3) The L2 data item is not required, but is used because multiple observations at the same TIME are in effect multivariate observations. Values of SIGMA corresponding to unused elements EPS(3) and EPS(4) are defined so that the EPS vector has the same length in both examples. See Remarks below. Although ADVAN2 is used, the same control stream could be implemented with ADVAN5 or ADVAN6 or ADVAN7, with minor changes.
$PROBLEM Parent drug, using ADVAN2
$INPUT ID TIME AMT UVOL DV CMT MDV EVID L2
$DATA dextroparent.dat IGNORE=#
$SUBROUTINES ADVAN2 TRANS1
$PK
K12=THETA(1)
KA=K12
MU_1=LOG(THETA(2))
V2=EXP(MU_1+ETA(1))
MU_2=LOG(THETA(3))
CLP=EXP(MU_2+ETA(2)) ; RENAL CL FOR PARENT
CLB=THETA(4) ; METABOLIC CL FOR METABOLIC
K23=CLP/V2+CLB/V2
K=K23
F0=CLP/(CLP+CLB)
S2=V2
S3=UVOL
$ERROR (EVERY EVENT)
ACMT=ABS(CMT) ; output compt. may have neg. value of CMT
IF(ACMT.EQ.2) Y=F*(1+EPS(1))
IF(ACMT.EQ.3) Y=F*(1+EPS(2))
$THETA
(0.01,0.8,6) ;KA
(0.01,43,1000);V2
(0.0001,20,190);CLP
(0.01,15,90);CLB
$OMEGA
0.05 0.05
$SIGMA
.01 .01
.01 .01 ; eps(3), eps(4) for consistent EPS with sim_metab_a6 example
$SIM (111111) ONLYSIM
$TABLE ID TIME AMT UVOL DV SIMP=PRED CMT MDV
EVID L2 NOAPPEND FILE=sim_parent_a2.tab
A fragment of the data file is as follows:
# ID TIME AMT UVOL DV CMT MDV EVID L2
2.000 0.000 30000.000 0.000 0.000 1.000 1.000 1.000 1
2.000 0.171 0.000 0.000 0.000 3.000 1.000 2.000 1
2.000 2.000 0.000 0.000 0.000 2.000 0.000 0.000 1
2.000 2.000 0.000 93.100 0.000 -3.000 0.000 0.000 1
2.000 2.000 0.000 0.000 0.000 3.000 1.000 2.000 2
2.000 3.150 0.000 0.000 0.000 2.000 0.000 0.000 2
2.000 3.150 0.000 133.600 0.000 -3.000 0.000 0.000 2
There is a dose at TIME=0. At TIME=.171, the urine compartment (which
is the default compartment for output) is turned on. At TIME=2, there
are observations of both plasma and urine. The urine amount is reset
to 0 and the compartment turned back on for the start of a new urine
collection. At TIME 3.15 there are observations of plasma and urine.
In the model, different EPS variables are used for the two types of
observations. In the table, values of PRED (SIMP) are computed with
ETA and EPS equal to 0. The DV values are computed with simulated ETA
and EPS, and would be used as the DV in subsequent analysis runs.
EXAMPLE 2 sim_metab_a6.ctl
The data dextrometab.dat includes measurements of metabolite in
plasma and urine. The model also include predictions of the metabo-
lite in both plasma (CMT=4) and urine (CMT=5; the default compartment
for output). Different compartments are used for parent and metabo-
lite in urine, but the same value of UVOL applies to both. The F0
parameter is not needed because the model for K23 now uses only CLP,
the clearance of the parent drug. No analytic ADVAN model has two
output compartments, so a general ADVAN (ADVAN6 or ADVAN5 or ADVAN7)
must be used.
Note that Compartment 3 is no longer the default for output. Instead
of being computed by PREDPP using mass balance, it is computed explic-
itly by the ADVAN routine itself (e.g., using a differential equation
when ADVAN6 is used). If compartment 3 is defined in $MODEL with
default attributes, e.g.,
COMP=(DEXURIN)
there is an error messsage from PREDPP:
SPECIFIED COMPARTMENT MAY NOT BE TURNED OFF WITH AN OBSERVATION RECORD
Instead, the compartment may be defined as an output-type compartment:
COMP=(DEXURIN INITIALOFF NODOSE)
Now observations with CMT=-3 are permitted and no change to the origi-
nal data records is needed.
$PROBLEM Parent drug and metabolite, using ADVAN6
$INPUT ID TIME AMT UVOL DV CMT MDV EVID L2
$DATA dextrometab.dat IGNORE=#
$SUBROUTINES ADVAN6 TRANS1 TOL=4
$MODEL
COMP=(DEPOT)
COMP=(PLASMA DEFOBS) ; PARENT IN PLASMA
COMP=(DEXURIN INITIALOFF NODOSE) ;PARENT IN URINE
COMP=METAB ;METABOLITE IN PLASMA
$PK
K12=THETA(1)
MU_1=LOG(THETA(2))
V2=EXP(MU_1+ETA(1))
MU_2=LOG(THETA(3))
CLP=EXP(MU_2+ETA(2)) ; RENAL CL FOR PARENT
CLB=THETA(4) ; METABOLIC CL FOR METABOLIC
CLMR=THETA(5) ; RENAL CL FOR METABOLITE
V4=1
K24=CLB/V2
K23=CLP/V2
; F0=CLP/(CLP+CLB) Omit F0 because parent and
; metab have different urine compts.
K45=CLMR/V4
S2=V2
S4=V4
S3=UVOL
S5=UVOL
$ERROR (EVERY EVENT)
ACMT=ABS(CMT) ; output compart. may have neg. value of CMT
IF(ACMT.EQ.2) Y=F*(1+EPS(1))
IF(ACMT.EQ.3) Y=F*(1+EPS(2))
IF(ACMT.EQ.4) Y=F*(1+EPS(3))
IF(ACMT.EQ.5) Y=F*(1+EPS(4))
$DES
DADT(1)=-K12*A(1)
DADT(2)=K12*A(1)-K23*A(2)-K24*A(2)
DADT(3)=K23*A(2)
DADT(4)=K24*A(2)-K45*A(4)
$THETA
(0.01,0.8,6) ;KA
(0.01,43,1000);V2
(0.0001,20,190);CLP
(0.01,15,90);CLB
(0.0001,5,90);CLMR
$OMEGA
0.05 0.05
$SIGMA
.01 .01 .01 .01
$SIM (111111) ONLYSIM
$TABLE ID TIME AMT UVOL DV SIMP=PRED CMT MDV EVID L2 NOAPPEND
FILE=sim_metab_a6.tab
A fragment of the data dextrometab.dat is:
# ID TIME AMT UVOL DV CMT MDV EVID L2
2.000 0.000 30000.000 0.000 0.000 1.000 1.000 1.000 1
2.000 0.171 0.000 0.000 0.000 3.000 1.000 2.000 1
2.000 0.171 0.000 0.000 0.000 5.000 1.000 2.000 1
2.000 2.000 0.000 0.000 0.000 2.000 0.000 0.000 1
2.000 2.000 0.000 0.000 0.000 4.000 0.000 0.000 1
2.000 2.000 0.000 93.100 0.000 -3.000 0.000 0.000 1
2.000 2.000 0.000 93.100 0.000 -5.000 0.000 0.000 1
2.000 2.000 0.000 0.000 0.000 3.000 1.000 2.000 2
2.000 2.000 0.000 0.000 0.000 5.000 1.000 2.000 2
2.000 3.150 0.000 0.000 0.000 2.000 0.000 0.000 2
2.000 3.150 0.000 0.000 0.000 4.000 0.000 0.000 2
2.000 3.150 0.000 133.600 0.000 -3.000 0.000 0.000 2
2.000 3.150 0.000 133.600 0.000 -5.000 0.000 0.000 2
Observations of metabolite in plasma (CMT=4) and metabolite in urine
(CMT=5) are present at the same values of TIME as observations of par-
ent in plasma (CMT=2) and parent in urine (CMT=3), respectively, and
have the same values of L2.
Remarks:
After the control streams are run, the predictions SIMP for compart-
ments 2 and 3 are the same in the two table files sim_parent_a2.tab
and sim_metab_a6.tab, as expected, because they are computed with
EPS=0, i.e., they are not simulated values. Because of the use of L2
and the same-length EPS vector in both examples, the simulated DV's
for compartments 2 and 3 are also the same. This is useful for illus-
trative and debugging purposes, but is not usually the case when the
EPS structure is different, or the L2 data item is not used.
To summarize:
With L2, during simulation SIMEPS generates a new set of values for
EPS only when L2 changes value.
Without L2, during simulation SIMEPS generates a new set of values for
EPS with every record, even when MDV is 1 or the PREDPP item EVID is
not 0.
(See simeps).
All control streams and data files are found in the NONMEM 7.4 exam-
ples\dextro directory.
dextroparent.dat (also displayed here)
sim_parent_a2.ctl (also displayed here)
sim_parent_a5.ctl
sim_parent_a6.ctl
dextrometab.dat (also displayed here)
sim_metab_a5.ctl
sim_metab_a6.ctl (also displayed here)REFERENCES: None
55 RECORDS=ID EXAMPLE
These examples illustrate the use of the RECORDS=ID option of the
$DATA record. (See $data).
Example 1 - Data checkout
During a checkout phase, one might want to see only the predictions
for the first individual record, and therefore, use a data set com-
prised only of the data records therein. With RECS=ID, it is not nec-
essary to count the number of records in the fist individual record.
$DATA filename RECS=ID
Example 2 - Data Analysis
One might want to analyze each subject's data independently of the
others. One way this can be done is by creating a single control
stream with as many problems as there are individual records in the
data, and including one each of the following series of $DATA records
in the series of problem specifications:
$DATA filename RECS=n1
$DATA filename RECS=n2 NOREWIND
$DATA filename RECS=n3 NOREWIND
etc.
where "nk" is the number of records in the kth individual record. It
would be necessary to count the exact number of records. Now, these
records can simply be coded
$DATA filename RECS=ID
$DATA filename RECS=ID NOREWIND
$DATA filename RECS=ID NOREWIND
etc.
Example 3 - Use of the INCLUDE record
With the INCLUDE record, the filename may now be followed by an inte-
ger n, whose default value is 1. Then NM-TRAN reads n copies of the
named file.
In the example above if the problem specification for all subjects
after the first are completely identical, a compact way of writing the
control stream is possible. The control stream would contain the
problem specification for the first individual, including the $SUBROU-
TINES record, abbreviated code, and the $DATA record for the first
subject, and would end with (if there are e.g. 12 subjects in all)
INCLUDE ctlfile2 11
The file ctlfile2 would contain one problem specification (including
$PROBLEM) for one subject, with no $SUBROUTINES record or abbreviated
code, but would include
$DATA filename RECS=ID NOREWIND
A fully worked-out example is supplied as recid2.exa (See recid2.exa).
An alternate approach for single-subject data is available with NONMEM
7. OMEGA diagonal values are fixed to a special value 1.0E+06.
See Guide Introduction_7 "Single-Subject Analysis using Population
with Unconstrained ETAs ".REFERENCES: Guide Introduction_7
56 RECORDS=ID EXAMPLE 2
This fully-worked out example shows how the RECORDS=ID option of $DATA and the INCLUDE statement may be used to analyze the Theophyline data set THEOPP as separate individuals. The NONMEM control file is as follows.
$PROB THEOPHYLLINE POPULATION DATA; Analysis of Individuals
; Modification of CONTROL5 control steam
$INPUT ID DOSE=AMT TIME CP=DV WT
$DATA THEOPP RECS=ID
;RECS=ID: Data set will be read until ID changes or end-of-file
$SUBROUTINES ADVAN2
$PK
;THETA(1)=MEAN ABSORPTION RATE CONSTANT (1/HR)
;THETA(2)=MEAN ELIMINATION RATE CONSTANT (1/HR)
;THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
;SCALING PARAMETER=VOLUME/WT SINCE DOSE IS WEIGHT-ADJUSTED
CALLFL=1
KA=THETA(1)
K=THETA(2)
CL=THETA(3)*WT
SC=CL/K/WT
$THETA (.1,3,5) (.008,.08,.5) (.004,.04,.9)
$OMEGA .4
;For single subject data OMEGA is residual variance.
$ERROR
Y=F+ERR(1)
;ERR must be used instead of EPS.
$EST MAXEVAL=450 PRINT=5
$COV SPECIAL
;SPECIAL is required to obtain the variance-covariance matrix for
;single-subject data.
$TABLE ID DOSE WT TIME NOPRINT ONEHEADER FILE=indest.tab
$TABLE ID KA K CL SC NOPRINT FIRSTONLY ONEHEADER NOAPPEND
FILE=indest.par
INCLUDE indest.txt 11
; INCLUDE: Inserts copies of the file named indest.txt for each
; additional individual.
The contents of file indest.txt are:
$PROB THEOPHYLLINE POPULATION DATA; Analysis of Individuals
$INPUT ID DOSE=AMT TIME CP=DV WT
$DATA THEOPP RECS=ID NOREWIND
;NOREWIND: data set will be read starting after the previous individual
$THETA (.1,3,5) (.008,.08,.5) (.004,.04,.9)
$OMEGA .4
;For single subject data OMEGA is residual variance
$EST MAXEVAL=450 PRINT=5
$COV SPECIAL
;SPECIAL is required to obtain the variance-covariance matrix for
;single-subject data
$TABLE ID DOSE WT TIME NOPRINT ONEHEADER FORWARD FILE=indest.tab
$TABLE ID KA K CL SC NOPRINT FIRSTONLY ONEHEADER FORWARD NOAPPEND
FILE=indest.par
An alternate approach for single-subject data is available with NONMEM
7. OMEGA diagonal values are fixed to a special value 1.0E+06.
See Guide Introduction_7 "Single-Subject Analysis using Population
with Unconstrained ETAs ".REFERENCES: Guide Introduction_7
57 REPETITION_1 EXAMPLE
This example illustrates the use of the repetition feature, and in particular, the use of the variables RPTI and RPTO. For a discussion of these reserved variables: (See Repetition Variables). We suppose that a compartment amount A2 at time t1 is given by a convolution, i.e. an integral with respect to s, of the product D(t1-s) and r(s), where D is a unit distribution function, and r is an input function. The integral is taken from 0 to t1. The amount in question can be computed by integrating the differential equation dA2/ds=D(t1-s)r(s) from time 0 to time t1. For A2 at a later time t2, one integrates a different differential equation: dA2/ds=D(t2-s)r(s) from time 0 to time t2. This equation differs from the first in so far as the constant t2, rather than t1, is used in D. Ordinarily, PREDPP will use a system of differential equations dA2/ds=f(A2,t) to obtain the value of a vector A2 of compartment amounts at times t1 and t2, by first integrating the system to time t1, and then by integrating the same or different system from time t1 to time t2. This scheme cannot be used with the two equations described above, because the second equation - with the constant t2 - applies throughout the time interval 0 to t1, as well as throughout the interval t1 to t2. PREDPP requires the event records to be ordered chronologically in the data set, and with its ordinary usage with differential equations (described above), these records are processed sequentially, starting with the first event record (which, for the purposes of this discussion, can be thought to have time 0). The records are passed sequentially to the PK routine. The differential equations used for the convolution requires something special. During a regular sequential pass of the event records to the PK routine, when a record R with time t - at which compartment amounts are desired - is initially processed and passed to the PK routine for the first time, the routine should set the output variable RPTO to -1. This signals to NONMEM that the regular pass should be interrupted and a second pass of the records should be initiated, starting with the first event record once again and ending with record R (thus the pass of these records is “repeated”). Only at the end of the second pass, after all event records up to and including record R have once again been passed, will the value for A2(t) be established. The amounts at times on preceeding records have already been established, and they remain unchanged as the pass of these records is repeated. After R has been passed for the second time, the first pass will be continued, the next event record will be passed for the first time, and RPTO can again be set to -1. Whenever an event record is passed for the first time, f(s) should be set to D(t-s)r(s), where t is the time on the record, so that during the second pass, the differential equation that is used is the desired one. When an event record is passed, if the input variable to the PK routine, RPTI, has the value 0, this signals that the record is not being passed as a part of a repeated pass, and so this can be used as a signal that the record is being passed for the first time. Suppose the kinetic system is a two-compartment system where drug absorption from a depot compartment into a sampling compartment occurs in a first-order manner, and drug elimination from the sampling compartment occurs according to a Weibull unit disposition function D(s)=exp(-K(s**c)). Suppose too that the superposition principle is assumed to hold. So drug amount is given by a convolution. Suppose the event records for an individual are simply structured, e.g. # TIME AMT DV EVID .0 10 . 1 rec 1 .25 . 2.44 0 rec 2 .5 . 5.24 0 rec 3 .75 . 5.57 0 rec 4 1 . 5.85 0 rec 5 Then the control stream should look like:
$PROB
$INPT TIME AMT DV EVID
$DATA data
$SUB ADVAN6 TOL=5
$MODEL COMP=(DEPOT DEFDOSE) COMP=(CENTRAL DEFOBS)
$INFN
IF (ICALL.EQ.0) RPTO=1 ;enables use of repetition feature
$PK
KA=THETA(1)*EXP(ETA(1))
K =THETA(2)*EXP(ETA(2))
C =THETA(3)
V =THETA(4)*EXP(ETA(4))
S2=V
IF (RPTI.EQ.0) TI=TIME
IF (NEWIND.EQ.2) RPTO=-1
$DES
DADT(1)=-KA*A(1)
D=EXP(-K*(TI-T)**C)
DADT(2)=D*KA*A(1)
etc.
Note that when TIME= 0.25, RPTO need not have been set to -1; after
the second event record has been passed for the first time, the kinet-
ics are advanced from time 0 to time 0.25 using TI=0.25.
(1) The above is supplied as repeat1.ctl and repeat1.dat in the exam-
ples directory
The code
IF (ICALL.EQ.0) RPTO=1 ;enables use of repetition feature
is needed to inform NONMEM that the repetition feature will be used.
The code
IF (NEWIND.EQ.2) RPTO=-1
initiates a pass starting with the first event record, which is the
default with RPTO=-1.
(2) Another example, repeat1t.ctl, demonstates that a repeated pass
may start with any designated record, not necessarily the first
record. The statement
IF (NEWIND.NE.2) RPTO=1
is used to indicate explicitly that the first record is a repetition
base.
(3) Example repeat1s.ctl uses the RPT_ data item for the same purpose.
This is repeat1s.ctl and repeat1s.dat in the examples directory. A
non-zero value of RPT_ marks a record as a repetition base. A state-
ment setting RPTO=-n (where n is the value of the RPT_ data item)
requests repetition starting with the marked record. When RPT_ is
used, it automatically enables use of the repetition feature and the
$INFN block may be omitted.
The data is augmented with the RPT_ item, as in this fragment:
# TIME AMT DV EVID RPT_
1 .0 10 . 1 1 rec 1
1 .25 . 2.44 0 . rec 2
and the control stream contains
$INPT ID TIME AMT DV EVID RPT_
(4) A different example of the use of RPT_ is given in repeatf.ctl and
repeatf.dat. Reserved variables RPTO=n and RPTO=-n are intended to
mark and repeat a sequence of more than 1 data records. If only one
record is to be repeated, RPT_ must be used. The example computes the
factorial of a value supplied as KRPT in the data record. A fragment
of the data is
ID TIME DV MDV RPT_ KRPT
1 1 1 0 . .
1 2 2 0 1 4
The code sets
IF (RPT_.EQ.1) RPTO=-1
RPTON=KRPT
NONMEM calls PRED KRPT+1 times for each record with RPT_=1. This
allows the factorial to be computed. This is intended as an illustra-
tion of what can be done using RPT_ data item. With NONMEM 7.3, fac-
torials can be computed more quickly using the DOWHILE feature and the
GAMLN function.REFERENCES: none
58 REPETITION_2 EXAMPLE
This example illustrates the use of the repetition feature, and in particular, the use of the PRED output control flag PRDFL. For a discussion of the reserved variables RPTO, RPTI and PRDFL: (See Repetition Variables). This example computes a variable RS to be output in a table, defined as the difference between the observation in the data record and the average observation in the individual record. User-defined tabled quantities are computed during a copying pass through the data set - signaled by COMACT=1 (See comact). In this example, during the copying pass, the data records in each individual record are actually passed twice; the average is computed during a first pass through the individual record - signaled by RPTI=0, and the difference is computed during a second (repeated) pass - signaled by RPTO=1.
$PRED
Y=...
IF (COMACT.EQ.0) PRDFL=1
IF (ICALL.EQ.0) THEN
RPTO=1 ;initiates use of repetition feature
PRDFL=1 ;initiates use of the PRED output control flag
N=0
AV=0
ENDIF
IF (COMACT.EQ.1) THEN
IF (RPTI.EQ.0.AND.NEWIND.LE.1) THEN
N=0
AV=0
ENDIF
IF (RPTI.EQ.0.AND.MDV.EQ.0) THEN
N=N+1
AV=AV+DV
ENDIF
IF (RPTI.EQ.1) THEN
RS=0
PRDFL=1
ENDIF
IF (RPTI.EQ.1.AND.MDV.EQ.0) RS=DV-AV
IF (NDREC.EQ.LIREC) THEN
RPTO=-1 ;RPTO=-1 is a signal initiating
;a second pass. As RPTON isn't also set,
;there are no further passes.
AV=AV/N
ENDIF
ENDIF
$TABLE RS NOAPPEND
The variable NDREC contains the index of the current data record
within the current individual record, and LIREC contains the index of
the last data record within the current individual record.
(See Record Counters: NIREC,NDREC)
(See Size of Individual Record).
If a data record does not contain an observation (i.e. MDV=1), this
code sets RS to 0.
If a variable is to be used only during a copying pass, but recur-
sively across data records (such as N and AV), then it should be
defined in an initialization block (e.g. at ICALL.EQ.0). Setting
PRDFL=1 during a noncopying pass, assures that during such passes,
PRED output will be accepted in the usual way: data record by data
record. Setting PRDFL=1 during a copying pass, but only with the sec-
ond pass through the data records, assures that PRED output (which
includes the RS variable) will be accepted data record by data record,
but only during the repeated pass.REFERENCES: none
59 OBJECTIVE FUNCTION VALUE INDIVIDUAL EXAMPLE
Note: With NONMEM 7, the additional output file root.phi contains the same information. (See additional_output_file). This example illustrates how individual contributions to the objective function may be obtained. (See Objective Function Value Individual). This example is a modification of CONTROL4 from the NONMEM distribution media. It is for NONMEM VI 2.0 and later.
$PROB THEOPHYLLINE POPULATION DATA
$INPUT ID DOSE TIME CP=DV WT
$DATA THEO
$PRED
;THETA(1)=MEAN ABSORPTION RATE CONSTANT (1/HR)
;THETA(2)=MEAN ELIMINATION RATE CONSTANT (1/HR)
;THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
;DOSE=WT-ADJUSTED DOSE (MG/KG)
;DS=NON-WT-ADJUSTED DOSE (MG)
IF (DOSE.NE.0) THEN
DS=DOSE*WT
W=WT
ENDIF
KA=THETA(1)+ETA(1)
KE=THETA(2)+ETA(2)
CL=THETA(3)*W+ETA(3)
D=EXP(-KE*TIME)-EXP(-KA*TIME)
E=CL*(KA-KE)
F=DS*KE*KA/E*D
Y=F+EPS(1)
IF (ICALL.EQ.3) THEN
WRITE (70,*) IIDX,CNTID
ENDIF
$THETA (.1,3,5) (.008,.08,.5) (.004,.04,.9)
$OMEGA BLOCK(3) 6 .005 .0002 .3 .006 .4
$SIGMA .4
$EST MAXEVAL=450 PRINT=5
Pairs of values are written, one pair per line, one pair for each
individual record. Each pair contains the ID data item and the con-
tribution the objective function.
1.000000000000000E+00 1.328986294492194E+01
2.000000000000000E+00 2.071754196741930E+01
...
They may also be displayed in a table, using
$ABBR COMRES=2
and code such as the following in the $ERROR or $PK block:
IF (COMACT.EQ.1) THEN
COM(1)=IIDX(NIREC)
COM(2)=CNTID(NIREC)
ENDIF
The following, for example, will produced a separate table for the
values:
$TABLE IID=COM(1) CNT=COM(2) FILE=comvals NOAPPEND NOPRINT FIRSTONLY
Note: With earlier versions than NONMEM 7.3, verbatim code is needed:
" COM(1)=IIDX(NIREC)
" COM(2)=CNTID(NIREC)REFERENCES: None
60 STOCHASTIC DIFFERENTIAL EQUATION EXAMPLES
A number of examples are given of the use of Stochastic Differential Equations. sde6 Example problem without sde sde7 Example problem with SDE, using expanded data set that have repeated lines to re-invoke PRED several times for each data record. sde8 Use of XVID to re-invoke PRED several times per data record wthout requiring the redundancy of sde7.csv. sde9 Alternative method of sde programming, using add-on sde.f90 They, along with the data files, can be found in the examples directory of the NONMEM 7 distribution medium. See INTRODUCTION TO NONMEM 7, Repeated Observation Records See INTRODUCTION TO NONMEM 7, Stochastic Differential Equation Plug-In
REFERENCES: Guide Introduction_7
61 SIMULATION_1 EXAMPLE
This is an example of simulation with population data. The presence of the $SIMULATION record causes NONMEM to generate DV data items. The data set in file SIMORIG may contain nulls in the DV column. Table file SIMDATA1 is created. It is similar to SIMORIG, except that the DV column contains the simulated observations. Because of the option NOHEADER of $TABLE, table header lines are not present, and, with the option NOPRINT, the table is not included in the NONMEM output report. SIMDATA1 can be used in subsequent runs as an input file for NONMEM analysis without further modification. These runs can involve different initial estimates, different models, etc.
$PROBLEM Simulation of population data
$INPUT ID TIME WT AMT DV
$DATA SIMORIG IGNORE=#
$SUBROUTINE ADVAN1
$PK
CL=THETA(1)*EXP(ETA(1))
V=THETA(2)*EXP(ETA(2))
K=CL/V
S1=V
$ERROR
Y=F+F*EPS(1)
$THETA .0625 10
$OMEGA .09 .05
$SIGMA .01
$SIMULATION (9215690) ; seed 1-7 digits
$TABLE ID TIME WT AMT NOPRINT FILE=SIMDATA1 NOHEADER
Contents of SIMORIG (for the first two individuals) follow.
# ID TIME WT AMT DV
1 0 80 100 .
1 1 80 . .
2 0 75 100 .
2 2 75 . .
Contents of SIMDATA1 follow, with floating-point notation converted to
fixed-point notation for the convenience of the reader. With versions
of NONMEM prior to 7.4, this help item contained an incorrect set of
values. The correct values are as follows.
1. 0. 80. 100. 0.000 10.000 0.000 0.000
1. 1. 80. 0. 13.646 9.938 3.708 1.531
2. 0. 75. 100. 0.000 10.000 0.000 0.000
2. 2. 75. 0. 8.379 9.876 -1.496 -0.625
The last four columns are automatically included by NONMEM in every
table (unless the NOAPPEND option is used). They are: DV (simulated
values), PRED (calculated using zeros for etas and epsilons), RES, and
WRES. The DV, RES, and WRES colunns are zero for dose event records.
The last three columns may be ignored (or the DROP option of $INPUT
may be used) when the file is subsequently used. Alternately, list DV
on the $TABLE record, and include the option NOAPPEND.
It is possible to display simulated values of PRED-defined variables.
With NM-TRAN, simulated etas are automatically displayable.
Add the following to the $ERROR block to capture simulated epsilons as
ERROR-defined items:
EP1=EPS(1)
One may change the $SIMULATION and $TABLE records as follows.
$SIMULATION (9215690) ONLYSIM
$TABLE ID TIME WT AMT CL V ETA1 ETA2 EP1
NOPRINT FILE=SIMDATA1 NOHEADER
Displayed values of CL, V, ETA1, ETA2, and EP1 are simulated values.
Without the ONLYSIM option of $SIMULATION, the displayed values are
typical values.
(See Simulation example 2, Simulation example 3).REFERENCES: Guide IV, section III.B.13 , IV.I; Guide V, section 12.4.8; Guide VI, section III.C , III.E , IV.B , IV.G.1
62 SIMULATION_2 EXAMPLE
This is an example of simulation with population data. It is a continuation of Simulation Example 1. First, during the Simulation Step, values of the covariate WT are simulated, as well as DV values, as PK-defined items SIMWT. Normallydistributed values of ETA(3) are used. Values of CL are calculated based on SIMWT. SIMWT replaces WT in the table file, and, because the ONLYSIM option is used, tabled values of SIMWT are simulated values. Second, this process is repeated 10 times, using the same model, by means of the SUBPROBLEMS option. The table file contains the concatenated output of all the subproblems, and is the simulated data for 20 individuals.
$PROBLEM Simulation of population data and weight
$INPUT ID TIME WT AMT DV
$DATA SIMORIG IGNORE=#
$SUBROUTINE ADVAN1
$PK
SIMWT=70+70*ETA(3)
CL=THETA(1)*SIMWT*EXP(ETA(1))
V=THETA(2)*EXP(ETA(2))
K=CL/V
S1=V
$ERROR
Y=F+F*EPS(1)
$THETA .0625 10
$OMEGA .09 .05 .04
$SIGMA .01
$SIMULATION (9215690) ONLYSIM SUBPROBLEMS=10
$TABLE ID TIME SIMWT AMT NOPRINT FILE=SIMDATA2 NOHEADER
SIMORIG is the same as with Simulation Example 1. Contents of SIMORIG
(for the first two individuals) follow.
# ID TIME WT AMT DV
1 0 80 100 .
1 1 80 . .
2 0 75 100 .
2 2 75 . .
Contents of SIMDATA2 follow, with floating-point notation converted to
fixed-point for the convenience of the reader.
1. 0. 77.486 100. 0.000 10.000 0.000 0.000
1. 1. 77.486 0. 6.825 6.457 0.369 0.000
2. 0. 72.882 100. 0.000 10.000 0.000 0.000
2. 2. 72.882 0. 3.515 4.169 -0.653 0.000REFERENCES: Guide IV, section III.B.13 , IV.I; Guide V, section 12.4.8; Guide VI, section III.C , III.E , IV.B , IV.G.1
63 SIMULATION_3 EXAMPLE
This is an example of simulation and estimation within the same NONMEM run. It is a continuation of Simulation Example 2. As in Simulation Example 2, values of the covariate WT are simulated during the Simulation Step. They also replace the original WT values in NONMEM’s internal copy of the data set. This transgeneration makes them available in the Estimation and Table Steps. There are two models for V, illustrating the use of diferent models in the Simulation and Estimation Steps: (simulated V) VS = THETA(2) * EXP(ETA(3)) (non-simulated V) VN = THETA(3) + ETA(4) There are two epsilons in the $ERROR block, illustrating the use of an of using an initial estimate for estimation different from the value used to simulate. Because theta(2), eta(2), eta(3) and eps(1) are used only in the Simulation Step, the FIXED attribute is specified for theta(2), omega(2,2), omega(3,3) and sigma(1,1) during the Estimation Step.
$PROBLEM Simulation of population data w/ covariates
$INPUT ID TIME WT AMT DV
$DATA SIMORIG IGNORE=#
$SUBROUTINE ADVAN1
$PK
IF (ICALL.EQ.4) THEN
WT=70+70*ETA(3)
ENDIF
CL=THETA(1)*WT*EXP(ETA(1))
IF (ICALL.EQ.4) THEN
V=THETA(2)*EXP(ETA(2)) ; Simulated V
ELSE
V=THETA(3)+ETA(4) ; Non-simulated V
ENDIF
K=CL/V
S1=V
$ERROR
IF (ICALL.EQ.4) THEN
Y=F+F*EPS(1)
ELSE
Y=F+F*EPS(2)
ENDIF
$THETA (0,.0625) (10 FIX) (0,12)
$OMEGA .09 (.05 FIX) (.04 FIX) 4
$SIGMA (.01 FIX) .03
$SIMULATION (9215690) SUBPROBLEMS=10
$ESTIMATION
$TABLE ID TIME WT AMT NOPRINT FILE=SIMDATA3 NOHEADERREFERENCES: Guide IV, section III.B.13 , IV.I; Guide V, section 12.4.8; Guide VI, section III.C , III.E , IV.B , IV.G.1
64 STIELTJES EXAMPLE
The following are examples of the specification of the Stieltjes method.
$EST PRINT=1 MET=1 STIELTJES
$EST MAX=9999 METH=1 LAPLACE LIKE PRINT=1 STIELTJES
$EST MAX=9999 METH=1 LAPLACE LIKE PRINT=1 STIELTJES GRID=(9,10,.1,.9)
$EST MAX=500 MET=COND LAPLACE LIKE STIELTJES GRID=(20,2,.1,.9)REFERENCES: Guide Introduction_7
65 SUPERPROBLEM_1 EXAMPLE
There are four problems: Problem 1: The (true) probability that with an individual chosen at random, after a 200 mg dose, the true concentration at 5 hrs is above 5, is computed using simulation. Problem 2: A population data set is simulated. Problem 3: Parameter estimates are obtained. Problem 4: The probability in question is estimated, using simulation and the parameter estimates obtained in problem 3. Problems 2-4 comprise a superproblem. The bias in the estimator defined by problem 4 is estimated by iterating this superproblem.
$PROB COMPUTE PERCENTAGE OF POP. WITH CP>5 AT DOSE=200, TIME=5
$INPUT ID TIME DOSE DV
$DATA dataB
$PRED
IF (ICALL.EQ.0) THEN
BIAS=0 ; initialize
TRUE=0
OPEN (50,FILE='BIAS') ; Give the output file a name
ENDIF
IF (ICALL.EQ.1) THEN
IF (IPROB.EQ.1.OR.IPROB.EQ.4) N=0 ; see rocm14
CALL SUPP (1,1)
ENDIF
KA=THETA(1)*EXP(ETA(1))
KE=THETA(2)*EXP(ETA(2))
V=THETA(3)*EXP(ETA(3))
A=EXP(-KE*TIME)
B=EXP(-KA*TIME)
C=KA-KE
D=A-B
E=KA*DOSE/(V*C)
F=E*D
Y=F+ERR(1)
IF (ICALL.EQ.4) THEN
IF (IPROB.EQ.1.OR.IPROB.EQ.4) THEN
IF (F.GT.5) N=N+1
IF (IREP.EQ.NREP) THEN ; see rocm10
PER=100.*N/NREP
IF (IPROB.EQ.1) TRUE=PER
IF (IPROB.EQ.4) THEN
BIAS=BIAS+(PER-TRUE)/TRUE
IF (S1IT.EQ.S1NIT) THEN
BIAS=BIAS/S1NIT
WRITE (50,*) BIAS
ENDIF
ENDIF
ENDIF
ENDIF
ENDIF
$THTA (.4,1.7,7) (.025,.102,.4) (10,29,80)
$OMEGA .04 .04 .04
$SIGMA 1.5
$SIM (5566898) ONLY SUB=1000
;
$SUPER SCOPE=3 ITERATIONS=10
$PROB SIMULATION
$INPUT ID TIME DOSE DV
$DATA dataA
$THTA (.4,1.7,7) (.025,.102,.4) (10,29,80)
$OMEGA .04 .04 .04
$SIGMA 1.5
$SIM (-1) ONLY
$TABLE ID TIME DOSE DV FILE=simulation NOHEADER NOPRINT NOFORWARD
;
$PROB ESTIMATION
$INPUT ID TIME DOSE DV
$DATA simulation (4E12.0) NRECS=500 NOOPEN
$THTA (.4,1.7,7) (.025,.102,.4) (10,29,80)
$OMEGA .04 .04 .04
$SIGMA 1.5
$EST PRINT=0 MSFO=msf
;
$PROB ESTIMATE PERCENTAGE OF POP. WITH CP>5 AT DOSE=200, TIME=5
$INPUT ID TIME DOSE DV
$DATA dataB
$MSFI msf
$SIM (-1) ONLY SUB=100 TRUE=FINAL
Comments:
(1) Variables whose values are to be retained across problems
(e.g. BIAS and TRUE) should be defined in initialization blocks.
In this instance these variables are defined at ICALL=0, before
any of the problems are implemented
(2) IREP and NREP are the number of the current replication and
number of requested replications for the Simulation Step, respec-
tively.
(3) S1IT and S1NIT are the number of the current iteration and
number of requested iterations for the (level 1) superproblem,
respectively.
(4) Problems 2 & 3 can, of course, be combined, making the con-
trol stream much simpler. But separation can sometimes be use-
ful, and here we illustrate how this can be done.
(5) The data set of problem 3 can be taken to be the internal
data set created with problem 2 simply by removing both the $TA-
BLE record of problem 2 and the $DATA record of problem 3. But
these records are included here to illustrate the ability to use
a table file as a data file in a subsequent problem of the same
run. When so doing, all the items in the $INPUT record of the
subsequent problem should be included in the $INPUT record of the
problem with the $TABLE record. (If a PRED-item in the table is
to be used in the subsequent problem as a data item, the PRED-
item name can be different from the data item name.) The FORMAT
and NOOPEN options should be included in the $DATA record of the
subsequent problem. The NRECS option is not needed, but if it is
also included, it should be set to the number of records in the
table file, i.e. to the number of data records in data set A
(assumed here to be 500).
Contents of dataA (for individual 1 only) follow:
1 .27 320. .
1 .52 320. .
1 1. 320. .
1 1.92 320. .
1 3.5 320. .
1 5.02 320. .
1 7.03 320. .
1 9. 320. .
1 12. 320. .
1 24.3 320. .
This data set cannot be run as-is. File dataA must include more sub-
jects. If not, the estimation step in problem 3 will not terminate
correctly. The final problem will produce the message:
THE SIMULATION TASK REQUIRES THAT A BETTER FINAL ESTIMATE
BE AVAILABLE IN THE MODEL SPECIFICATION FILE
OR
MULTIPLE SUB_PROBLEMS IS INCOMPATIBLE WITH OTHER INFORMATION
IN THE CONTROL STREAM AND INPUT MSF
This error message is related to an earlier message
PARAMETER ESTIMATE IS NEAR ITS BOUNDARY
THIS MUST BE ADDRESSED BEFORE THE COVARIANCE STEP CAN BE IMPLEMENTED
A larger version of dataA can be found in the examples directory of
the NONMEM distribution medium. The content of dataA above is repli-
cated 28 times. ID values 1 and 2 alternate so that there are 28 sub-
jects.
Contents of dataB follow:
1 5.0 200. .
Comments on NONMEM 7 changes:
In this example as distributed with earlier versions of NONMEM,
the WRITE statement specified unit 42:
WRITE (42,*) BIAS
With NONMEM 7, this must be changed to
WRITE (50,*) BIAS
This avoids a possible NONMEM error message
FORTRAN UNIT BEING OPENED OR CLOSED IN PROBLEM 1 IS OUTSIDE
ACCEPTABLE RANGE (50-2000)
A new statement is added:
OPEN (50,FILE='BIAS')
This causes the output to be written to a file named BIAS rather
than fort.50 or the like.REFERENCES: None
66 SUPERPROBLEM_2 EXAMPLE
Two problems, in order to effect a standard two-stage analysis: Problem 1: First stage analysis Problem 2: Obtain simple population average of PK parameter. Problem 2 comprises a superproblem. It is implemented once for CL and once for V, and can be implemented more times when there are more parameters of interest.
$PROBLEM TWO-STAGE STAGE 1
$DATA dataD
$INPUT ID DOSE TIME DV=CP
$PRED
IF (IPROB.EQ.1.AND.NEWIND.LE.1) D=DOSE
IF (IPROB.EQ.1) THEN
CL=THETA(1)*EXP(ETA(1))
V =THETA(2)*EXP(ETA(2))
K=CL/V
F=D*EXP(-K*TIME)/V
Y=F+F*EXP(ETA(3))*ERR(1)
ENDIF
IF (ICALL.EQ.4) THEN
; if S1IT >= 1 then problem 2 is being implemented.
; see comment 1 below.
IF (S1IT.GE.1) THEN
" IF (S1IT.EQ.1) DPV=DATREC(2)
" IF (S1IT.EQ.2) DPV=DATREC(3)
Y=DPV
RETURN
ENDIF
ENDIF
IF (S1IT.EQ.1) Y=THETA(1)+ERR(1)
IF (S1IT.EQ.2) Y=THETA(2)+ERR(2)
$THETA 5 50
$OMEGA 100 100 100
$SIGMA 1 FIX 0 FIX
$ESTIMATION MAXEVAL=0 POSTHOC METH=COND INTERACTION
$TABLE ID CL V FILE=stage2 NOHEADER NOPRINT FIRST NOFORWARD
;
$SUPER SCOPE=1 ITERATION=2
$PROBLEM TWO-STAGE STAGE 2
$INPT ID CLI VI DV
$DATA stage2 (4F12.0) NRECS=12 REWIND NOOPEN
$THETA 5 50
$OMEGA 0 FIX 0 FIX 0 FIX
$SIGMA 1 1
$SIMULATION (1)
$ESTIMATION ;population estimates are obtained
Comment:
(1) The simulation step in problem 2 is used to select the DV
data item from one of the several possibilities (CL or V) output
by problem 1.
(2) Using generated code, data item names that do not appear in
the $INPUT record of problem 1, such as CLI and VI, cannot be
used in abbreviated or verbatim code. DATREC is a vector-array
of data items that can be referenced in verbatim code.
Contents of dataD (for individual 1 only) follow.
1.000 4.020 0.000 0.091
1.000 0.000 0.250 0.084
1.000 0.000 0.570 0.105
1.000 0.000 1.120 0.057
1.000 0.000 2.020 0.038
1.000 0.000 3.820 0.060
1.000 0.000 5.100 0.046
1.000 0.000 9.050 0.017
1.000 0.000 7.030 0.024
1.000 0.000 12.120 0.018
1.000 0.000 24.370 0.002REFERENCES: None
67 TIME AFTER DOSE (TAD) EXAMPLE
This fully-worked out example shows how Time After Dose (TAD) may be computed in $PK abbreviated code. It works with transient and steadystate doses, and also with additional and lagged doses. It is based on suggestions from the NONMEM tips directory but is not identical to any of them. See Note 8 for a discussion of the tips. Two versions of the code are given. VERSION I Code to compute TADE (“TAD Effective”). The time that the dose enters the system is used, which is the time the dose was administered plus ALAG. VERSION II Code to compute TADA (“TAD Administered”). The time that the dose was administered is used, and ALAG is ignored. If PDxPop is used or any other software that looks for TAD in the table, specify TAD=TADA or TAD=TADE in the $TABLE record.
$PROB TIME AFTER DOSE (Based on RUN# 705 nonmem coding challenge #1)
; This example shows how Time After Dose can be computed.
; It contains two independent versions of the computation.
; Version I computes TADE (Time After Dose Effective)
; If there is absorption lag, TADE is the Time after the lagged dose,
; that is, the time the dose actually enters the sysem.
; Version II computes TADA (Time After Dose Administered)
; Absorption lag is ignored.
;
; As of Dec. 17, 2018
;
$INPUT C ID TIME DV AMT WT AGE CRCL SMK ADDL II EVID DROP
$DATA 706.csv IGNORE=C
$SUBROUTINE ADVAN2 TRANS2
$PK
;
ALAG1=0 ; ALAG1 is 0 for this example, but could be set to some other
; value. With ALAG1=0, either version can be used.
; However, VERSION I IS simpler and is preferred.
;
; VERSION I: CODE TO COMPUTE TADE = TAD EFFECTIVE (USE ALAG)
; PREDPP KEEPS TRACK INITIATING AND IMPLIED DOSES, AND COMPUTES DOSTIM
IF (NEWIND.LT.2.OR.EVID.EQ.3) THEN
TDOSE=-999 ; TIME OF MOST RECENT DOSE. -999 IF NO PREVIOUS DOSE.
TADE=0.0 ; TIME AFTER DOSE EFFECTIVE
ENDIF
IF (EVID.EQ.1.AND.ALAG1.EQ.0.OR.EVID.EQ.4.AND.ALAG1.EQ.0) TDOSE=TIME
IF (DOSTIM.GT.0) TDOSE=DOSTIM
IF (TDOSE.GT.-999) TADE=TIME-TDOSE
; If no ALAG and no ADDL doses, the above three lines become:
; IF (EVID.EQ.1.OR.EVID.EQ.4) TDOSE=TIME
; IF (TDOSE.GT.-999) TADE=TIME-TDOSE
;
; End of VERSION I
; VERSION II: CODE TO COMPUTE TADA = TAD ADMINISTERED (IGNORES ALAG)
; IF ALAG>0, PREDPP's DOSTIM CANNOT BE USED FOR TADA
; THE ABBREVIATED CODE MUST KEEP TRACK OF INITIATING AND IMPLIED DOSES.
; THIS CODE IS SIMILAR TO tip4-new-general-01-31-18.txt
IF (NEWIND.LT.2.OR.EVID.EQ.3) THEN
TDOSA=-999 ; TIME OF MOST RECENT DOSE. -999 IF NO PREVIOUS DOSE.
TADA=0.0 ; TIME AFTER DOSE ADMINISTERED
ENDIF
IF (DOSTIM==0) THEN ; IGNORE NON-EVENT DOSE TIMES (DOSTIM>0)
IF(EVID.EQ.1.OR.EVID.EQ.4) THEN ; NEW DOSE EVENT RECORD
TADA=0.0
DIV=II
TDOSA=TIME
TLAST=TDOSA+ADDL*II ; TLAST IS THE TIME OF THE FINAL IMPLIED DOSE
ENDIF
IF (TDOSA.GE.0.AND.EVID.NE.1.AND.EVID.NE.4) THEN ; THERE WAS AN EARLIER DOSE
IF (TIME>TLAST) THEN
TADA=TIME-TLAST ; CURRENT TIME IS PAST THE TIME OF THE FINAL IMPLIED DOSE
ELSE
DIFF=TIME-TDOSA
TADA=MOD(DIFF,DIV) ; COMPUTES TIME OF THE MOST RECENT IMPLIED DOSE
ENDIF
ENDIF ; NOT A DOSE
ENDIF ; END OF DOSTIM==0
;
; End of VERSION II
CL=THETA(1)*EXP(ETA(1))
TVV=THETA(2)
V=TVV*EXP(ETA(2))
TVKA=THETA(3)
KA=TVKA*EXP(ETA(3))
S2=V
$THETA
(0, 10)
(0, 100)
(0, 0.5)
$ERROR
Y=F+ERR(1)
IPRED=F
$OMEGA
0.04 ;[P] INTERIND VAR IN CL
0.04 ;[P] INTERIND VAR IN V
0.04 ;[P] INTERIND VAR IN KA
$SIGMA
0.2 ;[A] ADDITIVE COMPONENT
$EST MAXEVAL=0 PRINT=2 NOABORT ;POSTHOC
;
; If PDxPop is used or any other software that looks for TAD in the table, specify
; TAD=TADA or TAD=TADE in the $TABLE record
$TABLE ID TIME EVID TDOSA TADA TDOSE TADE IPRED NOAPPEND NOPRINT ONEHEADER
FILE=tadexa.tab FORMAT=SF11.7
The data for the first subject in 706.csv is:
C,Data Desc: test of TAD for ADDL,,,,,,,,,,
C,ID,TIME,DV,AMT,WT,AGE,CRCL,SMK,ADDL,II,EVID,TAD
0,1,0,0,0,58.4,51,4.49,1,0,0,2,0
0,1,0.1,0,500,58.4,51,4.49,1,1,24,1,0
0,1,1.1,1.67538,0,58.4,51,4.49,1,0,0,0,0
0,1,4.1,2.79283,0,58.4,51,4.49,1,0,0,0,0
0,1,10.1,0.583263,0,58.4,51,4.49,1,0,0,0,0
0,1,24.1,0.146255,0,58.4,51,4.49,1,0,0,0,0
0,1,25.1,1.67538,0,58.4,51,4.49,1,0,0,0,0
The original data set 704.csv is embedded in tip3-2-25-02.txt. Data
set 706.csv can be found in the tips directory. It was created from
704.csv by appending an extra column, TAD, whose values are zero.
$INFN code in the tips assign values to the TAD data item (see Note
8). The code in the present example does not use the TAD data item
(DROP is listed in $INPUT). Instead, TAD is a user-defined variable
listed in MODULE NMPRD4 from which its values may be displayed in a
table file.
Note 1.
It is assumed that TAD is displayed in tables for the purpose of
graphical display of the data. It is not a random variable, and
is not intended to be part of the prediction model.
Note 2.
If there are non-dose records prior to the first dose record of
the Individual record, TAD is set to zero for these records.
This is an arbitrary convention; the user could use different
values, especially if these conditions might signal an error in
the data set.
Note 3.
ALAG1 is used in the example, but a different dosing compartment
could be used, e.g., ALAG2 for doses into compartment 2. The
code would need modification to compute TAD for doses into two or
more compartments.
Note 4.
This code can handle random DOSTIM (Eta on ALAG or ETA on mod-
elled duration or rate.)
Note 5.
The code does not use INFN. The information is developed with
passes through the data set at ICALL=2, in which NONMEM calls
PREDPP. During a pass in which INFN is called, PREDPP is not
called. The abbreviated code would have to account for additional
and lagged dose times.
Note 6.
TDOS=-999 until the first dose. TDOS=0 when the first dose is at
TIME 0. This makes it possible to dispense with a flag variable
whose only function is to say "a dose has occurred". If there
are negative time values in the data set and any time value is
smaller than -999, then
Note 7.
The code is intended to compute TAD for multiple individual doses
or one initiating dose. (An initiating dose has ADDL>0,II>0).
Note 8.
Changes were made to tip3, tip4, and tip4update in the NONMEM
tips directory. See tips in
https://nonmem.iconplc.com/nonmem/tips
In each case, the name of the revised tip contains 01-31-18 and
supercedes the original. Although the tips were revised as part
of the nonmem 7.4.2 release they apply to previous versions of
nonmem.
tip3-revised-01-31-18.txt
The original version, tip3-2-25-02.txt, provides a data file
704.csv and a simple model. The revised version provides
improved code, similar to the code in tadexa.exa.
tip4-revised-01-31-18.txt
Uses Fortran suboutine INFN and equivalent $INFN code. The
datafile is 706.csv.
tip4update-revised-01-31-18.txt
Uses DOSTIM.
The original versions are still in the tips directory with
new names that include the original dates.
tip3-2-25-02.txt
tip4-3-04-02.txt
tip4update-10-11-13.txtREFERENCES: None
68 TDM EXAMPLE
This is an example of how F1 can be used to simplify the preparation of a data file when the dose is changed often, e.g., for Therapeutic drug monitoring (TDM). Suppose the dosing pattern is consistent (e.g., every 12 hours) but the amount changes often. Suppose also there are observation events at various times. In the following data set, there is only one dose record having AMT=1. ADDL and II are used to specify an (effectively infinite) number of additional doses, every 12 hours. When the dose amount changes, it is sufficient to supply a record giving the new dose amount in the DOSE data item. Such records have EVID=2. There is no need to interrupt the dosing pattern for an observation event, or to count the number of doses between observation events. Suppose a fragment of the dose records looks like this: ID DAT1=DROP TIME EVID DV DOSE AMT II ADDL 1 01.01.2008 20 1 0 4 1 12 99999 1 03.01.2008 8 0 6.1 0 0 0 0 1 04.01.2008 8 0 7.6 0 0 0 0 1 04.01.2008 20 2 0 3 0 0 0 1 06.01.2008 8 0 5.4 0 0 0 0 Here is a fragment of the control stream:
$INPUT ID DAT1=DROP TIME EVID DV DOSE AMT II ADDL
...
$PK (ADDITIONAL OR LAGGED)
...
IF (DOSE> 0) SAVEDOSE=DOSE
F1=SAVEDOSE
In $PK, the value of DOSE is saved, and used on subsequent calls to PK
as the value of F1. Although F1 is usually used to model the bio-
availability of doses to compartment 1, in this example it is the
actual dose amount.REFERENCES: Guide V, section 7.4.4
69 THREE_COMPARTMENT EXAMPLE
This example shows how to implement a three Compartment Linear Model with First-order Absorption using PREDPP, where elimination occurs from a central compartment and one of two peripheral compartments. As with any ADVAN, one or more doses could be input to any of the compartments. Steady state doses are also possible. The use of the DEFDOSE attribute on the following $MODEL record specifies that if the CMT data item is null on a given dose event record, then the dose(s) specified by that record are to be placed in the depot compartment. The use of the INITIALOFF attribute specifies that the depot compartment is to be off until a dose is placed in it. The use of the DEFOBS attribute on the $MODEL record specifies that if the CMT data item is null on a given observation event record, then the quantity F which is available as a right-hand quantity in $ERROR is to be the scaled drug amount in the central compartment. (See CMT PCMT data item).
$PROBLEM EXAMPLE OF A THREE COMPARTMENT MODEL WITH ABSORPTION
$INPUT ID TIME AMT DV CMT
$DATA datafile
$SUBROUTINES ADVAN7
$MODEL COMP = (DEPOT,DEFDOSE,INITIALOFF)
COMP = (CENTRAL,DEFOBS)
COMP = PERIPH1, COMP = PERIPH2
$PK
K12 = THETA(1)*EXP(ETA(1)) ; depot to central
K23 = THETA(2) ; central to periph1
K32 = THETA(3) ; periph1 to central
K30 = THETA(4)*EXP(ETA(2))
K24 = THETA(5) ; central to periph2
K42 = THETA(6) ; periph2 to central
K20 = THETA(7)*EXP(ETA(3)) ; elimination from central
S2 = THETA(8)*EXP(ETA(4)) ; scale for central
;There are other reasonable ways that etas can be assigned.
;Elaborate these expressions for the K's as is appropriate.
;Add $ERROR, parameter and task records as is necessary/desired.
Other parameterizations are possible, as long as the above set of
micro-constants is defined. E.g., the definitions of K20 and S2 can
be replaced by the following.
CL = THETA(7)*EXP(ETA(3))
V = THETA(8)*EXP(ETA(4))
K20 = CL/V
S2 = V
If the depot compartment is not needed, it can be omitted. The num-
berings of the compartments, the thetas and etas must then be
adjusted.REFERENCES: Guide VI, section VII.C.5 , VII.C.7
70 TNPRI EXAMPLES
These examples illustrate the use of the NONMEM utility routine TNPRI by the NM-TRAN $PRIOR record, and also by an equivalent user-supplied routine PRIOR. For a discussion of PRIOR: (See prior). For a discussion of the TNPRI: (See tnpri). Example 1 Use of PRIOR and TNPRI and a control stream for population data: This example obtains parameter estimates from population data, using a two-compartment PK model, and it incorporates a frequentist prior of transformed normal form for all of THETA, OMEGA and SIGMA. The prior information is found in model specification file msf1.
$PROB READ THE MODEL SPECIFICATION FILE
$DATA data
$INPUT ID DOSE TIME DV WT
$PRIOR TNPRI (PROBLEM 2)
$MSFI msf1 ONLYREAD
$PRED
;
; THETA(1)=MEAN ABSORPTION RATE CONSTANT - MEAN ELIM. RATE CONSTANT (l/HR)
; THETA(2)=MEAN ELIM. RATE CONSTANT (1/HR)
; THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
; DOSE=WEIGHT-ADJUSTED DOSE (MG/KG)
;
IF (NEWIND.NE.2) THEN
AMT=DOSE*WT
W=WT
ENDIF
T0=THETA(1)*EXP(ETA(1))
T2=THETA(2)*EXP(ETA(2))
T1=T2+T0
T3=THETA(3)*W*EXP(ETA(3))
A=EXP(-T2*TIME)
B=EXP(-T1*TIME)
C=T1-T2
D=A-B
E=T3*C
Y=AMT*T1*T2/E*D+EPS(1)
$PROB POPULATION DATA WITH PRIOR ON THETA, OMEGA AND SIGMA
$INPUT ID DOSE TIME DV WT
$DATA data REWIND
$THETA (0,4,5) (0,.09,.5) (.004,.01,.9)
$OMEGA BLOCK (3) .7 .04 .05 .02 .06 .08
$SIGMA .4
$EST
Instead of a $PRIOR record, the following may be used:
$SUBROUTINE PRIOR=prior
The prior routine is as follows:
SUBROUTINE PRIOR (ICALL,CNT)
USE SIZES, ONLY: DPSIZE,ISIZE
USE NMPRD_INT, ONLY: IPROB
REAL(KIND=DPSIZE) :: CNT
REAL(KIND=DPSIZE) :: PLEV
INTEGER(KIND=ISIZE) :: ICALL
IF (IPROB==2) THEN
PLEV=0.
ITYP=0
NSAM=0
ISS=0
IFND=0
MODE=0
IVAR=0
CALL TNPRI(IFND,MODE,ITYP,PLEV,NSAM,ISS,IVAR,CNT)
ENDIF
RETURN
END
Example 2
Use of PRIOR and TNPRI and a control stream for population data:
This example obtains parameter estimates from population data, using a
two-compartment PK model, and it incorporates a frequentist prior of
transformed normal form for all of THETA and OMEGA, but for only the
first element of a diagonal SIGMA, which is a prior-specific parame-
ter. Prior information for these parameters is found in model speci-
fication file msf1.
$PROB READ THE MODEL SPECIFICATION FILE
$DATA data
$INPUT ID DOSE TIME DV WT
$PRIOR TNPRI (PROBLEM 2)
$MSFI msf1 ONLYREAD
$PRED
;
; THETA(1)=MEAN ABSORPTION RATE CONSTANT - MEAN ELIM. RATE CONSTANT (l/HR)
; THETA(2)=MEAN ELIM. RATE CONSTANT (1/HR)
; THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
; DOSE=WEIGHT-ADJUSTED DOSE (MG/KG)
;
IF (NEWIND.NE.2) THEN
AMT=DOSE*WT
W=WT
ENDIF
T0=THETA(1)*EXP(ETA(1))
T2=THETA(2)*EXP(ETA(2))
T1=T2+T0
T3=THETA(3)*W*EXP(ETA(3))
A=EXP(-T2*TIME)
B=EXP(-T1*TIME)
C=T1-T2
D=A-B
E=T3*C
Y=AMT*T1*T2/E*D+EPS(2)
$PROB POPULATION DATA WITH PRIOR ON THETA AND OMEGA
$INPUT ID DOSE TIME DV WT
$DATA data REWIND
$THETA (0,4,5) (0,.09,.5) (.004,.01,.9)
$OMEGA BLOCK (3) .7 .04 .05 .02 .06 .08
$SIGMA .4 FIX .5
$EST
Instead of a $PRIOR record, use the same $SUBROUTINE record and the
same subroutine prior as in Example 1.
Example 3
Use of PRIOR and TNPRI and a control stream for population data:
This example simulates the THETA and OMEGA parameters of the same two-
compartment population PK model, as well as the population data from
this model. As with Example 2, the first element of SIGMA is prior-
specific, and so its value need not be simulated. The value .5 is
used for the value of the second element of SIGMA. The default vari-
ance-covariance matrix from the prior problem is used for the hyperpa-
rameter variance-covariance matrix (IVAR=1).
$PROB READ THE MODEL SPECIFICATION FILE
$DATA data
$INPUT ID DOSE TIME DV WT
$PRIOR TNPRI (PROBLEM 2) PLEV=.9999 IVAR=1
$MSFI msf1 ONLYREAD
$PRED
;
; THETA(1)=MEAN ABSORPTION RATE CONSTANT - MEAN ELIM. RATE CONSTANT (l/HR)
; THETA(2)=MEAN ELIM. RATE CONSTANT (1/HR)
; THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
; DOSE=WEIGHT-ADJUSTED DOSE (MG/KG)
;
IF (NEWIND.NE.2) THEN
AMT=DOSE*WT
W=WT
ENDIF
T0=THETA(1)*EXP(ETA(1))
T2=THETA(2)*EXP(ETA(2))
T1=T2+T0
T3=THETA(3)*W*EXP(ETA(3))
A=EXP(-T2*TIME)
B=EXP(-T1*TIME)
C=T1-T2
D=A-B
E=T3*C
Y=AMT*T1*T2/E*D+EPS(2)
$PROB POPULATION DATA WITH PRIOR ON THETA AND OMEGA
$INPUT ID DOSE TIME DV WT
$DATA data REWIND
$THETA (0,4,5) (0,.09,.5) (.004,.01,.9)
$OMEGA BLOCK (3) .7 .04 .05 .02 .06 .08
$SIGMA .4 FIX .5
$SIM (547676) ONLY TRUE=PRIOR
Instead of a $PRIOR record, the following may be used:
$SUBROUTINE PRIOR=prior
The prior routine is as follows:
SUBROUTINE PRIOR (ICALL,CNT)
USE SIZES, ONLY: DPSIZE,ISIZE
USE NMPRD_INT, ONLY: IPROB
REAL(KIND=DPSIZE) :: CNT
REAL(KIND=DPSIZE) :: PLEV
INTEGER(KIND=ISIZE) :: ICALL
IF (IPROB==2) THEN
PLEV=.9999
IVAR=1
ITYP=0
NSAM=0
ISS=0
IFND=0
MODE=0
CALL TNPRI(IFND,MODE,ITYP,PLEV,NSAM,ISS,IVAR,CNT)
ENDIF
RETURN
ENDREFERENCES: none
71 TRANSIT COMPARTMENT EXAMPLES
A number of examples are provided for the implementation of dose superimposition into a transit compartment. BACKGROUND: These examples are based on the following two papers: Implementation of a transit compartment model for describing drug absorption in pharmacokinetic studies J Pharmacokinet Pharmacodyn (2007) 34:711-726 DOI 10.1007/s10928-007-9066-0 Radojka M. Savic, Daniel M. Jonker, Thomas Kerbusch, Mats. O. Karlsson Implementation of dose superimposition to introduce multiple doses for a mathematical absorption model (transit compartment model) J Pharmacokinet Pharmacodyn (2012) 39:251-262 DOI 10.1007/s10928-012-9247 Jun Shen, Alison Boeckmann & Andrew Vick The first paper defines what a transit compartment is, and gives a model for a single dose into such a compartment. The second paper gives several ways that multiple dosing with dose superimposition can be modelled with NONMEM 7. Files on ftp site and NONMEM 7 The files for the second paper are provided on the IDS ftp site: ftp://nonmem.iconplc.com/Public/nonmem/transit_compartment/ The same files are also provided with NONMEM 7 in directory examples_compartment The general solution (Appendix 4; sumdosetf.ctl) is of particular interest in that it gives an example of the use of vectors (e.g., VECTRA(n), VECTRB(n)) and user-written functions (FUNCA, FUNCB) in abbreviated code. NEW FILES FOR NONMEM 7: The following files may be found in the examples directory. sumdosetn.ctl, sumdosetf.dat With NONMEM 7, it is also possible to code the general solution directly in abbreviated code, using the DO WHILE statement, user-declared arrays, and subscripted variables. Note that data file sumdosetf.csv is identical to data file sall6.csv on the ftp site. As noted in the INTRODUCTION TO NONMEM 7, this includes abbreviated code for summation such as: $abbr declare dosetime(100),dose(100) $abbr declare dowhile i
$DES
INPT=0
I=1
DOWHILE (I<=NDOSE)
IPT=0
IF (T>=dosetime(I)) &
IPT=DOSE(I)*(T-dosetime(I))**NN*EXP(-KTR*(T-dosetime(I)))
INPT=INPT+IPT
I=I+1
ENDDO
ssmultidose.ctl, ssmultidose.dat
This is an example of a steady state reached using multiple tran-
sient bolus doses. The abbreviated code is the same as sum-
dosetn.ctl, but the data file ssmultidose.dat has nine identical
bolus doses, at 12 hour intervals, leading to steady state (to
the number of significant digits printed).
ssonedose.ctl, ssonedose.dat
This is an example of a steady state computed with a single dose
record having SS=1 and II=12. The data file is ssonedose.dat.
The $PK block saves the values of AMT and II from the dose
record. Abbreviated code in the $DES block contains a DO WHILE
loop that adds drug from each implied dose, going back in time
till the amount from the implied dose is negligible.REFERENCES: Guide Introduction_7
72 USER-DEFINED RESERVED FUNCTION EXAMPLE
If you wish to define your own function, and have the information about its proper use of arguments be conveyed upon its execution, so the compiler may detect errors, then one method is to package the definition of the function in a USE module. The function must be listed in an include file whose name starts with the characters nonmem_reserved Caution: NMTRAN will permit the use of such a function but will not compute eta deriatives with respect to the arguments. Do not use such functions to compute any variable that affects the objective function! Here is an example. Myfuncmodule.f90 defines the functions mymin and mymax: MODULE MYFUNCS contains function mymin(a,b,c,d,e) integer mymin integer a,b,c,d,e mymin=min(a,b,c,d,e) end function function mymax(a,b,c,d,e) integer mymax integer a,b,c,d,e mymax=max(a,b,c,d,e) end function END MODULE MYFUNCS Nonmem_reserved_myfunc is the include file that declares its use: ” USE myfuncs, only: mymin,mymax The following control stream file uses the function:
$PROB THEOPHYLLINE POPULATION DATA
$INPUT ID DOSE=AMT TIME CP=DV WT
$DATA THEOPP
$SUBROUTINES ADVAN2 OTHER=myfuncmodule
$PK
;THETA(1)=MEAN ABSORPTION RATE CONSTANT (1/HR)
;THETA(2)=MEAN ELIMINATION RATE CONSTANT (1/HR)
;THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
;SCALING PARAMETER=VOLUME/WT SINCE DOSE IS WEIGHT-ADJUSTED
include "nonmem_reserved_myfunc"
CALLFL=1
KA=THETA(1)+ETA(1)
K=THETA(2)+ETA(2)
CL=THETA(3)*WT+ETA(3)
SC=CL/K/WT
I=mymin(1,2,3,4,5.0)
print *,'I ',I
$THETA (.1,3,5) (.008,.08,.5) (.004,.04,.9)
$OMEGA BLOCK(3) 6 .005 .0002 .3 .006 .4
$ERROR
Y=F+EPS(1)
$SIGMA .4
If you use the wrong argument type (real instead of integer), or per-
haps use the wrong number of arguments, the compiler will readily flag
this.REFERENCES: Guide Introduction_7
73 YLO EXAMPLE
In this example, blood concentrations below a known quantification limit (BQL) have been eliminated from the analysis. A one-sided interval with lower boundary equal to the quantification limit is defined. The conditional likelihood for an observation is “conditioned” on values of the population parameters and also on the condition that the observation is inside the one-sided interval. (See pr_y, YLO YUP) (See YLO YUP Probability: PR_Y). Here is an NM-TRAN control stream.
$PROB YLO example
$INPUT ID TIME DV AMT EVID MDV
$DATA ylow.dat IGNORE=C
$SUBROUTINES ADVAN=1 TRANS=2
$PK
CL=THETA(1)*EXP(ETA(1))
V=THETA(2)*EXP(ETA(2))
K=CL/V
S1=V
$ERROR
; Limit of Quantification (LOQ) is 1
; DV's below LOQ are given MDV=1 in dataset
; PR_Y (PRB) is probability that DV is > LOQ
; PR_Y (PRB) is provided even if MDV is 1
YLO=1
PRB=PR_Y
Y = F+ERR(1)
$THETA
(0,0.693); CL
(0,24); V
$OMEGA
0.09
0.09
$SIGMA
.5
$ESTIMATION
MAXEVAL=9999
NSIG=4
METHOD=COND INTERACTION LAPLACIAN
PRINT=5
$TABLE ID TIME AMT MDV EVID PRB
NOPRINT
ONEHEADER
FILE=ylo.tab
Here is a segment of the dataset for a subject with last DV < LOQ
ID TIME DV AMT EVID MDV
. . . . . .
. . . . . .
. . . . . .
2 0 0 100 1 1
2 0.333 2.47 0 0 0
2 2 3.26 0 0 0
2 6 2.32 0 0 0
2 15.917 2.46 0 0 0
2 24 0 100 1 1
2 48 0 100 1 1
2 72 0 100 1 1
2 95.917 4.04 0 0 0
2 96 0 100 1 1
2 96.08 8.85 0 0 0
2 120 3.16 0 0 0
C Last observation record (below) has DV < LOQ (LOQ = 1) and MDV = 1
C This observation does not contribute to the -2LL estimate. However,
C probability that DV is > LOQ can be tabled using PRB=PR_Y.
2 168 0.163 0 0 1REFERENCES: None