NONMEM model library

Pharmacokinetics

Two-compartment IV population PK

The structural PK model most exposure work rests on. Allometric weight on clearance and volume, between-subject variability on CL and V1, proportional residual error.

The data, before any model

Covariates are summarised once per subject. Summarising per record would count heavily sampled subjects several times over and describe a different population from the one in the study.

Every random effect against every covariate, ranked. Continuous covariates get a rank correlation and categorical ones the gap in median eta between groups, in standard deviations, so the two are on the same scale. The ranking only suggests what to try next. Whether a covariate belongs in the model is decided by fitting it.

CYP genotype against the random effect on CL is by far the strongest signal, at 1.26 standard deviations between groups, and it isn't in the model yet. That's the next run to make. The other 1 above the threshold are weaker and probably the same effect showing up again: the empirical Bayes estimates for the two random effects are correlated, so a real shift in one appears faintly in the other.

Estimates against the truth

Fitted by adaptive Gauss-Hermite quadrature, 5 nodes, objective function 1514.6, 10.3 s. 7 of 7 parameters have a 95% confidence interval containing the value the data was simulated from.

Against the simpler model

Does the distribution phase exist, or would one compartment describe these data just as well?

Refitting as a one-compartment model and re-estimating everything else: dropping it costs 3522.5 objective function on 2 degrees of freedom, p < 0.0001. The null is on the boundary of the parameter space, so the chi-square reference is conservative and the actual p-value is smaller than the one quoted.

Goodness of fit

Shrinkage: CL 0%, V1 9%. Residual (epsilon) shrinkage is 7%. Empirical Bayes estimates are a compromise between a subject's own data and the population, so a subject with little information ends up close to the population value and the spread of the etas understates OMEGA. All low, so the individual predictions below really are individual and the etas are usable for spotting covariate relationships.

The control stream

pk_2cmt_iv.mod passes every static check

$PROBLEM Two-compartment IV infusion, allometric weight, IIV on CL and V1
; Data: data/pk_2cmt_iv.csv, simulated by nmlib.simulate.simulate_pk_2cmt
; The workhorse structural PK model: everything downstream in this library
; that needs an exposure driver assumes something of this shape.
;
; Four covariates are carried on the data. Only WT is in the model, and
; that is the point of the covariate screen on the dashboard: AGE and SEX
; have no effect to find, and CYP -- a metaboliser genotype that really
; does lower clearance -- is deliberately left out so the screen has
; something to catch. Adding it is the obvious next run:
;
;   CL = THETA(1)*(WT/70)**0.75*EXP(ETA(1))*(1 + THETA(5)*CYP)
;
; Until then the OMEGA on CL is carrying the genotype split as if it were
; unexplained between-subject variability, which is exactly what an
; unmodelled covariate does to a population model.

$INPUT ID TIME AMT RATE DV MDV EVID CMT WT AGE SEX CYP
$DATA ../data/pk_2cmt_iv.csv IGNORE=@

$SUBROUTINE ADVAN3 TRANS4

$PK
; Allometry on a 70 kg reference. Exponents are fixed at the conventional
; 0.75 for flows and 1 for volumes; estimate them instead only if the data
; can carry it, which a single study rarely can.
WTCL = (WT/70)**0.75
WTV  = (WT/70)**1

CL = THETA(1)*WTCL*EXP(ETA(1))
V1 = THETA(2)*WTV *EXP(ETA(2))
Q  = THETA(3)*WTCL
V2 = THETA(4)*WTV

S1 = V1

$ERROR
IPRED = F
Y     = IPRED*(1 + EPS(1))

$THETA
(0, 5.0)    ; 1 CL  (L/h) at 70 kg
(0, 15.0)   ; 2 V1  (L)   at 70 kg
(0, 3.0)    ; 3 Q   (L/h) at 70 kg
(0, 25.0)   ; 4 V2  (L)   at 70 kg

$OMEGA
0.09        ; 1 IIV on CL  (var on log scale, ~30% CV)
0.06        ; 2 IIV on V1  (~25% CV)

