This blog post shows some of the new features that came into mrgsolve
starting with version 0.10.7 (December 2020). The purpose
is to illustrate what is possible rather than detailed documentation.
There were other small bugs and gremlins fixed as well, but we are
focused on bigger, user-facing changes here.
Big Idea: The feature helps give you more
flexibility for getting data out of your model and into your simulated
output. You have the ability to select output items when you run
mread() in addition to specifying them in the model
code.
In your mrgsolve model, you can request that the values of different
variables are written into the simulated output. This is done through
the capture block or by declaring variables as type
capture:
mod <- mread("popexample.mod"). Building popexample_mod ... done.blocks(mod, MAIN, TABLE, CAPTURE).
. Model file: popexample.mod
.
. $MAIN
. double CL = exp(log(TVCL) + 0.75*log(WT/70) + ECL);
. double V = exp(log(TVV) + log(WT/70) + EV );
. double KA = exp(log(TVKA) + EKA);
.
. $TABLE
. double IPRED = CENT/V;
. capture DV = IPRED*exp(EPS(1));
.
. $CAPTURE
. CLHere we have derived several variables (CL,
V, KA, IPRED, DV)
and some of them have been marked for “capture” (DV and
CL).
We can check what will be captured:
outvars(mod). $cmt
. [1] "GUT" "CENT"
.
. $capture
. [1] "CL" "DV"Starting with mrgsolve 0.10.7, you can specify additional variables
for capture when you read in the model with mread():
mod <- mread("popexample.mod", capture = "ETA(1), WT, IPRED, V"). Building popexample_mod ... done.Here, I’ve asked for the value of one of the ETAs, a parameter
(WT) and a couple of derived quantities (IPRED
and V):
outvars(mod). $cmt
. [1] "GUT" "CENT"
.
. $capture
. [1] "CL" "DV" "ETA_1" "WT" "IPRED" "V"And now when I simulate from the model, I’ll get all of these quantities back:
mrgsim(mod, ev(amt = 100)). Model: popexample_mod
. Dim: 482 x 10
. Time: 0 to 240
. ID: 1
. ID time GUT CENT CL DV ETA_1 WT IPRED V
. 1: 1 0.0 0.00 0.00 0.3599 0.000 -1.022 70 0.000 18.93
. 2: 1 0.0 100.00 0.00 0.3599 0.000 -1.022 70 0.000 18.93
. 3: 1 0.5 75.40 24.48 0.3599 1.293 -1.022 70 1.293 18.93
. 4: 1 1.0 56.86 42.70 0.3599 2.256 -1.022 70 2.256 18.93
. 5: 1 1.5 42.87 56.21 0.3599 2.970 -1.022 70 2.970 18.93
. 6: 1 2.0 32.33 66.17 0.3599 3.496 -1.022 70 3.496 18.93
. 7: 1 2.5 24.37 73.46 0.3599 3.881 -1.022 70 3.881 18.93
. 8: 1 3.0 18.38 78.73 0.3599 4.159 -1.022 70 4.159 18.93This is temporary until the model is re-compiled again
mread("popexample.mod") %>% outvars. Building popexample_mod ... (waiting) ...
. done.. $cmt
. [1] "GUT" "CENT"
.
. $capture
. [1] "CL" "DV"mrgsolve warns you when non-numeric columns are dropped from the data set when you run the model. For example:
library(mrgsolve)
mod <- modlib("popex") %>% update(outvars = "IPRED")In this model, we have the following parameters:
param(mod).
. Model parameters (N=4):
. name value . name value
. TVCL 1 | TVV 24
. TVKA 0.5 | WT 70Let’s create an input data set with some character values: two parameter columns and one non-parameter column:
data <- expand.ev(amt = 100, WT = "500", TVV = "ABC", name = "KYLE")When we simulation, we get some messages for dropped columns:
mrgsim(mod, data). [data-set] dropped non-numeric: WT. [data-set] dropped non-numeric: TVV. Model: popex
. Dim: 482 x 3
. Time: 0 to 240
. ID: 1
. ID time IPRED
. 1: 1 0.0 0.000
. 2: 1 0.0 0.000
. 3: 1 0.5 2.302
. 4: 1 1.0 3.419
. 5: 1 1.5 3.949
. 6: 1 2.0 4.189
. 7: 1 2.5 4.285
. 8: 1 3.0 4.311Previously, mrgsolve would have warned you about all 3 columns; now it only warns you about columns that are dropped and would have been relevant to the simulation.
This example is made up to illustrate the behavior. But in typical usage when we do things right, we might have extra character columns that we don’t care about. In this case we get no warnings:
data <- expand.ev(
amt = 100, WT = 70, TVV = 20,
name = "KYLE",
state = "minnesota"
)
mrgsim(mod, data). Model: popex
. Dim: 482 x 3
. Time: 0 to 240
. ID: 1
. ID time IPRED
. 1: 1 0.0 0.000
. 2: 1 0.0 0.000
. 3: 1 0.5 2.111
. 4: 1 1.0 3.631
. 5: 1 1.5 4.707
. 6: 1 2.0 5.451
. 7: 1 2.5 5.946
. 8: 1 3.0 6.255Remember, we can still ask for non-numeric items to be brought into
the data set with the recover argument:
mrgsim(mod, data, recover = "name,state"). Model: popex
. Dim: 482 x 5
. Time: 0 to 240
. ID: 1
. ID time IPRED name state
. 1: 1 0.0 0.0000 KYLE minnesota
. 2: 1 0.0 0.0000 KYLE minnesota
. 3: 1 0.5 0.6898 KYLE minnesota
. 4: 1 1.0 1.3054 KYLE minnesota
. 5: 1 1.5 1.8535 KYLE minnesota
. 6: 1 2.0 2.3404 KYLE minnesota
. 7: 1 2.5 2.7715 KYLE minnesota
. 8: 1 3.0 3.1520 KYLE minnesota$ERRORThere is a new block alias called $ERROR that is just an
alias for $TABLE
code <- '
$PARAM CL = 1, V = 20, KA = 1
$PKMODEL cmt = "CENT"
$ERROR
capture CP = CENT/V;
'
mod <- mcode("error", code). Building error ... done.mrgsim(mod, ev(amt = 100)). Model: error
. Dim: 26 x 4
. Time: 0 to 24
. ID: 1
. ID time CENT CP
. 1: 1 0 0.00 0.000
. 2: 1 0 100.00 5.000
. 3: 1 1 95.12 4.756
. 4: 1 2 90.48 4.524
. 5: 1 3 86.07 4.304
. 6: 1 4 81.87 4.094
. 7: 1 5 77.88 3.894
. 8: 1 6 74.08 3.704Sometimes, you find yourself with a model object but the shared object hasn’t been loaded into the session. This could happen when you are running simulations in parallel and the model doesn’t get loaded in the worker R process. For the post, I’ll just phony up a model that happens to not be loaded:
This code loads and then unloads the model, and then shows that it is not loaded:
mod <- modlib("pk1"). Loading model from cache.mrgsolve:::unloadso.mrgmod(mod). unloaded pk1-so-2fe31ccfab94.somrgsolve:::funset(mod). $symbols
. name loaded
. 1 _model_pk1_main__ FALSE
. 2 _model_pk1_ode__ FALSE
. 3 _model_pk1_table__ FALSE
. 4 _model_pk1_config__ FALSE
.
. $shlib
. package version compiled
. 1 pk1-so-2fe31ccfab94 1.0.0.9000 TRUEWe can still simulate from this model; mrgsolve will determine that the model isn’t loaded and it will try to load it for you.
mrgsim(mod). Model: pk1
. Dim: 25 x 5
. Time: 0 to 24
. ID: 1
. ID time EV CENT CP
. 1: 1 0 0 0 0
. 2: 1 1 0 0 0
. 3: 1 2 0 0 0
. 4: 1 3 0 0 0
. 5: 1 4 0 0 0
. 6: 1 5 0 0 0
. 7: 1 6 0 0 0
. 8: 1 7 0 0 0It’s possible that the model can’t be loaded (e.g. due to missing
.so file). In that case, mrgsolve will give you the usual
error message that the model needs to be compiled.
Version 0.10.7 also adds a new plugin with the ability to calculate time after dose in any compartment.
We write a model using the tad plugin to track time
after dose in compartments one and two. We
create tadose objects to track this and we can call the
tad() method on these objects, passing in the
self data item.
code <- '
[plugin] tad
[ global ]
mrg::tadose tad_cmt_1(1);
mrg::tadose tad_cmt_2(2);
[ pkmodel ] cmt = "GUT,CENT", depot = TRUE
[ param ] CL = 1, V = 20, KA = 1
[ main ]
capture tad1 = tad_cmt_1.tad(self);
capture tad2 = tad_cmt_2.tad(self);
'
mod <- mcode("tad", code, soloc = '.'). Building tad ... done.data <- c(
ev(amt = 100, cmt = 1, time = 1),
ev(amt = 200, cmt = 2, time = 3)
)
mrgsim(mod, data). Model: tad
. Dim: 27 x 6
. Time: 0 to 24
. ID: 1
. ID time GUT CENT tad1 tad2
. 1: 1 0 0.000 0.00 -1 -1
. 2: 1 1 0.000 0.00 -1 -1
. 3: 1 1 100.000 0.00 0 -1
. 4: 1 2 36.788 61.41 1 -1
. 5: 1 3 13.534 81.00 2 -1
. 6: 1 3 13.534 281.00 2 0
. 7: 1 4 4.979 275.61 3 1
. 8: 1 5 1.832 265.22 4 2Note that time after dose is -1 until a dose is administered.
Recall also that time after dose can be calculated more simply if
there is only one dose type by passing the tad
argument:
data1 <- filter(data, cmt ==1) %>% mutate(time = 3)
mod %>% mrgsim(data1, tad = TRUE). Model: tad
. Dim: 26 x 7
. Time: 0 to 24
. ID: 1
. ID time tad GUT CENT tad1 tad2
. 1: 1 0 -3 0.000 0.00 -1 -1
. 2: 1 1 -2 0.000 0.00 -1 -1
. 3: 1 2 -1 0.000 0.00 -1 -1
. 4: 1 3 0 0.000 0.00 -1 -1
. 5: 1 3 0 100.000 0.00 0 -1
. 6: 1 4 1 36.788 61.41 1 -1
. 7: 1 5 2 13.534 81.00 2 -1
. 8: 1 6 3 4.979 85.36 3 -1This is a little nicer because it will fill in negative
tad values for you.
This one is a little harder to demonstrate. mrgsolve has always been
able to simulate time-varying covariates, but it hasn’t always inserted
a hard discontinuity when parameter values change in the data set over
time. Starting in 0.10.7, when parameters are read from the
data-set and they change from record to record, we now
re-initialize the ode solver so a hard discontinuity is created in the
simulation. This will create more accurate simulations around the time
that the covariate value changes. This results in no change when running
1- and 2-compartment models which are handled with the closed form
equations.
mrgsolve: mrgsolve.github.io | metrum research group: metrumrg.com