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
CL

Here 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 1.321 0.0000     0  0     0 0
2:   1  0.0 100.00  0.00 1.321 0.0000     0  0     0 0
3:   1  0.5  85.19 14.57 1.321 0.7019     0  0     0 0
4:   1  1.0  72.57 26.53 1.321 1.2779     0  0     0 0
5:   1  1.5  61.82 36.27 1.321 1.7473     0  0     0 0
6:   1  2.0  52.66 44.15 1.321 2.1265     0  0     0 0
7:   1  2.5  44.86 50.44 1.321 2.4296     0  0     0 0
8:   1  3.0  38.22 55.40 1.321 2.6684     0  0     0 0

This is temporary until the model is re-compiled again

mread("popexample.mod") %>% outvars

Building popexample_mod ... done.

$cmt
[1] "GUT"  "CENT"

$capture
[1] "CL" "DV"

3 Smarter warnings when non-numeric columns are dropped

Pay Attention

This behavior has been updated in more recent mrgsolve versions. mrgsolve will now generate an error when a parameter is in the data set but it is not numeric.

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)

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

try(mrgsim(mod, data))

Error in valid_data_set(data, x, verbose = verbose) : 
  Found input data that cannot be used for simulation
✖ data set column: WT (character)
✖ data set column: TVV (character)

Previously, 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.0000
2:   1  0.0 0.0000
3:   1  0.5 0.7273
4:   1  1.0 1.2449
5:   1  1.5 1.6108
6:   1  2.0 1.8672
7:   1  2.5 2.0445
8:   1  3.0 2.1645

Remember, 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 2.190 KYLE minnesota
4:   1  1.0 3.230 KYLE minnesota
5:   1  1.5 3.700 KYLE minnesota
6:   1  2.0 3.886 KYLE minnesota
7:   1  2.5 3.933 KYLE minnesota
8:   1  3.0 3.912 KYLE minnesota

4 `$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.704

5 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-180fc23fad9f8.so

mrgsolve:::funset(mod)

$symbols
                 name loaded
1   _model_pk1_main__  FALSE
2    _model_pk1_ode__  FALSE
3  _model_pk1_table__  FALSE
4  _model_pk1_event__  FALSE
5 _model_pk1_config__  FALSE

$shlib
               package version compiled
1 pk1-so-180fc23fad9f8   1.5.2     TRUE

We 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  0

It’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    2

Note 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   -1

This 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.

1 Introduction

2 Specify capture at compile time

3 Smarter warnings when non-numeric columns are dropped

4 $ERROR

5 Try loading the model at runtime if it isn’t loaded

6 Time after dose in specific compartment

7 Improved handling of time-varying covariates

4 `$ERROR`