This chapter details the different components of a model in
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.
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.
SEX), flags, other numeric data that we
commonly call “parameter” (e.g.
The parameter list is declared in the code block
$PARAM. While there may be
$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
See also: 2.2.4,
?param in the
R help system after
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
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
$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
updated so that the symbol
$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”.
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:
initial values must be supplied for each compartment. The main difference
$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
function. There are several different ways to set the initial conditions in a
model; section 11.2 illustrates several of these.
See also: section 11.2 and
?init in the
R help system after
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:
delta are passed to
add is any arbitrary vector of additional times to simulate.
The simulation time grid in a model object may be queried with the
function or by printing the model object to the
tgrid object has
add attributes. This object
is independent of the model object.
tgrid objects may be created and combined
to create complex sampling designs.
See section 11.4 for examples and usage.
mrgsolve uses the
DLSODA solver 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:
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)
. Loading model from cache.
out <- mrgsim_e(mod, ev(amt = 1000)) %>% filter(V < 0) out
. # A tibble: 6 × 8 . ID time expos T I V logV logChange . <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> . 1 1 95 1000 5187630. -2.74e-10 -1.79e-12 NaN NaN . 2 1 97.5 1000 5298189. -3.18e-10 -2.09e-12 NaN NaN . 3 1 120 1000 6252828. -1.55e-10 -1.01e-12 NaN NaN . 4 1 122. 1000 6354547. -2.90e-10 -1.90e-12 NaN NaN . 5 1 125 1000 6455422. -1.97e-10 -1.30e-12 NaN NaN . 6 1 128. 1000 6555459. -7.81e-11 -5.21e-13 NaN NaN
1E-30 will prevent this.
. # A tibble: 6 × 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 120 1000 6252828. 1.01e-15 6.62e-18 -17.2 -22.9 . 4 1 122. 1000 6354547. 3.94e-16 2.58e-18 -17.6 -23.3 . 5 1 125 1000 6455422. 1.54e-16 1.01e-18 -18.0 -23.7 . 6 1 128. 1000 6555459. 5.99e-17 3.92e-19 -18.4 -24.1
Relative tolerance parameter. Adjust this value lower when you want more precision around the calculation of state variables as the system advances.
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:
- At current T (=R1), MXSTEP (=I1) steps DLSODA taken on this call before reaching TOUT , I = In above message 2000 , R = In above message 0.0004049985 - ISTATE (=I1) illegal. DLSODA, I = In above message -1 - Run aborted.. apparent infinite loop. DLSODA(function (x, data, idata = null_idata, carry.out = character(0), : Error in error from XERRWD
You 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.
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
in hours, reducing
0.1 will prevent the solver from taking a step
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.
A flag to enable printing messages to the R console when the solver switches between non-stiff and stiff solving modes. Rarely used.
There are four
C++ functions that
mrgsolve creates and manages:
TABLE. Each function is created from an entire code block in
the model specification file. The user is responsible for writing correct
code in each of these blocks. mrgsolve will parse these blocks and augment this
code with the necessary elements to create the
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
$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 to
$TABLE- after advancing the system
But the order in which they are coded will not affect model compilation or the simulation result.
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
See 2.2.13 for details.
MAIN function gets called at least once before the the solver advances
from the current time (
T1) to the next time (
T2). In the
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
In addition to getting called once per record, the
MAIN function may be
called several times prior to starting the simulation run. The
function is also called whenever the user queries the compartment list.
See 2.2.7 for details.
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.
See 2.2.9 for details.
TABLE function is called after the solver advances in time. The
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.
See 2.2.11 for details.
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).
The user may use the
revar() function to query both
The matrices are specified in
$OMEGA blocks in the model specification file.
OMEGA may be queried or updated with the