library(mrgsolve)
library(dplyr)
mod <- modlib("1005")1 mwrite
mrgsolve has used mread() since the start to read, parse, and compile a model. The 1.5.1 release introduces functions to write the resulting model back to a file.
Let’s illustrate with a model from the model library inside the mrgsolve package.
We are using a model that was estimated in NONMEM. So the $NMXML code block is in play; this imports THETAs, OMEGAs, and SIGMAs from the root.xml file generated by NONMEM.
blocks(mod, NMXML)
Model file: 1005.cpp
$NMXML
project = system.file("nonmem", package = "mrgsolve")
run = "@cppstem"The $NMXML (and $NMEXT) blocks are very convenient for making sure that mrgsolve is using the same estimates that NONMEM found during estimation.
But what if we wanted to share this model with a collaborator? We can’t just send this file around with out packaging up the completed NONMEM run as well.
We might also want to make updates to this model. For example,
mod <- param(mod, THETA1 = 8.51, SEX = 1)
mod <- update(
mod,
end = 168,
rtol = 1e-5,
ss_rtol = 1e-4,
outvars = "IPRED"
)This model is now set up for my simulations and I want to save this model back to a file. With mrgsolve 1.5.1, we can write the model back to a file with mwrite_yaml() (yaml structured code) or mwrite_cpp() (native’ mrgsolve format).
mwrite_yaml(mod, file = "1005-updated.mod")
See the code in yaml format
1005-updated.mod
source: mrgsolve::mwrite
mrgsolve: 1.5.2
format: yaml
version: 1.0
model: '1005'
prob:
- 1005 phase1 2 CMT like 1004 but diff. initial on V3
- Run
- file.show(system.file("nonmem", "1005", "1005.ctl", package = "mrgsolve"))
- for equivalent NONMEM control stream.
param:
SEX: 1.0
WT: 70.0
THETA1: 8.51
THETA2: 22.79098949
THETA3: 0.07143366
THETA4: 3.47450589
THETA5: 113.27671224
THETA6: 1.02435433
THETA7: 1.19211818
init:
GUT: 0.0
CENT: 0.0
PERIPH: 0.0
capture:
- CL
- Q
- V2
- V3
- KA
- ETA_1 = ETA(1)
- ETA_2 = ETA(2)
- ETA_3 = ETA(3)
- IPRED
omega:
data:
matrix1:
row1: 0.21387884
row2:
- 0.1207702
- 0.09451047
row3:
- -0.01162777
- -0.03720637
- 0.04656315
labels:
matrix1:
- '.'
- '.'
- '.'
names: '...'
sigma:
data:
matrix1:
row1: 0.04917071
row2:
- 0.0
- 0.20176885
labels:
matrix1:
- '.'
- '.'
names: '...'
envir: []
plugin: base
update:
start: 0.0
end: 168.0
delta: 1.0
add: []
atol: 1.0e-08
rtol: 1.0e-05
ss_atol: 1.0e-08
ss_rtol: 0.0001
maxsteps: 20000.0
hmax: 0.0
hmin: 0.0
mxhnil: 2.0
ixpr: 0.0
mindt: 2.22044605e-15
digits: -1.0
tscale: 1.0
outvars: IPRED
set: []
code:
- $PKMODEL
- depot = TRUE
- ncmt = 2
- ' '
- $PK
- double CL = THETA(1)*exp(ETA(1)) * pow(THETA(6),SEX) * pow(WT/70.0,THETA(7));
- double V2 = THETA(2)*exp(ETA(2));
- double KA = THETA(3)*exp(ETA(3));
- double Q = THETA(4);
- double V3 = THETA(5);
- double S2 = V2;
- ' '
- $ERROR
- double F = CENT/S2;
- double Y = F*(1+EPS(1)) + EPS(2);
- double IPRED = F;
- ' 'Now, we have the updated model in a format suitable for transport and that reflects all the changes we’ve made to the model object so far.
To read the model back into R, we call mread_yaml()
mod <- mread_yaml("1005-updated.mod", model = "updated")Building updated_mod ... done.mod
--------------- source: updated.mod ---------------
project: /private/var/fol.../T/RtmpUPuIsx
shared object: updated.mod-so-181d45aa25350
time: start: 0 end: 168 delta: 1
add: <none>
compartments: GUT CENT PERIPH [3]
parameters: SEX WT THETA1 THETA2 THETA3 THETA4
THETA5 THETA6 THETA7 [9]
captures: CL Q V2 V3 KA ETA_1 ETA_2 ETA_3 IPRED [9]
omega: 3x3
sigma: 2x2
solver: atol: 1e-08 rtol: 1e-05 maxsteps: 20k
------------------------------------------------------and we have the model object back in hand, ready to simulate. Notice that we’ve read the model back with no knowledge of the estimated NONMEM model. All of those estimates are now hard-coded into the model file. Of course, if the model gets re-estimated and the estimates change, this derived file will be out of date with that run and you’ll have to repeat this process to get updated estimated into the model.
Note that you can also write the model out directly to native mrgsolve format.
mwrite_cpp(mod, file = "1005-updated.txt")In this case, load the model back with mread() as you would any other model
mod <- mread("1005-updated.txt")Building 1005-updated_txt ... done.There is also a translator form the yaml format to the native mrgsolve format. Call yaml_to_cpp() to do this conversion
x <- yaml_to_cpp("1005-updated.mod", model = "1005-converted")
print(basename(x))[1] "1005-converted.mod"2 Addition to evtools plugin
The 1.5.1 mrgsolve release also adds a new function to the evtools plugin that allows you to replace the amount in a specific compartment with a new amount. The function is evt::replace() and it works just like evt::bolus() but the compartment is zeroed-out before putting the new “dose” into the compartment.
The replace functionality uses EVID = 8 that mrgsolve already knows how to handle in event objects and data sets. The addition to the evtools plugin allows you to easily access this functionality from within your model. To zero-out the first compartment and reset the second compartment to 100 at 24 hours into the simulation, the code might look like this
$PLUGIN evtools
$PK
if(TIME==24) {
evt::replace(self, 0, 1);
evt::replace(self, 100, 2);
}3 Refactored as_data_set()
This release also expands the flexibility of as_data_set(), which can combine event objects or data frames into a single data set for simulation.
For example, if we have two event objects, we can combine them into a single object which is a data frame
e1 <- ev(amt = 100)
e2 <- ev(amt = 200)
as_data_set(e1, e2) ID time amt cmt evid
1 1 0 100 1 1
2 2 0 200 1 1We can also make this data set with upper case names using the evd() constructor
e3 <- evd(amt = 100)
e4 <- evd(amt = 200)
as_data_set(e3, e4) ID TIME AMT CMT EVID
1 1 0 100 1 1
2 2 0 200 1 1Or by calling uctran() on the result
as_data_set(e1, e2) %>% uctran() ID TIME AMT CMT EVID
1 1 0 100 1 1
2 2 0 200 1 1The 1.5.1 release gives you more flexibility in combining objects. If you have objects with lower and upper case names, mrgsolve will follow the pattern in the first object
as_data_set(e1, e3, e2, e4) ID time amt cmt evid
1 1 0 100 1 1
2 2 0 100 1 1
3 3 0 200 1 1
4 4 0 200 1 1The 1.5.1 release also lets you combine data frames or a mixture of data frame and event objects
d1 <- mutate(e1, amt = 600) %>% ev_rep(1:2)
d3 <- mutate(e3, amt = 50) %>% ev_rep(1:4)
as_data_set(d3, e3, e4, d1) ID TIME AMT CMT EVID
1 1 0 50 1 1
2 2 0 50 1 1
3 3 0 50 1 1
4 4 0 50 1 1
5 5 0 100 1 1
6 6 0 200 1 1
7 7 0 600 1 1
8 8 0 600 1 14 nm-vars convenience functions
The nm-vars plugin gives you syntax in your mrgsolve model that is similar to what you can use in your NONMEM model. Using this plugin makes it much easier to translate between the two formats.
New in 1.5.1, we have aliases for additional functions
DEXP()LOG10()COS()SIN()
These are just upper case versions of functions that are already in C++.
5 Error when KA is (about) equal to K10
The $PKMODEL block allows you to write mrgsolve models which have analytical solutions to them (you don’t need to solve ODEs). The math behind these models breaks down when KA is equal to K10 (CL/V). Starting with version 1.5.1, mrgsolve will give an error when this is the case.