$SIGMA
0.0225      ; 1 proportional residual error (~15%)

$ESTIMATION METHOD=1 INTERACTION MAXEVAL=9999 SIG=3 PRINT=5 NOABORT
$COVARIANCE UNCONDITIONAL
$TABLE ID TIME DV IPRED PRED CWRES CL V1 Q V2 WT ETA(1) ETA(2)
       ONEHEADER NOPRINT FILE=sdtab_pk

Oncology

Tumour growth inhibition (Claret)

Exponential growth with a drug kill term that decays over time. Without the decay the model can't produce regrowth while treatment continues, which is what's seen in practice and what drives survival predictions.

The data, before any model

Estimates against the truth

Fitted by adaptive Gauss-Hermite quadrature, 3 nodes, objective function 6985.0, 6.1 s. 8 of 8 parameters have a 95% confidence interval containing the value the data was simulated from.

Against the simpler model

Without the resistance term the model can't regrow a tumour while treatment continues. Does the data show regrowth?

Refitting as a plain kill model, with no loss of drug effect and re-estimating everything else: dropping it costs 437.6 objective function on 1 degree of freedom, p < 0.0001. The null is on the boundary of the parameter space, so the chi-square reference is conservative and the actual p-value is smaller than the one quoted.

Goodness of fit

Shrinkage: Y0 2%, KL 7%, KD 45%. Residual (epsilon) shrinkage is 10%. Empirical Bayes estimates are a compromise between a subject's own data and the population, so a subject with little information ends up close to the population value and the spread of the etas understates OMEGA. KD is high enough to matter. Those estimates have been pulled a long way back towards the population, so the individual-prediction panel below makes that parameter look better than it is, and an eta-versus-covariate plot on it would tell you very little.

The control stream

tgi_claret.mod passes every static check

$PROBLEM Claret tumour growth inhibition with resistance
; Data: data/tgi_claret.csv, simulated by nmlib.simulate.simulate_tgi_claret
;
; dY/dt = KL*Y - KD*EXPO*exp(-LAMBDA*t)*Y
;
; KL is exponential growth, KD the drug kill rate, and LAMBDA the resistance
; term that lets the kill effect fade with time. Without LAMBDA the model
; cannot reproduce regrowth while treatment continues, which is the
; behaviour the model exists to describe.

$INPUT ID TIME DV MDV EVID EXPO ARM
$DATA ../data/tgi_claret.csv IGNORE=@

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(TUMOUR, DEFDOSE)

$PK
Y0     = THETA(1)*EXP(ETA(1))
KL     = THETA(2)*EXP(ETA(2))
KD     = THETA(3)*EXP(ETA(3))
LAMBDA = THETA(4)

A_0(1) = Y0

$DES
; EXPO is the subject's average exposure, constant within a subject here.
KILL       = KD*EXPO*EXP(-LAMBDA*T)
DADT(1)    = KL*A(1) - KILL*A(1)

$ERROR
IPRED = A(1)
Y     = IPRED*(1 + EPS(1))

$THETA
(0, 50.0)    ; 1 Y0     baseline tumour size (mm)
(0, 0.006)   ; 2 KL     growth rate (1/day)
(0, 0.0004)  ; 3 KD     kill rate per unit exposure
(0, 0.015)   ; 4 LAMBDA resistance rate (1/day)

$OMEGA
0.18         ; 1 IIV on Y0
0.22         ; 2 IIV on KL
0.30         ; 3 IIV on KD

$SIGMA
0.0144       ; 1 proportional residual error (~12%)

$ESTIMATION METHOD=1 INTERACTION MAXEVAL=9999 SIG=3 PRINT=5 NOABORT
$COVARIANCE UNCONDITIONAL
$TABLE ID TIME DV IPRED PRED CWRES EXPO ARM Y0 KL KD
       ONEHEADER NOPRINT FILE=sdtab_tgi

PK/PD

Indirect response, inhibition of production

