library(tidyverse)
library(mrgsolve)1 Introduction
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.
2 Specify capture at compile time
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 ... (waiting) ...
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 1.057 0.000 0 0 0 0
2: 1 0.0 100.00 0.00 1.057 0.000 0 0 0 0
3: 1 0.5 83.53 16.28 1.057 0.675 0 0 0 0
4: 1 1.0 69.77 29.53 1.057 1.224 0 0 0 0
5: 1 1.5 58.28 40.26 1.057 1.669 0 0 0 0
6: 1 2.0 48.69 48.88 1.057 2.026 0 0 0 0
7: 1 2.5 40.67 55.74 1.057 2.311 0 0 0 0
8: 1 3.0 33.97 61.16 1.057 2.535 0 0 0 0This is temporary until the model is re-compiled again
mread("popexample.mod") %>% outvarsBuilding popexample_mod ... done.$cmt
[1] "GUT" "CENT"
$capture
[1] "CL" "DV"3 Smarter warnings when non-numeric columns are dropped
mrgsolve warns you when non-numeric columns are dropped from the data set when you run the model. For example
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 70 Let’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 column: WT (character)[data-set] dropped column: TVV (character)Model: popex
Dim: 482 x 3
Time: 0 to 240
ID: 1
ID time IPRED
1: 1 0.0 0.0000
2: 1 0.0 0.0000
3: 1 0.5 0.4051
4: 1 1.0 0.7667
5: 1 1.5 1.0887
6: 1 2.0 1.3745
7: 1 2.5 1.6272
8: 1 3.0 1.8498Previously, 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 1.215
4: 1 1.0 2.020
5: 1 1.5 2.547
6: 1 2.0 2.882
7: 1 2.5 3.087
8: 1 3.0 3.204Remember, 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.000 KYLE minnesota
2: 1 0.0 0.000 KYLE minnesota
3: 1 0.5 3.058 KYLE minnesota
4: 1 1.0 4.876 KYLE minnesota
5: 1 1.5 5.923 KYLE minnesota
6: 1 2.0 6.493 KYLE minnesota
7: 1 2.5 6.767 KYLE minnesota
8: 1 3.0 6.859 KYLE minnesota4 $ERROR
There 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.7045 Try loading the model at runtime if it isn’t loaded
Sometimes, 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")Building pk1 ... done.mrgsolve:::unloadso.mrgmod(mod)unloaded pk1-so-59812455e037.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-59812455e037 1.0.8.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.
6 Time after dose in specific compartment
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.
7 Improved handling of time-varying covariates
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.