# babelmixr2, nlmixr2 and Monolix

By Matt Fidler and the nlmixr2 Development Team in babelmixr2

December 5, 2022

As with `NONMEM`

, it is important to be able to compare `nlmixr2`

to
industry standard software like Monolix. With that , in mind, I am
proud to announce the first `nlmixr2`

to `Monolix`

translator in
`babelmixr2`

.

As with NONMEM, while this has been done before, the *method* whereby we are converting
between the two is novel and has some surprising advantages.

# How to use `Monolix`

with `nlmixr2`

To use `Monolix`

in nlmixr, you do not need to change your data or your
`nlmixr2`

dataset. `babelmixr2`

will do the heavy lifting here.

Lets take the classic warfarin example to start the comparison with…

The model we use in the `nlmixr2`

vignettes is:

```
library(babelmixr2)
# The nonmem translation requires the package pmxTools as well.
# You do not need to load it, simply have it available for use.
pk.turnover.emax3 <- function() {
ini({
tktr <- log(1)
tka <- log(1)
tcl <- log(0.1)
tv <- log(10)
##
eta.ktr ~ 1
eta.ka ~ 1
eta.cl ~ 2
eta.v ~ 1
prop.err <- 0.1
pkadd.err <- 0.1
##
temax <- logit(0.8)
tec50 <- log(0.5)
tkout <- log(0.05)
te0 <- log(100)
##
eta.emax ~ .5
eta.ec50 ~ .5
eta.kout ~ .5
eta.e0 ~ .5
##
pdadd.err <- 10
})
model({
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)
##
DCP = center/v
PD=1-emax*DCP/(ec50+DCP)
##
effect(0) = e0
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
cp ~ prop(prop.err) + add(pkadd.err)
effect ~ add(pdadd.err) | pca
})
}
```

Next you have to either figure out the command to run `Monolix`

, or
simply install the `lixoftConnectors`

package
(https://monolix.lixoft.com/monolix-api/lixoftconnectors_installation/).
If the `lixoftConnectors`

package is available, it will run Monolix
from R using that package.

If the way you run Monolix is a command, you can also set that command for the session by `options("babelmixr2.monolix"="monolixRunCommand")`

or you can use it directly based on
`monolixControl(runCommand="nmfe743")`

. This could also be a
function if you prefer (but I will not cover using the function here).

Lets assume you have `lixoftConnectors`

or something similar setup appropriately. To run the
`nlmixr2`

model using `Monolix`

you simply can run it directly:

```
fit <- nlmixr(pk.turnover.emax3, nlmixr2data::warfarin, "monolix",
monolixControl(modelName="pk.turnover.emax3"))
```

`## ℹ assuming monolix is running because 'pk.turnover.emax3-monolix.txt' is present`

`## → 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...`

`## ✔ done`

`## → Calculating residuals/tables`

`## ✔ done`

`## → compress origData in nlmixr2 object, save 27560`

`## ℹ monolix parameter history needs expoted charts, please export charts`

This fit issues an informational tidbit -

- monolix parameter history needs expoted charts, please export charts

This will automatically be generated as well when `lixoftConnectors`

package is generated and you have a recent version of Monolix. If you
don’t have that information then the important parameter history plots will
not be imported and you cannot see those plots.

Just like with the NONMEM translation, the `monolixControl()`

has
`modelName`

which helps control the output directory of `Monolix`

(if
not specified `babelmixr2`

tries to guess based on the model name
based on the input).

Printing this out this `nlmixr2`

fit you see:

`fit`

Of particular interest is the comparison between Monolix predictions and nlmixr predictions:

```
IPRED relative difference compared to Monolix IPRED: 0.09%; 95% percentile: (0.01%,0.49%); rtol=0.000941
PRED relative difference compared to Monolix PRED: 0.04%; 95% percentile: (0%,0.2%); rtol=0.000428
IPRED absolute difference compared to Monolix IPRED: atol=0.00911; 95% percentile: (0.000493, 0.0928)
PRED absolute difference compared to Monolix PRED: atol=0.000428; 95% percentile: (3.14e-07, 0.203)
```

In this case, I believe that these also imply the models are
predicting the same thing. Note that the model predictions are not as
close as they were with NONMEM because Monolix does not use the
`lsoda`

ODE solver. Hence this small deviation is expected, but still
gives a validated Monolix model.

As in the case of `NONMEM`

, this gives some things that are not
available to Monolix, like adding conditional weighted residuals:

`fit <- addCwres(fit)`

`## → Calculating residuals/tables`

`## ✔ done`

Which will add nlmixr’s CWRES as well as adding the nlmixr2 `FOCEi`

objective function

```
OBJF AIC BIC Log-likelihood Condition Number
FOCEi 1335.312 2261.007 2340.427 -1111.503 2203.836
monolix 1522.704 2448.398 2527.819 -1205.199 2203.836
```

Because you now have an objective function compared based on the same assumptions, you could compare the performance of Monolix and NONMEM based on objective function.

To be fair, objective function values must always be used with caution. How the model performs and predicts the data is far more valuable.

Also since it is a `nlmixr2`

object it would be easy to perform a VPC
too (the same is true for NONMEM models):

```
v1s <- vpcPlot(fit, show=list(obs_dv=TRUE), scales="free_y") +
ylab("Warfarin Cp [mg/L] or PCA") +
xlab("Time [h]")
v2s <- vpcPlot(fit, show=list(obs_dv=TRUE), pred_corr = TRUE, scales="free_y") +
ylab("Prediction Corrected Warfarin Cp [mg/L] or PCA") +
xlab("Time [h]")
v1s
```

`v2s`

## Note about data

The input dataset expected to be compatible with `rxode2`

or
`nlmixr2`

. This dataset is then converted to Monolix format:

The combination of

`CMT`

and Dose type creates a unique`ADM`

variable.The

`ADM`

definition is saved in the monolix model file`babelmixr2`

creates a macro describing the compartment, ie`compartment(cmt=#, amount=stateName)`

`babelmixr2`

also creates a macro for each type of dosing:Bolus/infusion uses

`depot()`

and adds modeled lag time (`Tlag`

) or bioavailibilty (`p`

) if specifiedModeled rate uses

`depot()`

with`Tk0=amtDose/rate`

.`babelmixr2`

also adds modeled lag time (`Tlag`

) or bioavailibilty (`p`

) if specifiedModeled duration uses

`depot()`

with`Tk0=dur`

, also add adds modeled lag time (`Tlag`

) or bioavailibilty (`p`

) if specified Turning off a compartment uses empty macro

- Posted on:
- December 5, 2022

- Length:
- 5 minute read, 1000 words

- Categories:
- babelmixr2

- Tags:
- new-version Monolix