Dayneka and Jusko's model I. The drug acts on the turnover of a biomarker rather than on the biomarker itself, so the response lags exposure and washout depends on the biomarker's own loss rate rather than the drug's half-life.

The data, before any model

Estimates against the truth

Fitted by adaptive Gauss-Hermite quadrature, 5 nodes, objective function 11711.2, 216.8 s. 7 of 7 parameters have a 95% confidence interval containing the value the data was simulated from.

Against the simpler model

In a direct model the response tracks concentration exactly. Is there enough lag in the data to justify the indirect structure?

Refitting as a direct effect model and re-estimating everything else: the two have the same number of parameters, and the full model fits 872.4 objective function better. The models aren't nested, so there's no likelihood ratio test, but the objective functions can be compared directly.

Goodness of fit

Shrinkage: KOUT 0%, IC50 44%. Residual (epsilon) shrinkage is 4%. Empirical Bayes estimates are a compromise between a subject's own data and the population, so a subject with little information ends up close to the population value and the spread of the etas understates OMEGA. IC50 is high enough to matter. Those estimates have been pulled a long way back towards the population, so the individual-prediction panel below makes that parameter look better than it is, and an eta-versus-covariate plot on it would tell you very little.

The control stream

pkpd_idr_inhibition.mod passes every static check

$PROBLEM Indirect response, drug inhibiting production of a biomarker
; Data: data/pkpd_idr_inhibition.csv
;   simulated by nmlib.simulate.simulate_idr_inhibition
;
; dR/dt = KIN*(1 - IMAX*C/(IC50 + C)) - KOUT*R
;
; This is Dayneka and Jusko's model I. The response lags exposure because
; the biomarker has its own turnover; washout after stopping is governed by
; KOUT, not by the drug's half-life. Baseline is the untreated steady state
; KIN/KOUT, so it is not a separate parameter.
;
; The design matters as much as the structure here. The drug is given daily
; for a week and then stopped, and the biomarker's half-life (about 3.5
; days) is roughly five times the drug's (about 17 hours). That separation
; is what makes the two clocks distinguishable: the response is still
; falling after the concentration has reached steady state, and still
; recovering weeks after the last dose. With a biomarker that turned over
; as fast as the drug cleared, this model and a direct effect model would
; fit the same data equally well and there would be nothing to estimate.

$INPUT ID TIME AMT DV MDV EVID CMT DOSE
$DATA ../data/pkpd_idr_inhibition.csv IGNORE=@

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(CENTRAL, DEFDOSE)   ; 1 drug
COMP=(RESPONSE)           ; 2 biomarker

$PK
CL   = THETA(1)
V    = THETA(2)
KIN  = THETA(3)
KOUT = THETA(4)*EXP(ETA(1))
IMAX = THETA(5)
IC50 = THETA(6)*EXP(ETA(2))

S1 = V

; The biomarker starts at its untreated steady state.
A_0(2) = KIN/KOUT

$DES
CONC    = A(1)/V
INH     = 1 - IMAX*CONC/(IC50 + CONC)
DADT(1) = -(CL/V)*A(1)
DADT(2) = KIN*INH - KOUT*A(2)

$ERROR
IPRED = A(2)
Y     = IPRED*(1 + EPS(1))

$THETA
(0, 4.0)          ; 1 CL   (L/h), driving PK, fixed from the PK analysis
(0, 100.0)        ; 2 V    (L), likewise
(0, 0.80)         ; 3 KIN  production rate (units/h)
(0, 0.008)        ; 4 KOUT loss rate constant (1/h); baseline KIN/KOUT = 100
(0, 0.80, 1)      ; 5 IMAX maximum fractional inhibition, bounded at 1
(0, 8.0)          ; 6 IC50 concentration at half maximum inhibition (mg/L)

$OMEGA
0.09              ; 1 IIV on KOUT
0.22              ; 2 IIV on IC50

$SIGMA
0.01              ; 1 proportional residual error (~10%)

