mod <- update(mod, rtol = 1e-5)1 Model components
This chapter details the different components of a model in mrgsolve. Each component listed here is maintained within a “model object”. This is an updatable S4 object in R that contains all of the basic information required to properly configure and simulate from the model.
1.1 Parameter list
The parameter list is an updatable set of name-value pairs. Referencing the name of an item in the parameter list will substitute the current value associated with that name. While the name “parameter” may have a certain connotation in the modeling world, in mrgsolve a “parameter” could be any category of numeric data: covariates (e.g. WT, AGE, SEX), flags, other numeric data that we commonly call “parameter” (e.g. CL or VC).
The parameter list is declared in the code block $PARAM. While there may be multiple $PARAM blocks in a model, these are condensed to a single parameter list stored in the model object. The names and numbers of all parameters in the model must be declared at the time that the model is compiled. Also, a default value for each parameter must be declared at model compile time, but the value of each parameter may be updated in one of several ways.
The parameters in a model object can be queried or updated with the param() function.
See also: Section 2.2.4, ?param in the R help system after loading mrgsolve.
1.1.1 Central role of parameters in planning simulations
The data items in the parameter list are more than just values associated with a name. When an name is added to the parameter list, that name becomes a key word that mrgsolve will start to recognize in input data sets or when manipulating the model object.
For example, when you want to include a covariate in the model, say weight (WT), you’ll include a column in the data set called WT that will indicate the weight of this or that patient. It is crucial that you also list WT in $PARAM with some default value. It helps if that value is sensible too. When mrgsolve receives the data set prior to simulating, the WT column is matched up with the WT parameter name. As mrgsolve works its way through the input data set (from person to person or from time to time), the value of WT is updated so that the symbol WT in $MAIN or $ODE or $TABLE always points to the value of WT. If the WT name is not in the parameter list, it won’t matter if it is in the data set or not. Only listing a name in $PARAM gets it “into the game”.
Understanding the parameter update mechanism is very important for planning complicated simulations with mrgsolve. Please see the information in Section 4.1 and in Section 12.3.
1.2 Compartment list
Like the parameter list, the compartment list is a series of name-value pairs. The compartment list defines the number, names, and initial values of each compartment in the model. The names, numbers, and order of the compartment in a model is established at the time of model compile and changes to the compartment list require re-compilation of the model.
Compartments are declared in one of two code blocks: $INIT and $CMT. Nominal initial values must be supplied for each compartment. The main difference between $INIT and $CMT is that $CMT assumes a default initial value of 0 for each compartment; thus only compartment names are entered. When using $INIT, both names and values must be explicitly stated for each compartment.
The initial values for each compartment can be queried with the init() function. There are several different ways to set the initial conditions in a model; Section 12.2 illustrates several of these.
See also: Section 12.2 and ?init in the R help system after loading mrgsolve.
1.3 Simulation time grid
The mrgsolve model object stores the parameters for the series of time points to be output for a a simulation. This is the default output time grid that will be used if not over-ridden by another mechanism.
The elements of the simulation time grid are: start, end, delta and add. start, end, delta are passed to seq() as from, to, and by, respectively. add is any arbitrary vector of additional times to simulate.
The simulation time grid in a model object may be queried with the stime() function or by printing the model object to the R console.
See also Section 4.2 for discussion of the simulation time grid and input data sets and Section 1.3.1 and Section 12.4 for using time grid objects .
1.3.1 tgrid objects
A tgrid object has start, end, delta and add attributes. This object is independent of the model object. tgrid objects may be created and combined to create complex sampling designs.
See Section 12.4 for examples and usage.
1.4 Solver settings
mrgsolve uses the DLSODA solver like the one from ODEPACK. Several of the settings for that solver are stored in the model object and passed to the solver when the problem is started. Settings include: atol, rtol, maxsteps, hmax, hmin, ixpr, mxhnil.
Solver settings can be changed through the update() method for your mrgsolve model object. For example, to change rtol to 1e-5, write
where mod is your mrgsolve model object.
1.4.1 atol
Absolute tolerance parameter. Adjust this value lower when you see state variables (compartments) that are becoming very small and possibly turning negative. For example:
mod <- modlib("viral1", end = 144)
out <- mrgsim_e(mod, ev(amt = 1000)) %>% filter(V < 0)
out. # A tibble: 12 × 8
. ID time expos T I V logV logChange
. <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
. 1 1 95 1000 5187630. -1.79e-11 -1.16e-13 NaN NaN
. 2 1 97.5 1000 5298189. -1.92e-10 -1.26e-12 NaN NaN
. 3 1 100 1000 5407830. -2.44e-10 -1.60e-12 NaN NaN
. 4 1 102. 1000 5516561. -2.75e-10 -1.80e-12 NaN NaN
. 5 1 105 1000 5624390. -1.78e-10 -1.17e-12 NaN NaN
. 6 1 108. 1000 5731324. -8.56e-11 -5.60e-13 NaN NaN
. 7 1 110 1000 5837370. -1.24e-11 -8.12e-14 NaN NaN
. 8 1 122. 1000 6354547. -1.67e-11 -1.10e-13 NaN NaN
. 9 1 125 1000 6455421. -1.91e-11 -1.25e-13 NaN NaN
. 10 1 128. 1000 6555459. -2.11e-11 -1.38e-13 NaN NaN
. 11 1 130 1000 6654666. -1.23e-11 -8.07e-14 NaN NaN
. 12 1 132. 1000 6753050. -2.75e-12 -1.82e-14 NaN NaNAdjusting atol to 1E-20 or 1E-30 will prevent this.
mrgsim_e(mod, ev(amt = 1000), atol = 1E-20) %>% filter(time %in% out$time). # A tibble: 12 × 8
. ID time expos T I V logV logChange
. <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
. 1 1 95 1000 5187630. 1.33e-11 8.72e-14 -13.1 -18.7
. 2 1 97.5 1000 5298189. 5.13e-12 3.36e-14 -13.5 -19.2
. 3 1 100 1000 5407830. 1.98e-12 1.30e-14 -13.9 -19.6
. 4 1 102. 1000 5516561. 7.66e-13 5.01e-15 -14.3 -20.0
. 5 1 105 1000 5624390. 2.96e-13 1.94e-15 -14.7 -20.4
. 6 1 108. 1000 5731324. 1.15e-13 7.51e-16 -15.1 -20.8
. 7 1 110 1000 5837370. 4.44e-14 2.91e-16 -15.5 -21.2
. 8 1 122. 1000 6354547. 3.94e-16 2.58e-18 -17.6 -23.3
. 9 1 125 1000 6455422. 1.54e-16 1.01e-18 -18.0 -23.7
. 10 1 128. 1000 6555459. 5.99e-17 3.92e-19 -18.4 -24.1
. 11 1 130 1000 6654666. 2.34e-17 1.53e-19 -18.8 -24.5
. 12 1 132. 1000 6753050. 9.14e-18 5.98e-20 -19.2 -24.91.4.2 rtol
Relative tolerance parameter. Adjust this value lower when you want more precision around the calculation of state variables as the system advances.
1.4.3 maxsteps
This is the maximum number of steps the solver will take when advancing from one time to the next. If the solver can’t make it in maxsteps it will stop and give an error message like this:
DLSODA- At current T (=R1), MXSTEP (=I1) steps
taken on this call before reaching TOUT
In above message, I =
[1] 2000
In above message, R =
[1] 0.0004049985
DLSODA- ISTATE (=I1) illegal.
In above message, I =
[1] -1
DLSODA- Run aborted.. apparent infinite loop.
Error in (function (x, data, idata = null_idata, carry.out = character(0), :
error from XERRWDYou might see this when you have to integrate along time between records in a data set. There isn’t necessarily a problem, but the solver might have to advance over many doses to get to the next record and it only has a limited number of steps it can take between those records before it stops with this error.
When you see this, increase maxsteps to 50000 or larger.
But keep in mind that sometimes the solver can’t make it to the next record because there are issues with the model. It might take thousands of steps to make it 24 hours down the road. In that case, go back to the model code and look for problems in how it is coded.
1.4.4 hmax
The maximum step size. By default, the solver will take steps of different sizes based on what is happening in the simulation. Setting hmax tells the solver not to take a step larger than that value. So in a model where time is in hours, reducing hmax to 0.1 will prevent the solver from taking a step larger than 0.1 hours as it tries to advance to the next time. The will slow down the simulation a bit. But sometimes helpful when the solver starts taking large steps. We don’t recommend using this routinely; for most applications, it should be reserved for troubleshooting situations. If your model doesn’t give the results that you want without setting hmax, we’d recommend a new setup where this isn’t needed.
1.4.5 hmin
The minimum step size. Only set this if you know what you’re doing.
1.4.6 ixpr
A flag to enable printing messages to the R console when the solver switches between non-stiff and stiff solving modes. Rarely used.
1.4.7 mxhnil
The maximum number of messages printed when the model is solving. If you have a lot of messages, keep working on your model code.
1.5 Functions
There are four C++ functions that mrgsolve creates and manages: PREAMBLE, MAIN, ODE, TABLE. Each function is created from an entire code block in the model specification file. The user is responsible for writing correct C++ code in each of these blocks. mrgsolve will parse these blocks and augment this code with the necessary elements to create the C++ function.
These functions may be specified in any order in the model specification file, but there is a specific calling order for these functions. Recognizing and understanding this calling order will help understand how the different pieces of the model specification fit together.
Just prior to starting the problem, mrgsolve calls $PREAMBLE. Then, during advance from time T1 to T2, first $MAIN is called, then $ODE is called repeatedly as the solver finds the values of state variables at T2, and, once the solution is found, $TABLE is called to calculate derived quantities at T2 and to specify variables that should be included in the model output. So, it is helpful to write model specification files in the order:
$PREAMBLE- called only once just prior to processing the first record of the data set$MAIN- before advancing the system$ODE- the system advances toT2$TABLE- after advancing the system
But the order in which they are coded will not affect model compilation or the simulation result.
1.5.1 The $PREAMBLE function
The PREAMBLE function gets called only once, just prior to processing the first record of the data set. This function is composed of C++ code and is used to initialize variables and get them set up prior to starting on the problem.
See Section 2.2.15 for details.
1.5.2 The $MAIN function
The MAIN function gets called at least once before the the solver advances from the current time (T1) to the next time (T2). In the MAIN function, the user may:
- Set initial conditions for any compartment
- Derive new variables to be used in the model
- Write covariate models
- Add between-subject variability to quantities to structural model parameters (e.g.
CLorVC).
In addition to getting called once per record, the MAIN function may be called several times prior to starting the simulation run. The MAIN function is also called whenever the user queries the compartment list.
mrgsolve allows you access compartment amounts in the MAIN function. But it is important to remember the calling order of these model functions (Section 1.5): because MAIN is called before the system advances, the compartment amounts inside this function will reflect the pre-advance values. This is in contrast to accessing compartment amounts inside the TABLE function, which will reflect the values after the system advances. While there are some use cases where it is useful to check the pre-advance compartment amounts in MAIN, most applications should interact with compartment amounts after the system advances in TABLE.
The MAIN block contains code that is equivalent to the code you’d write in NONMEM’s $PK block. You can use either $MAIN or $PK in your mrgsolve model.
See Section 2.2.8 and Section 2.2.9 for details.
1.5.3 The $ODE function
The ODE function is where the user writes the model differential equations. Any derived quantity that depends on a state variable and is used to advanced the system must be calculated inside $ODE. But, this function is called repeatedly during the simulation run, so any calculation that can be moved out of $ODE (for example: to $MAIN) should be.
The ODE block contains code that is equivalent to the code you’d write in NONMEM’s $DES block. You can use either $ODE or $DES in your mrgsolve model.
See Section 2.2.10 and Section 2.2.11 for details.
1.5.4 The $TABLE function
The TABLE function is called after the solver advances in time. The purpose of TABLE is to allow the user to interact with the values of the state variables after advancing, potentially derive new variables, and to insert different outputs into the table of simulated results.
The TABLE block contains code that is equivalent to the code you’d write in NONMEM’s $ERROR block. You can use either $TABLE or $ERROR in your mrgsolve model.
See Section 2.2.12 and Section 2.2.13 for details.
1.6 Random effect variances
The mrgsolve model object keeps track of a arbitrary number of block matrices that are used to simulate variates from multivariate normal distributions. Users can specify OMEGA matrices for simulating between-subject random effects (one draw per individual) or SIGMA matrices for simulating within-subject random effects (one draw per observation).
See Chapter 6 for complete details on how to work with OMEGA and SIGMA in your mrgsolve model.
1.6.1 OMEGA
The matrices are specified in $OMEGA blocks in the model specification file.
OMEGA may be queried or updated with the omat() function.
1.6.2 SIGMA
The matrices are specified in $SIGMA blocks in the model specification file.
SIGMA may be queried or updated by the smat() function.