One of the things that can be helpful is reading NONMEM with rounding errors.
For this example we will take some of the test examples from
babelmixr2
to translate to show examples of rounding errors
that can be read in to help diagnose what could be happening with the
model.
library(nonmem2rx)
library(babelmixr2)
#> Loading required package: nlmixr2
#> Loading required package: nlmixr2data
Step 1: Have a NONMEM model with rounding errors (and
.phi
or other information about the etas
)
The first step is to load a model with rounding errors using
nonmem2rx()
:
# Unzip example with rounding error
# included, but can be accessed with nlmixr2
#
# unzip(system.file("tests/testthat/pk.turnover.emax3-nonmem.zip", package="babelmixr2"))
# Load the model with `nonmem2rx`:
mod <- nonmem2rx("pk.turnover.emax3-nonmem/pk.turnover.emax3.nmctl")
#> ℹ getting information from 'pk.turnover.emax3-nonmem/pk.turnover.emax3.nmctl'
#> ℹ reading in xml file
#> ℹ done
#> ℹ reading in phi file
#> ℹ done
#> ℹ reading in lst file
#> ℹ abbreviated list parsing
#> ℹ done
#> ℹ done
#> ℹ splitting control stream by records
#> ℹ done
#> ℹ Processing record $INPUT
#> ℹ Processing record $MODEL
#> ℹ Processing record $THETA
#> ℹ Processing record $OMEGA
#> ℹ Processing record $SIGMA
#> ℹ Processing record $PROBLEM
#> ℹ Processing record $DATA
#> ℹ Processing record $SUBROUTINES
#> ℹ Processing record $PK
#> ℹ Processing record $DES
#> ℹ Processing record $ERROR
#> ℹ Processing record $ESTIMATION
#> ℹ Ignore record $ESTIMATION
#> ℹ Processing record $COVARIANCE
#> ℹ Ignore record $COVARIANCE
#> ℹ Processing record $TABLE
#> ℹ change initial estimate of `theta1` to `6.24053043162953e-07`
#> ℹ change initial estimate of `theta2` to `-3.00642760553675e-06`
#> ℹ change initial estimate of `theta3` to `-2.00405074386117`
#> ℹ change initial estimate of `theta4` to `2.05188410700476`
#> ℹ change initial estimate of `theta5` to `0.0985804613565218`
#> ℹ change initial estimate of `theta6` to `0.511625249037084`
#> ℹ change initial estimate of `theta7` to `6.4184983102259`
#> ℹ change initial estimate of `theta8` to `0.140763261319656`
#> ℹ change initial estimate of `theta9` to `-2.9534704318737`
#> ℹ change initial estimate of `theta10` to `4.57045413136592`
#> ℹ change initial estimate of `theta11` to `3.71714384851537`
#> ℹ change initial estimate of `eta1` to `0.558129815059436`
#> ℹ change initial estimate of `eta2` to `0.558402321309217`
#> ℹ change initial estimate of `eta3` to `0.0785849119252598`
#> ℹ change initial estimate of `eta4` to `0.0508226905750953`
#> ℹ change initial estimate of `eta5` to `5e-05`
#> ℹ change initial estimate of `eta6` to `0.18426809257979`
#> ℹ change initial estimate of `eta7` to `0.0083631531443303`
#> ℹ change initial estimate of `eta8` to `0.00274561514766752`
#> ℹ read in nonmem input data (for model validation): /home/runner/work/nonmem2rx/nonmem2rx/vignettes/articles/pk.turnover.emax3-nonmem/pk.turnover.emax3.csv
#> ℹ ignoring lines that begin with a letter (IGNORE=@)'
#> ℹ applying names specified by $INPUT
#> ℹ renaming 'dvid' to 'nmdvid'
#> ℹ done
#> using C compiler: ‘gcc (Ubuntu 11.3.0-1ubuntu1~22.04.1) 11.3.0’
#> ℹ read in nonmem IPRED data (for model validation): /home/runner/work/nonmem2rx/nonmem2rx/vignettes/articles/pk.turnover.emax3-nonmem/pk.turnover.emax3.pred
#> ℹ done
#> ℹ changing most variables to lower case
#> ℹ done
#> ℹ replace theta names
#> ℹ done
#> ℹ replace eta names
#> ℹ done
#> ℹ renaming compartments
#> ℹ done
#> using C compiler: ‘gcc (Ubuntu 11.3.0-1ubuntu1~22.04.1) 11.3.0’
#> ℹ solving ipred problem
#> ℹ done
#> ℹ solving pred problem
#> ℹ done
mod
#> ── rxode2-based free-form 4-cmt ODE model ──────────────────────────────────────
#> ── Initalization: ──
#> Fixed Effects ($theta):
#> tktr tka tcl tv prop.err
#> 6.240530e-07 -3.006428e-06 -2.004051e+00 2.051884e+00 9.858046e-02
#> pkadd.err temax tec50 tkout te0
#> 5.116252e-01 6.418498e+00 1.407633e-01 -2.953470e+00 4.570454e+00
#> pdadd.err
#> 3.717144e+00
#>
#> Omega ($omega):
#> eta.ktr eta.ka eta.cl eta.v eta.emax eta.ec50
#> eta.ktr 0.5581298 0.0000000 0.00000000 0.00000000 0e+00 0.0000000
#> eta.ka 0.0000000 0.5584023 0.00000000 0.00000000 0e+00 0.0000000
#> eta.cl 0.0000000 0.0000000 0.07858491 0.00000000 0e+00 0.0000000
#> eta.v 0.0000000 0.0000000 0.00000000 0.05082269 0e+00 0.0000000
#> eta.emax 0.0000000 0.0000000 0.00000000 0.00000000 5e-05 0.0000000
#> eta.ec50 0.0000000 0.0000000 0.00000000 0.00000000 0e+00 0.1842681
#> eta.kout 0.0000000 0.0000000 0.00000000 0.00000000 0e+00 0.0000000
#> eta.e0 0.0000000 0.0000000 0.00000000 0.00000000 0e+00 0.0000000
#> eta.kout eta.e0
#> eta.ktr 0.000000000 0.000000000
#> eta.ka 0.000000000 0.000000000
#> eta.cl 0.000000000 0.000000000
#> eta.v 0.000000000 0.000000000
#> eta.emax 0.000000000 0.000000000
#> eta.ec50 0.000000000 0.000000000
#> eta.kout 0.008363153 0.000000000
#> eta.e0 0.000000000 0.002745615
#>
#> States ($state or $stateDf):
#> Compartment Number Compartment Name
#> 1 1 DEPOT
#> 2 2 GUT
#> 3 3 CENTER
#> 4 4 EFFECT
#> ── Model (Normalized Syntax): ──
#> function() {
#> description <- c("translated from babelmixr2", "; comments show mu referenced model in ui$getSplitMuModel")
#> validation <- c("IPRED relative difference compared to Nonmem IPRED: 0%; 95% percentile: (0%,0%); rtol=6.13e-06",
#> "IPRED absolute difference compared to Nonmem IPRED: 95% percentile: (3.12e-06, 0.000497); atol=6.17e-05",
#> "PRED relative difference compared to Nonmem PRED: 0%; 95% percentile: (0%,0%); rtol=6.18e-06",
#> "PRED absolute difference compared to Nonmem PRED: 95% percentile: (3.79e-07,0.00313) atol=6.18e-06")
#> ini({
#> tktr <- 6.24053043162953e-07
#> label("1 - tktr")
#> tka <- -3.00642760553675e-06
#> label("2 - tka")
#> tcl <- -2.00405074386117
#> label("3 - tcl")
#> tv <- 2.05188410700476
#> label("4 - tv")
#> prop.err <- c(0, 0.0985804613565218)
#> label("5 - prop.err")
#> pkadd.err <- c(0, 0.511625249037084)
#> label("6 - pkadd.err")
#> temax <- 6.4184983102259
#> label("7 - temax")
#> tec50 <- 0.140763261319656
#> label("8 - tec50")
#> tkout <- -2.9534704318737
#> label("9 - tkout")
#> te0 <- 4.57045413136592
#> label("10 - te0")
#> pdadd.err <- c(0, 3.71714384851537)
#> label("11 - pdadd.err")
#> eta.ktr ~ 0.558129815059436
#> eta.ka ~ 0.558402321309217
#> eta.cl ~ 0.0785849119252598
#> eta.v ~ 0.0508226905750953
#> eta.emax ~ 5e-05
#> eta.ec50 ~ 0.18426809257979
#> eta.kout ~ 0.0083631531443303
#> eta.e0 ~ 0.00274561514766752
#> })
#> model({
#> cmt(DEPOT)
#> cmt(GUT)
#> cmt(CENTER)
#> cmt(EFFECT)
#> mu_1 <- tktr
#> mu_2 <- tka
#> mu_3 <- tcl
#> mu_4 <- tv
#> mu_5 <- temax
#> mu_6 <- tec50
#> mu_7 <- tkout
#> mu_8 <- te0
#> ktr <- exp(mu_1 + eta.ktr)
#> ka <- exp(mu_2 + eta.ka)
#> cl <- exp(mu_3 + eta.cl)
#> v <- exp(mu_4 + eta.v)
#> emax <- ((1) - (0)) * (1/(1 + exp(-(mu_5 + eta.emax)))) +
#> (0)
#> ec50 <- exp(mu_6 + eta.ec50)
#> kout <- exp(mu_7 + eta.kout)
#> e0 <- exp(mu_8 + eta.e0)
#> rxini.rxddta4. <- e0
#> EFFECT(0) <- rxini.rxddta4.
#> dcp <- CENTER/v
#> rxdz001 <- (ec50 + dcp)
#> if (rxdz001 >= 0 && rxdz001 <= 1e-06) {
#> rxdz001 <- 1e-06
#> }
#> if (rxdz001 >= -1e-06 && rxdz001 < 0) {
#> rxdz001 <- -1e-06
#> }
#> pd <- 1 - emax * dcp/rxdz001
#> kin <- e0 * kout
#> d/dt(DEPOT) <- -ktr * DEPOT
#> d/dt(GUT) <- ktr * DEPOT - ka * GUT
#> d/dt(CENTER) <- ka * GUT - cl/v * CENTER
#> d/dt(EFFECT) <- kin * pd - kout * EFFECT
#> cp <- CENTER/v
#> f <- DEPOT
#> rxe_dcp <- CENTER/v
#> rxdze001 <- (ec50 + rxe_dcp)
#> if (rxdze001 >= 0 && rxdze001 <= 1e-06) {
#> rxdze001 <- 1e-06
#> }
#> if (rxdze001 >= -1e-06 && rxdze001 < 0) {
#> rxdze001 <- -1e-06
#> }
#> rxe_pd <- 1 - emax * rxe_dcp/rxdze001
#> rxe_kin <- e0 * kout
#> rxe_cp <- CENTER/v
#> rx_pf1 <- rxe_cp
#> rx_pf2 <- EFFECT
#> rx_ip1 <- rx_pf1
#> rx_p1 <- rx_ip1
#> w1 <- sqrt((pkadd.err)^2 + (rx_pf1)^2 * (prop.err)^2)
#> if (w1 == 0)
#> w1 <- 1
#> rx_ip2 <- rx_pf2
#> rx_p2 <- rx_ip2
#> w2 <- sqrt((pdadd.err)^2)
#> if (w2 == 0)
#> w2 <- 1
#> ipred <- rx_ip1
#> w <- w1
#> if (nmdvid == 2) {
#> ipred <- rx_ip2
#> w <- w2
#> }
#> y <- ipred + w * eps1
#> })
#> }
#> ── nonmem2rx translation notes ($notes): ──
#> • some etas defaulted to non-mu referenced, possible parsing error: eta.emax as a work-around try putting the mu-referenced expression on a simple line
#> • some etas defaulted to non-mu referenced, possible parsing error: eta5 as a work-around try putting the mu-referenced expression on a simple line
#> • some NONMEM input has tied times; they are offset by a small offset
#> • is.na() applied to non-(list or vector) of type 'language'
#> • 'dvid' variable has special meaning in rxode2, renamed to 'nmdvid', rename/copy in your data too
#> • $MODEL NCOMPARTMENTS/NEQUILIBRIUM/NPARAMETERS statement(s) ignored
#> ── nonmem2rx extra properties: ──
#>
#> Sigma ($sigma):
#> eps1
#> eps1 1
#>
#> other properties include: $nonmemData, $etaData, $dfSub, $dfObs
#> captured NONMEM table outputs: $predData, $ipredData
#> NONMEM/rxode2 comparison data: $iwresCompare, $predCompare, $ipredCompare
#> NONMEM/rxode2 composite comparison: $predAtol, $predRtol, $ipredAtol, $ipredRtol, $iwresAtol, $iwresRtol
Step 2: Convert rxode2
model to model with
endpoints/residuals specified like nlmixr2
This example has rounding errors and isn’t a fully qualified
nlmixr2
model (even though it was generated from
nlmixr2
).
We can use the model and create an equivalent model with
as.nonmem2rx()
mod2 <- function() {
ini({
tktr <- 6.24053043162953e-07
label("1 - tktr")
tka <- -3.00642760553675e-06
label("2 - tka")
tcl <- -2.00405074386117
label("3 - tcl")
tv <- 2.05188410700476
label("4 - tv")
prop.err <- c(0, 0.0985804613565218)
label("5 - prop.err")
pkadd.err <- c(0, 0.511625249037084)
label("6 - pkadd.err")
temax <- 6.4184983102259
label("7 - temax")
tec50 <- 0.140763261319656
label("8 - tec50")
tkout <- -2.9534704318737
label("9 - tkout")
te0 <- 4.57045413136592
label("10 - te0")
pdadd.err <- c(0, 3.71714384851537)
label("11 - pdadd.err")
eta.ktr ~ 0.558129815059436
eta.ka ~ 0.558402321309217
eta.cl ~ 0.0785849119252598
eta.v ~ 0.0508226905750953
eta.emax ~ 5e-05
eta.ec50 ~ 0.18426809257979
eta.kout ~ 0.0083631531443303
eta.e0 ~ 0.00274561514766752
})
model({
cmt(DEPOT)
cmt(GUT)
cmt(CENTER)
cmt(EFFECT)
ktr <- exp(tktr + eta.ktr)
ka <- exp(tka + eta.ka)
cl <- exp(tcl + eta.cl)
v <- exp(tv + eta.v)
emax <- expit(temax + eta.emax)
ec50 <- exp(tec50 + eta.ec50)
kout <- exp(tkout + eta.kout)
e0 <- exp(te0 + eta.e0)
EFFECT(0) <- e0
dcp <- CENTER/v
pd <- 1 - emax * dcp/(ec50 + dcp)
kin <- e0 * kout
d/dt(DEPOT) <- -ktr * DEPOT
d/dt(GUT) <- ktr * DEPOT - ka * GUT
d/dt(CENTER) <- ka * GUT - cl/v * CENTER
d/dt(EFFECT) <- kin * pd - kout * EFFECT
eff <- EFFECT
dcp ~ add(pkadd.err)+prop(prop.err)
eff ~ add(pdadd.err)
})
}
In the model above, the following modifications were made:
The code for protecting for rounding errors was removed
Endpoints/residual specifications were added
Duplicate code in NONMEM’s
$ERROR
block was removed from imported modelAll the
rx
andnlmixr
prefixed items were removed from the model. It is possible that not removing these variables can causenlmixr2
conversion to fail. It is best practice to remove them completely.
Also, it is good practice to make sure the model parses correctly before trying to validate/convert the model. If you find that your model doesn’t parse correctly it definitely won’t validate (and the error may not be as easy to track-down).
The first step is to validate if the translation is correct. This is done below:
new <- as.nonmem2rx(mod, mod2)
#> ℹ merging 'dvid' with nlmixr2 'cmt' definition
#> using C compiler: ‘gcc (Ubuntu 11.3.0-1ubuntu1~22.04.1) 11.3.0’
#> ℹ solving ipred problem
#> ℹ done
#> ℹ solving pred problem
#> ℹ done
print(new)
#> ── rxode2-based free-form 4-cmt ODE model ──────────────────────────────────────
#> ── Initalization: ──
#> Fixed Effects ($theta):
#> tktr tka tcl tv prop.err
#> 6.240530e-07 -3.006428e-06 -2.004051e+00 2.051884e+00 9.858046e-02
#> pkadd.err temax tec50 tkout te0
#> 5.116252e-01 6.418498e+00 1.407633e-01 -2.953470e+00 4.570454e+00
#> pdadd.err
#> 3.717144e+00
#>
#> Omega ($omega):
#> eta.ktr eta.ka eta.cl eta.v eta.emax eta.ec50
#> eta.ktr 0.5581298 0.0000000 0.00000000 0.00000000 0e+00 0.0000000
#> eta.ka 0.0000000 0.5584023 0.00000000 0.00000000 0e+00 0.0000000
#> eta.cl 0.0000000 0.0000000 0.07858491 0.00000000 0e+00 0.0000000
#> eta.v 0.0000000 0.0000000 0.00000000 0.05082269 0e+00 0.0000000
#> eta.emax 0.0000000 0.0000000 0.00000000 0.00000000 5e-05 0.0000000
#> eta.ec50 0.0000000 0.0000000 0.00000000 0.00000000 0e+00 0.1842681
#> eta.kout 0.0000000 0.0000000 0.00000000 0.00000000 0e+00 0.0000000
#> eta.e0 0.0000000 0.0000000 0.00000000 0.00000000 0e+00 0.0000000
#> eta.kout eta.e0
#> eta.ktr 0.000000000 0.000000000
#> eta.ka 0.000000000 0.000000000
#> eta.cl 0.000000000 0.000000000
#> eta.v 0.000000000 0.000000000
#> eta.emax 0.000000000 0.000000000
#> eta.ec50 0.000000000 0.000000000
#> eta.kout 0.008363153 0.000000000
#> eta.e0 0.000000000 0.002745615
#>
#> States ($state or $stateDf):
#> Compartment Number Compartment Name
#> 1 1 DEPOT
#> 2 2 GUT
#> 3 3 CENTER
#> 4 4 EFFECT
#> ── Multiple Endpoint Model ($multipleEndpoint): ──
#> variable cmt dvid*
#> 1 dcp ~ … cmt='dcp' or cmt=5 dvid='dcp' or dvid=1
#> 2 eff ~ … cmt='eff' or cmt=6 dvid='eff' or dvid=2
#> * If dvids are outside this range, all dvids are re-numered sequentially, ie 1,7, 10 becomes 1,2,3 etc
#>
#> ── μ-referencing ($muRefTable): ──
#> theta eta level
#> 1 tktr eta.ktr id
#> 2 tka eta.ka id
#> 3 tcl eta.cl id
#> 4 tv eta.v id
#> 5 temax eta.emax id
#> 6 tec50 eta.ec50 id
#> 7 tkout eta.kout id
#> 8 te0 eta.e0 id
#>
#> ── Model (Normalized Syntax): ──
#> function() {
#> description <- c("translated from babelmixr2", "; comments show mu referenced model in ui$getSplitMuModel")
#> validation <- c("IPRED relative difference compared to Nonmem IPRED: 0%; 95% percentile: (0%,0%); rtol=6.13e-06",
#> "IPRED absolute difference compared to Nonmem IPRED: 95% percentile: (3.12e-06, 0.000497); atol=6.17e-05",
#> "PRED relative difference compared to Nonmem PRED: 0%; 95% percentile: (0%,0%); rtol=6.18e-06",
#> "PRED absolute difference compared to Nonmem PRED: 95% percentile: (3.79e-07,0.00313) atol=6.18e-06")
#> ini({
#> tktr <- 6.24053043162953e-07
#> label("1 - tktr")
#> tka <- -3.00642760553675e-06
#> label("2 - tka")
#> tcl <- -2.00405074386117
#> label("3 - tcl")
#> tv <- 2.05188410700476
#> label("4 - tv")
#> prop.err <- c(0, 0.0985804613565218)
#> label("5 - prop.err")
#> pkadd.err <- c(0, 0.511625249037084)
#> label("6 - pkadd.err")
#> temax <- 6.4184983102259
#> label("7 - temax")
#> tec50 <- 0.140763261319656
#> label("8 - tec50")
#> tkout <- -2.9534704318737
#> label("9 - tkout")
#> te0 <- 4.57045413136592
#> label("10 - te0")
#> pdadd.err <- c(0, 3.71714384851537)
#> label("11 - pdadd.err")
#> eta.ktr ~ 0.558129815059436
#> eta.ka ~ 0.558402321309217
#> eta.cl ~ 0.0785849119252598
#> eta.v ~ 0.0508226905750953
#> eta.emax ~ 5e-05
#> eta.ec50 ~ 0.18426809257979
#> eta.kout ~ 0.0083631531443303
#> eta.e0 ~ 0.00274561514766752
#> })
#> model({
#> cmt(DEPOT)
#> cmt(GUT)
#> cmt(CENTER)
#> cmt(EFFECT)
#> ktr <- exp(tktr + eta.ktr)
#> ka <- exp(tka + eta.ka)
#> cl <- exp(tcl + eta.cl)
#> v <- exp(tv + eta.v)
#> emax <- expit(temax + eta.emax)
#> ec50 <- exp(tec50 + eta.ec50)
#> kout <- exp(tkout + eta.kout)
#> e0 <- exp(te0 + eta.e0)
#> EFFECT(0) <- e0
#> dcp <- CENTER/v
#> pd <- 1 - emax * dcp/(ec50 + dcp)
#> kin <- e0 * kout
#> d/dt(DEPOT) <- -ktr * DEPOT
#> d/dt(GUT) <- ktr * DEPOT - ka * GUT
#> d/dt(CENTER) <- ka * GUT - cl/v * CENTER
#> d/dt(EFFECT) <- kin * pd - kout * EFFECT
#> eff <- EFFECT
#> dcp ~ add(pkadd.err) + prop(prop.err)
#> eff ~ add(pdadd.err)
#> })
#> }
#> ── nonmem2rx extra properties: ──
#> other properties include: $nonmemData, $etaData, $thetaMat, $dfSub, $dfObs
#> captured NONMEM table outputs: $predData, $ipredData
#> NONMEM/rxode2 comparison data: $iwresCompare, $predCompare, $ipredCompare
#> NONMEM/rxode2 composite comparison: $predAtol, $predRtol, $ipredAtol, $ipredRtol, $iwresAtol, $iwresRtol
Step 3: Convert new model to nlmixr2 fit with
as.nlmixr2()
Once the translation is complete, and it is validated, you can convert the fit to a full nlmixr2 fit object.
fit <- as.nlmixr2(new)
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of full model...
#> ✔ done
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling EBE model...
#> using C compiler: ‘gcc (Ubuntu 11.3.0-1ubuntu1~22.04.1) 11.3.0’
#> ✔ done
#> rxode2 2.0.13.9000 using 1 threads (see ?getRxThreads)
#> no cache: create with `rxCreateCache()`
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 21592
#> → compress parHist in nlmixr2 object, save 4928
# Once it is loaded remove the directory (we don't need the files any
# more for this example)
#
# In this example, we don't remove, but note where it can be removed
#
# unlink("pk.turnover.emax3-nonmem", recursive = TRUE)
Note this object does not rerun an estimation, but
rather imports everything it knows about the fit to rerun the
nlmixr2
table steps to calculate the things you need.
Step 4: Explore the data/fit as if it came from
nlmixr2
With this information you can see about the properties of the model:
print(fit)
#> ── nlmixr² nonmem2rx reading NONMEM ver 7.4.3 ──
#>
#> OBJF AIC BIC Log-likelihood
#> nonmem2rx 439.2156 1364.91 1444.331 -663.4551
#>
#> ── Time (sec $time): ──
#>
#> setup table compress NONMEM as.nlmixr2
#> elapsed 0.025754 0.043 0.014 320.27 3.947
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#>
#> Parameter Est. Back-transformed BSV(CV% or SD) Shrink(SD)%
#> tktr 1 - tktr 6.24e-07 1 86.5 59.8%
#> tka 2 - tka -3.01e-06 1 86.5 59.8%
#> tcl 3 - tcl -2 0.135 28.6 1.34%
#> tv 4 - tv 2.05 7.78 22.8 6.44%
#> prop.err 5 - prop.err 0.0986 0.0986
#> pkadd.err 6 - pkadd.err 0.512 0.512
#> temax 7 - temax 6.42 0.998 3.00 100.%
#> tec50 8 - tec50 0.141 1.15 45.0 6.06%
#> tkout 9 - tkout -2.95 0.0522 9.16 32.4%
#> te0 10 - te0 4.57 96.6 5.24 18.1%
#> pdadd.err 11 - pdadd.err 3.72 3.72
#>
#> No correlations in between subject variability (BSV) matrix
#> Full BSV covariance ($omega) or correlation ($omegaR; diagonals=SDs)
#> Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink
#> Censoring ($censInformation): No censoring
#> Minimization message ($message):
#>
#>
#> WARNINGS AND ERRORS (IF ANY) FOR PROBLEM 1
#>
#> (WARNING 2) NM-TRAN INFERS THAT THE DATA ARE POPULATION.
#>
#>
#> 0MINIMIZATION TERMINATED
#> DUE TO ROUNDING ERRORS (ERROR=134)
#> NO. OF FUNCTION EVALUATIONS USED: 1088
#> NO. OF SIG. DIGITS UNREPORTABLE
#> 0PARAMETER ESTIMATE IS NEAR ITS BOUNDARY
#>
#> IPRED relative difference compared to Nonmem IPRED: 0%; 95% percentile: (0%,0%); rtol=5.09e-06
#> PRED relative difference compared to Nonmem PRED: 0%; 95% percentile: (0%,0%); rtol=5.29e-06
#> IPRED absolute difference compared to Nonmem IPRED: 95% percentile: (2.2e-06, 0.000454); atol=3.03e-05
#> PRED absolute difference compared to Nonmem PRED: 95% percentile: (4.72e-07,0.00361); atol=5.29e-06
#> there are solving errors during optimization (see '$prderr')
#> nonmem2rx model file: 'pk.turnover.emax3-nonmem/pk.turnover.emax3.nmctl'
#>
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 483 × 35
#> ID TIME CMT DV PRED RES IPRED IRES IWRES eta.ktr eta.ka eta.cl
#> <fct> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0.5 dcp 0 1.16 -1.16 0.444 -0.444 -0.864 -0.506 -0.506 0.699
#> 2 1 1 dcp 1.9 3.37 -1.47 1.45 0.446 0.840 -0.506 -0.506 0.699
#> 3 1 2 dcp 3.3 7.51 -4.21 3.96 -0.660 -1.03 -0.506 -0.506 0.699
#> # ℹ 480 more rows
#> # ℹ 23 more variables: eta.v <dbl>, eta.emax <dbl>, eta.ec50 <dbl>,
#> # eta.kout <dbl>, eta.e0 <dbl>, dcp <dbl>, eff <dbl>, DEPOT <dbl>, GUT <dbl>,
#> # CENTER <dbl>, EFFECT <dbl>, ktr <dbl>, ka <dbl>, cl <dbl>, v <dbl>,
#> # emax <dbl>, ec50 <dbl>, kout <dbl>, e0 <dbl>, pd <dbl>, kin <dbl>,
#> # tad <dbl>, dosenum <dbl>
In this model you can see:
- High shrinkage in
temax
,ktr
,ka
and moderate shrinkage inkout
By removing these parameters, you could possibly get a successful NONMEM run.
If you want, you can use model piping to remove these parameters as follows:
mod3 <- fit %>%
model(ktr <- exp(tktr)) %>%
model(ka <- exp(tka)) %>%
model(emax <- expit(temax)) %>%
model(kout <- exp(tkout))
#> ! remove between subject variability `eta.ktr`
#> ! remove between subject variability `eta.ka`
#> ! remove between subject variability `eta.emax`
#> ! remove between subject variability `eta.kout`
mod3
#> ── rxode2-based free-form 4-cmt ODE model ──────────────────────────────────────
#> ── Initalization: ──
#> Fixed Effects ($theta):
#> tktr tka tcl tv prop.err
#> 6.240530e-07 -3.006428e-06 -2.004051e+00 2.051884e+00 9.858046e-02
#> pkadd.err temax tec50 tkout te0
#> 5.116252e-01 6.418498e+00 1.407633e-01 -2.953470e+00 4.570454e+00
#> pdadd.err
#> 3.717144e+00
#>
#> Omega ($omega):
#> eta.cl eta.v eta.ec50 eta.e0
#> eta.cl 0.07858491 0.00000000 0.0000000 0.000000000
#> eta.v 0.00000000 0.05082269 0.0000000 0.000000000
#> eta.ec50 0.00000000 0.00000000 0.1842681 0.000000000
#> eta.e0 0.00000000 0.00000000 0.0000000 0.002745615
#>
#> States ($state or $stateDf):
#> Compartment Number Compartment Name
#> 1 1 DEPOT
#> 2 2 GUT
#> 3 3 CENTER
#> 4 4 EFFECT
#> ── Multiple Endpoint Model ($multipleEndpoint): ──
#> variable cmt dvid*
#> 1 dcp ~ … cmt='dcp' or cmt=5 dvid='dcp' or dvid=1
#> 2 eff ~ … cmt='eff' or cmt=6 dvid='eff' or dvid=2
#> * If dvids are outside this range, all dvids are re-numered sequentially, ie 1,7, 10 becomes 1,2,3 etc
#>
#> ── μ-referencing ($muRefTable): ──
#> theta eta level
#> 1 tcl eta.cl id
#> 2 tv eta.v id
#> 3 tec50 eta.ec50 id
#> 4 te0 eta.e0 id
#>
#> ── Model (Normalized Syntax): ──
#> function() {
#> description <- c("translated from babelmixr2", "; comments show mu referenced model in ui$getSplitMuModel")
#> validation <- c("IPRED relative difference compared to Nonmem IPRED: 0%; 95% percentile: (0%,0%); rtol=6.13e-06",
#> "IPRED absolute difference compared to Nonmem IPRED: 95% percentile: (3.12e-06, 0.000497); atol=6.17e-05",
#> "PRED relative difference compared to Nonmem PRED: 0%; 95% percentile: (0%,0%); rtol=6.18e-06",
#> "PRED absolute difference compared to Nonmem PRED: 95% percentile: (3.79e-07,0.00313) atol=6.18e-06")
#> ini({
#> tktr <- 6.24053043162953e-07
#> label("1 - tktr")
#> tka <- -3.00642760553675e-06
#> label("2 - tka")
#> tcl <- -2.00405074386117
#> label("3 - tcl")
#> tv <- 2.05188410700476
#> label("4 - tv")
#> prop.err <- c(0, 0.0985804613565218)
#> label("5 - prop.err")
#> pkadd.err <- c(0, 0.511625249037084)
#> label("6 - pkadd.err")
#> temax <- 6.4184983102259
#> label("7 - temax")
#> tec50 <- 0.140763261319656
#> label("8 - tec50")
#> tkout <- -2.9534704318737
#> label("9 - tkout")
#> te0 <- 4.57045413136592
#> label("10 - te0")
#> pdadd.err <- c(0, 3.71714384851537)
#> label("11 - pdadd.err")
#> eta.cl ~ 0.0785849119252598
#> eta.v ~ 0.0508226905750953
#> eta.ec50 ~ 0.18426809257979
#> eta.e0 ~ 0.00274561514766752
#> })
#> model({
#> cmt(DEPOT)
#> cmt(GUT)
#> cmt(CENTER)
#> cmt(EFFECT)
#> ktr <- exp(tktr)
#> ka <- exp(tka)
#> cl <- exp(tcl + eta.cl)
#> v <- exp(tv + eta.v)
#> emax <- expit(temax)
#> ec50 <- exp(tec50 + eta.ec50)
#> kout <- exp(tkout)
#> e0 <- exp(te0 + eta.e0)
#> EFFECT(0) <- e0
#> dcp <- CENTER/v
#> pd <- 1 - emax * dcp/(ec50 + dcp)
#> kin <- e0 * kout
#> d/dt(DEPOT) <- -ktr * DEPOT
#> d/dt(GUT) <- ktr * DEPOT - ka * GUT
#> d/dt(CENTER) <- ka * GUT - cl/v * CENTER
#> d/dt(EFFECT) <- kin * pd - kout * EFFECT
#> eff <- EFFECT
#> dcp ~ add(pkadd.err) + prop(prop.err)
#> eff ~ add(pdadd.err)
#> })
#> }
Since you have babelmixr2
you could use
babelmixr2
to run the model for you in NONMEM (if you have
it setup to run NONMEM), or even run the model from nlmixr2
itself. In this example we will choose to use nlmixr2
(since the babelmixr2
example runs NONMEM and shows a
reduction in shrinkage)
fit2 <- nlmixr(mod3, new$nonmemData, "focei", foceiControl(print=0))
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of full model...
#> ✔ done
#> → calculate jacobian
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate sensitivities
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate ∂(f)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate ∂(R²)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling inner model...
#> using C compiler: ‘gcc (Ubuntu 11.3.0-1ubuntu1~22.04.1) 11.3.0’
#> ✔ done
#> → finding duplicate expressions in FD model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in FD model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling EBE model...
#> using C compiler: ‘gcc (Ubuntu 11.3.0-1ubuntu1~22.04.1) 11.3.0’
#> ✔ done
#> → compiling events FD model...
#> using C compiler: ‘gcc (Ubuntu 11.3.0-1ubuntu1~22.04.1) 11.3.0’
#> ✔ done
#> Theta reset (zero/bad gradient values); Switch to bobyqa
#> Error in .foceiFitInternal(.env) : Starting values violate bounds
#> Restart 1
#> [====|====|====|====|====|====|====|====|====|====] 0:06:25
#> done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 21592
#> → compress parHist in nlmixr2 object, save 158656
fit2
#> ── nlmixr² FOCEi (outer: bobyqa) ──
#>
#> OBJF AIC BIC Log-likelihood
#> FOCEi 3180.774 4098.469 4161.169 -2034.234
#>
#> ── Time (sec fit2$time): ──
#>
#> setup optimize covariance table compress other
#> elapsed 0.00962 0.000213 0.000215 0.069 0.014 401.767
#>
#> ── Population Parameters (fit2$parFixed or fit2$parFixedDf): ──
#>
#> Parameter Est. Back-transformed BSV(CV%) Shrink(SD)%
#> tktr 1 - tktr 12 1.65e+05
#> tka 2 - tka 11.9 1.4e+05
#> tcl 3 - tcl -0.561 0.57 2.80 97.5%
#> tv 4 - tv 2.54 12.6 2.25 97.1%
#> prop.err 5 - prop.err 5.96e-07 5.96e-07
#> pkadd.err 6 - pkadd.err 6.42 6.42
#> temax 7 - temax 1.98e+04 1
#> tec50 8 - tec50 -11.9 7.04e-06 4.78 97.5%
#> tkout 9 - tkout -27.6 1.01e-12
#> te0 10 - te0 3.5 33 0.610 99.7%
#> pdadd.err 11 - pdadd.err 44.8 44.8
#>
#> Covariance Type (fit2$covMethod): Boundary issue; Get SEs with `getVarCov()`: "prop.err"
#> No correlations in between subject variability (BSV) matrix
#> Full BSV covariance (fit2$omega) or correlation (fit2$omegaR; diagonals=SDs)
#> Distribution stats (mean/skewness/kurtosis/p-value) available in fit2$shrink
#> Information about run found (fit2$runInfo):
#> • parameter estimate near boundary; covariance not calculated: "prop.err" use 'getVarCov' to calculate anyway
#> • ETAs were reset to zero during optimization; (Can control by foceiControl(resetEtaP=.))
#> • initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> • Hessian reset during optimization; (can control by foceiControl(resetHessianAndEta=.))
#> • thetas were reset during optimization because of a zero gradient
#> Censoring (fit2$censInformation): No censoring
#> Minimization message (fit2$message):
#> Normal exit from bobyqa
#>
#> ── Fit Data (object fit2 is a modified tibble): ──
#> # A tibble: 483 × 33
#> ID TIME CMT DV PRED RES WRES IPRED IRES IWRES CPRED CRES
#> <fct> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0.5 dcp 0 7.75 -7.75 -1.21 7.74 -7.74 -1.21 7.75 -7.75
#> 2 1 1 dcp 1.9 7.58 -5.68 -0.883 7.57 -5.67 -0.883 7.58 -5.68
#> 3 1 2 dcp 3.3 7.24 -3.94 -0.613 7.23 -3.93 -0.612 7.24 -3.94
#> # ℹ 480 more rows
#> # ℹ 21 more variables: CWRES <dbl>, eta.cl <dbl>, eta.v <dbl>, eta.ec50 <dbl>,
#> # eta.e0 <dbl>, DEPOT <dbl>, GUT <dbl>, CENTER <dbl>, EFFECT <dbl>,
#> # ktr <dbl>, ka <dbl>, cl <dbl>, v <dbl>, emax <dbl>, ec50 <dbl>, kout <dbl>,
#> # e0 <dbl>, pd <dbl>, kin <dbl>, tad <dbl>, dosenum <dbl>
With these modifications the shrinkage is also reduced (just like in the NONMEM case)
Step 5: Get the covariance of the model
Another thing that can be helpful when the fit has been imported into
a nlmixr2
fit is to get the variance/covariance matrix.
This can be especially helpful to diagnose things to help simplify your
model
getVarCov(fit)
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of full model...
#> ✔ done
#> → calculate jacobian
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate sensitivities
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate ∂(f)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling inner model...
#> using C compiler: ‘gcc (Ubuntu 11.3.0-1ubuntu1~22.04.1) 11.3.0’
#> ✔ done
#> → finding duplicate expressions in FD model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in FD model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling EBE model...
#> using C compiler: ‘gcc (Ubuntu 11.3.0-1ubuntu1~22.04.1) 11.3.0’
#> ✔ done
#> → compiling events FD model...
#> using C compiler: ‘gcc (Ubuntu 11.3.0-1ubuntu1~22.04.1) 11.3.0’
#> ✔ done
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|=
#> Warning in foceiFitCpp_(.ret): S matrix non-positive definite
#> ===] 0:00:09
#> Warning in foceiFitCpp_(.ret): R matrix non-positive definite but corrected by
#> R = sqrtm(R%*%R)
#> Warning in foceiFitCpp_(.ret): using R matrix to calculate covariance, can
#> check sandwich or S matrix with $covRS and $covS
#> Warning in foceiFitCpp_(.ret): gradient problems with covariance; see
#> $scaleInfo
#> → compress origData in nlmixr2 object, save 21592
#> Updated original fit object fit
#> tktr tka tcl tv temax
#> tktr 0.027507185 -0.0229743954 -0.0118847050 0.0024605889 -0.0015303604
#> tka -0.022974395 0.0845683475 -0.1209728069 -0.0042645578 0.0003075893
#> tcl -0.011884705 -0.1209728069 0.6242221602 -0.0041134945 0.0095313445
#> tv 0.002460589 -0.0042645578 -0.0041134945 0.0108631250 0.0008375679
#> temax -0.001530360 0.0003075893 0.0095313445 0.0008375679 44.4431461692
#> tec50 -0.034910826 0.0225962330 0.1150552735 -0.0012409312 0.0491652327
#> tkout -0.002987253 0.0081951496 0.0041489065 0.0073475010 -0.1357478339
#> te0 -0.004760930 0.0039403914 0.0007457058 -0.0036847801 -0.0262804555
#> tec50 tkout te0
#> tktr -0.034910826 -0.002987253 -0.0047609302
#> tka 0.022596233 0.008195150 0.0039403914
#> tcl 0.115055273 0.004148906 0.0007457058
#> tv -0.001240931 0.007347501 -0.0036847801
#> temax 0.049165233 -0.135747834 -0.0262804555
#> tec50 0.640371952 0.062631972 0.0150208540
#> tkout 0.062631972 0.074705994 -0.0291146671
#> te0 0.015020854 -0.029114667 0.0751101557
fit
#> ── nlmixr² nonmem2rx reading NONMEM ver 7.4.3 ──
#>
#> OBJF AIC BIC Log-likelihood
#> nonmem2rx 439.2156 1364.91 1444.331 -663.4551
#>
#> ── Time (sec fit$time): ──
#>
#> setup table compress NONMEM as.nlmixr2 covariance
#> elapsed 0.025754 0.043 0.014 320.27 3.947 19.5
#>
#> ── Population Parameters (fit$parFixed or fit$parFixedDf): ──
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI)
#> tktr 1 - tktr 6.24e-07 0.166 2.66e+07 1 (0.722, 1.38)
#> tka 2 - tka -3.01e-06 0.291 9.67e+06 1 (0.566, 1.77)
#> tcl 3 - tcl -2 0.79 39.4 0.135 (0.0287, 0.634)
#> tv 4 - tv 2.05 0.104 5.08 7.78 (6.34, 9.55)
#> prop.err 5 - prop.err 0.0986 0.0986
#> pkadd.err 6 - pkadd.err 0.512 0.512
#> temax 7 - temax 6.42 6.67 104 0.998 (0.0013, 1)
#> tec50 8 - tec50 0.141 0.8 568 1.15 (0.24, 5.52)
#> tkout 9 - tkout -2.95 0.273 9.25 0.0522 (0.0305, 0.0891)
#> te0 10 - te0 4.57 0.274 6 96.6 (56.4, 165)
#> pdadd.err 11 - pdadd.err 3.72 3.72
#> BSV(CV% or SD) Shrink(SD)%
#> tktr 86.5 59.8%
#> tka 86.5 59.8%
#> tcl 28.6 1.34%
#> tv 22.8 6.44%
#> prop.err
#> pkadd.err
#> temax 3.00 100.%
#> tec50 45.0 6.06%
#> tkout 9.16 32.4%
#> te0 5.24 18.1%
#> pdadd.err
#>
#> Covariance Type (fit$covMethod): |r|
#> No correlations in between subject variability (BSV) matrix
#> Full BSV covariance (fit$omega) or correlation (fit$omegaR; diagonals=SDs)
#> Distribution stats (mean/skewness/kurtosis/p-value) available in fit$shrink
#> Censoring (fit$censInformation): No censoring
#> Minimization message (fit$message):
#>
#>
#> WARNINGS AND ERRORS (IF ANY) FOR PROBLEM 1
#>
#> (WARNING 2) NM-TRAN INFERS THAT THE DATA ARE POPULATION.
#>
#>
#> 0MINIMIZATION TERMINATED
#> DUE TO ROUNDING ERRORS (ERROR=134)
#> NO. OF FUNCTION EVALUATIONS USED: 1088
#> NO. OF SIG. DIGITS UNREPORTABLE
#> 0PARAMETER ESTIMATE IS NEAR ITS BOUNDARY
#>
#> IPRED relative difference compared to Nonmem IPRED: 0%; 95% percentile: (0%,0%); rtol=5.09e-06
#> PRED relative difference compared to Nonmem PRED: 0%; 95% percentile: (0%,0%); rtol=5.29e-06
#> IPRED absolute difference compared to Nonmem IPRED: 95% percentile: (2.2e-06, 0.000454); atol=3.03e-05
#> PRED absolute difference compared to Nonmem PRED: 95% percentile: (4.72e-07,0.00361); atol=5.29e-06
#> there are solving errors during optimization (see '$prderr')
#> nonmem2rx model file: 'pk.turnover.emax3-nonmem/pk.turnover.emax3.nmctl'
#>
#> ── Fit Data (object fit is a modified tibble): ──
#> # A tibble: 483 × 35
#> ID TIME CMT DV PRED RES IPRED IRES IWRES eta.ktr eta.ka eta.cl
#> <fct> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0.5 dcp 0 1.16 -1.16 0.444 -0.444 -0.864 -0.506 -0.506 0.699
#> 2 1 1 dcp 1.9 3.37 -1.47 1.45 0.446 0.840 -0.506 -0.506 0.699
#> 3 1 2 dcp 3.3 7.51 -4.21 3.96 -0.660 -1.03 -0.506 -0.506 0.699
#> # ℹ 480 more rows
#> # ℹ 23 more variables: eta.v <dbl>, eta.emax <dbl>, eta.ec50 <dbl>,
#> # eta.kout <dbl>, eta.e0 <dbl>, dcp <dbl>, eff <dbl>, DEPOT <dbl>, GUT <dbl>,
#> # CENTER <dbl>, EFFECT <dbl>, ktr <dbl>, ka <dbl>, cl <dbl>, v <dbl>,
#> # emax <dbl>, ec50 <dbl>, kout <dbl>, e0 <dbl>, pd <dbl>, kin <dbl>,
#> # tad <dbl>, dosenum <dbl>
Note this covariance step is not 100% successful since it is not
r, s
. However, it can give some insights on which
parameters are estimated well. In this case you can see that the
emax
parameter is more poorly estimated than other
parameters, which means fixing the parameter or reducing other
parameters may help your estimate progress in NONMEM
.