$ESTIMATION METHOD=1 INTERACTION MAXEVAL=9999 SIG=3 PRINT=5 NOABORT
$COVARIANCE UNCONDITIONAL
$TABLE ID TIME DV IPRED PRED CWRES DOSE KOUT IC50
       ONEHEADER NOPRINT FILE=sdtab_idr

Survival

Time to event, Weibull hazard

Exposure on the hazard of a first event, with a Weibull baseline so the hazard can rise or fall with time. Most dropout and safety analyses depend on the shape of the survival curve.

The data, before any model

Estimates against the truth

Fitted by exact maximum likelihood (no random effects to integrate out), objective function 5047.0, 0.0 s. 3 of 3 parameters have a 95% confidence interval containing the value the data was simulated from.

There's no goodness-of-fit panel or predictive check for this model. Both rely on residuals, and a likelihood model doesn't have any: the data is a zero or a one, so there's nothing to subtract a prediction from. The figure above does that job instead, comparing observed outcomes grouped by exposure with what the model predicts.

Against the simpler model

Does the hazard change with time, or is a single rate enough?

Refitting as a constant hazard (exponential) and re-estimating everything else: dropping it costs 20.5 objective function on 1 degree of freedom, p < 0.0001.

The control stream

tte_weibull.mod passes every static check

$PROBLEM Time to event, Weibull baseline hazard with an exposure effect
; Data: data/tte_weibull.csv, simulated by nmlib.simulate.simulate_tte_weibull
;
; h(t) = LAMBDA*SHAPE*(LAMBDA*t)^(SHAPE-1) * exp(BETA*EXPO)
;
; Fitted on the likelihood directly, not on a residual: the contribution of
; a subject is the survivor function at their last time if they were
; censored, and survivor times hazard if they had the event. SHAPE above 1
; means the hazard rises with time, below 1 that it falls.
;
; Data layout: one record at TIME=0 opening the interval with DV=0, then one
; record at the event or censoring time carrying DV=1 or DV=0.
;
; The frailty term ETA(1) is written in but fixed at zero. With one event
; per subject there is nothing in the data to separate a subject who was
; always at high risk from one who happened to have an early event, so an
; estimated frailty variance here is not identified by the design; the data
; is simulated without one for the same reason. The term stays in the
; stream because it is where a frailty would go given repeated events, and
; freeing it is a one-character change.

$INPUT ID TIME DV EVID MDV EXPO
$DATA ../data/tte_weibull.csv IGNORE=@

$SUBROUTINE ADVAN13 TOL=9

$MODEL
COMP=(CUMHAZ)

$PK
LAMBDA = THETA(1)*EXP(ETA(1))
SHAPE  = THETA(2)
BETA   = THETA(3)

; A(1) accumulates the cumulative hazard from time zero.
A_0(1) = 0

$DES
; DEL keeps the power term finite at T=0 when SHAPE < 1.
DEL     = 1E-12
HAZ     = LAMBDA*SHAPE*(LAMBDA*(T+DEL))**(SHAPE-1)*EXP(BETA*EXPO)
DADT(1) = HAZ

$ERROR
DEL    = 1E-12
CUMHAZ = A(1)
SURV   = EXP(-CUMHAZ)
HAZNOW = LAMBDA*SHAPE*(LAMBDA*(TIME+DEL))**(SHAPE-1)*EXP(BETA*EXPO)

; Censored record: probability of surviving to here.
; Event record: density, survivor times the hazard at the event time.
IF (DV.EQ.0) Y = SURV
IF (DV.EQ.1) Y = SURV*HAZNOW

$THETA
(0, 0.0025)       ; 1 LAMBDA scale (1/day)
(0, 1.35)         ; 2 SHAPE  Weibull shape, >1 means hazard rises with time
(-3, -0.025, 3)   ; 3 BETA   log hazard ratio per unit exposure

$OMEGA
0 FIX             ; 1 frailty on LAMBDA, held at zero -- see below

$ESTIMATION METHOD=1 LAPLACE LIKELIHOOD MAXEVAL=9999 SIG=3 PRINT=5 NOABORT
$COVARIANCE UNCONDITIONAL
$TABLE ID TIME DV EXPO CUMHAZ SURV
       ONEHEADER NOPRINT FILE=sdtab_tte

Exposure-response

Binary response, logistic on exposure

The shape most exposure-safety and exposure-response analyses take when the endpoint is yes or no, with between-subject variability on the logit.

The data, before any model

Estimates against the truth

Fitted by adaptive Gauss-Hermite quadrature, 15 nodes, objective function 2515.9, 0.5 s. 3 of 3 parameters have a 95% confidence interval containing the value the data was simulated from.

There's no goodness-of-fit panel or predictive check for this model. Both rely on residuals, and a likelihood model doesn't have any: the data is a zero or a one, so there's nothing to subtract a prediction from. The figure above does that job instead, comparing observed outcomes grouped by exposure with what the model predicts.

Against the simpler model

Six visits per subject make this estimable. Is there really variation between subjects?

Refitting as plain logistic regression, every subject alike and re-estimating everything else: dropping it costs 40.6 objective function on 1 degree of freedom, p < 0.0001. The null is on the boundary of the parameter space, so the chi-square reference is conservative and the actual p-value is smaller than the one quoted.

The control stream

logistic_binary.mod passes every static check

$PROBLEM Binary response, logistic regression on exposure with IIV
; Data: data/logistic_binary.csv, simulated by nmlib.simulate.simulate_logistic
;
; logit(P) = BASE + SLOPE*EXPO + ETA
;
; Fitted on the likelihood: Y is the probability of the observed outcome,
; P for a responder and 1-P otherwise. This is the shape most exposure-safety
; and exposure-response analyses take when the endpoint is yes or no.
; LAPLACE is required because the likelihood is not normal.

$INPUT ID TIME DV MDV EVID EXPO
$DATA ../data/logistic_binary.csv IGNORE=@

$PRED
BASE  = THETA(1)
SLOPE = THETA(2)

LOGIT = BASE + SLOPE*EXPO + ETA(1)
P     = EXP(LOGIT)/(1 + EXP(LOGIT))

; Contribution of this record to the likelihood.
IF (DV.EQ.1) Y = P
IF (DV.EQ.0) Y = 1 - P

$THETA
(-10, -1.20, 10)  ; 1 BASE  log odds at zero exposure
(-1, 0.055, 1)    ; 2 SLOPE change in log odds per unit exposure

$OMEGA
0.36              ; 1 IIV on the logit (SD 0.6)

$ESTIMATION METHOD=1 LAPLACE LIKELIHOOD MAXEVAL=9999 SIG=3 PRINT=5 NOABORT
$COVARIANCE UNCONDITIONAL
$TABLE ID EXPO DV P ETA(1)
       ONEHEADER NOPRINT FILE=sdtab_logistic

Reference collection

88 published models, catalogued

Control streams collected from published analyses and filed by topic. The five models above were written for this library; these are other people’s and are kept in a private repository. The table shows what can be read off each stream automatically (structure, estimation method, techniques used and what the static checks find), so you can search the collection by technique without opening any files. Titles link to the stream if you have access.

Techniques in use

Detected from the code rather than the file name, so a technique written in an unusual way may be missed. Select one to list the streams that use it.

• Log-transformed both sides30
• Residual error via THETA (SIGMA 1 FIX)29
• Allometric scaling23
• Emax / sigmoid exposure-response23
• Transit compartments20
• Inter-occasion variability14
• Count / categorical likelihood13
• M3 / BLQ likelihood11
• Time to event8
• Indirect response / turnover5
• Markov dependence5
• Target-mediated disposition4
• Age maturation (Hill)3
• Box-Cox transformed ETA3
• Effect compartment3
• Michaelis-Menten elimination3
• $PRIOR2
• Circadian rhythm2
• Mixture model2

Browse

88 of 88 shown. Parameters are THETA / ETA / EPS as declared. 5 streams carry a static-check flag, each a compartment declared in $MODEL that no equation in $DES assigns. That's usually a leftover, and NONMEM runs it without complaint.