This takes the uncertainty in the model parameter estimates and to simulate a number of theoretical studies. Each study simulates a realization of the parameters from the uncertainty in the fixed parameter estimates. In addition the omega and sigma matrices are simulated from the uncertainty in the Omega/Sigma matrices based on the number of subjects and observations the model was based on.
nlmixrSim(object, ...)
# S3 method for nlmixrFitData
rxSolve(
object,
params = NULL,
events = NULL,
inits = NULL,
scale = NULL,
method = c("liblsoda", "lsoda", "dop853", "indLin"),
transitAbs = NULL,
atol = 1e-08,
rtol = 1e-06,
maxsteps = 70000L,
hmin = 0,
hmax = NA_real_,
hmaxSd = 0,
hini = 0,
maxordn = 12L,
maxords = 5L,
...,
cores,
covsInterpolation = c("locf", "linear", "nocb", "midpoint"),
addCov = FALSE,
matrix = FALSE,
sigma = NULL,
sigmaDf = NULL,
sigmaLower = -Inf,
sigmaUpper = Inf,
nCoresRV = 1L,
sigmaIsChol = FALSE,
sigmaSeparation = c("auto", "lkj", "separation"),
sigmaXform = c("identity", "variance", "log", "nlmixrSqrt", "nlmixrLog",
"nlmixrIdentity"),
nDisplayProgress = 10000L,
amountUnits = NA_character_,
timeUnits = "hours",
stiff,
theta = NULL,
thetaLower = -Inf,
thetaUpper = Inf,
eta = NULL,
addDosing = FALSE,
stateTrim = Inf,
updateObject = FALSE,
omega = NULL,
omegaDf = NULL,
omegaIsChol = FALSE,
omegaSeparation = c("auto", "lkj", "separation"),
omegaXform = c("variance", "identity", "log", "nlmixrSqrt", "nlmixrLog",
"nlmixrIdentity"),
omegaLower = -Inf,
omegaUpper = Inf,
nSub = 1L,
thetaMat = NULL,
thetaDf = NULL,
thetaIsChol = FALSE,
nStud = 1L,
dfSub = 0,
dfObs = 0,
returnType = c("rxSolve", "matrix", "data.frame", "data.frame.TBS", "data.table",
"tbl", "tibble"),
seed = NULL,
nsim = NULL,
minSS = 10L,
maxSS = 1000L,
infSSstep = 12,
strictSS = TRUE,
istateReset = TRUE,
subsetNonmem = TRUE,
maxAtolRtolFactor = 0.1,
from = NULL,
to = NULL,
by = NULL,
length.out = NULL,
iCov = NULL,
keep = NULL,
indLinPhiTol = 1e-07,
indLinPhiM = 0L,
indLinMatExpType = c("expokit", "Al-Mohy", "arma"),
indLinMatExpOrder = 6L,
drop = NULL,
idFactor = TRUE,
mxhnil = 0,
hmxi = 0,
warnIdSort = TRUE,
warnDrop = TRUE,
ssAtol = 1e-08,
ssRtol = 1e-06,
safeZero = TRUE,
sumType = c("pairwise", "fsum", "kahan", "neumaier", "c"),
prodType = c("long double", "double", "logify"),
sensType = c("advan", "autodiff", "forward", "central"),
linDiff = c(tlag = 1.5e-05, f = 1.5e-05, rate = 1.5e-05, dur = 1.5e-05, tlag2 =
1.5e-05, f2 = 1.5e-05, rate2 = 1.5e-05, dur2 = 1.5e-05),
linDiffCentral = c(tlag = TRUE, f = TRUE, rate = TRUE, dur = TRUE, tlag2 = TRUE, f2 =
TRUE, rate2 = TRUE, dur2 = TRUE),
resample = NULL,
resampleID = TRUE
)
# S3 method for nlmixrFitData
simulate(object, nsim = 1, seed = NULL, ...)
# S3 method for nlmixrFitData
solve(a, b, ...)nlmixr object
Other arguments sent to rxSolve
a numeric named vector with values for every parameter in the ODE system; the names must correspond to the parameter identifiers used in the ODE specification;
an eventTable object describing the input
(e.g., doses) to the dynamic system and observation sampling
time points (see eventTable());
a vector of initial values of the state variables (e.g., amounts in each compartment), and the order in this vector must be the same as the state variables (e.g., PK/PD compartments);
a numeric named vector with scaling for ode
parameters of the system. The names must correspond to the
parameter identifiers in the ODE specification. Each of the
ODE variables will be divided by the scaling factor. For
example scale=c(center=2) will divide the center ODE
variable by 2.
The method for solving ODEs. Currently this supports:
"liblsoda" thread safe lsoda. This supports parallel
thread-based solving, and ignores user Jacobian specification.
"lsoda" -- LSODA solver. Does not support parallel thread-based
solving, but allows user Jacobian specification.
"dop853" -- DOP853 solver. Does not support parallel thread-based
solving nor user Jacobain specification
"indLin" -- Solving through inductive linearization. The RxODE dll
must be setup specially to use this solving routine.
boolean indicating if this is a transit compartment absorption
a numeric absolute tolerance (1e-8 by default) used by the ODE solver to determine if a good solution has been achieved; This is also used in the solved linear model to check if prior doses do not add anything to the solution.
a numeric relative tolerance (1e-6 by default) used
by the ODE solver to determine if a good solution has been
achieved. This is also used in the solved linear model to check
if prior doses do not add anything to the solution.
maximum number of (internally defined) steps allowed during one call to the solver. (5000 by default)
The minimum absolute step size allowed. The default value is 0.
The maximum absolute step size allowed. When
hmax=NA (default), uses the average difference +
hmaxSd*sd in times and sampling events. The hmaxSd is a user
specified parameter and which defaults to zero. When
hmax=NULL RxODE uses the maximum difference in times in
your sampling and events. The value 0 is equivalent to infinite
maximum absolute step size.
The number of standard deviations of the time difference to add to hmax. The default is 0
The step size to be attempted on the first step. The
default value is determined by the solver (when hini = 0)
The maximum order to be allowed for the nonstiff (Adams) method. The default is 12. It can be between 1 and 12.
The maximum order to be allowed for the stiff (BDF) method. The default value is 5. This can be between 1 and 5.
Number of cores used in parallel ODE solving. This
is equivalent to calling setRxThreads()
specifies the interpolation method for
time-varying covariates. When solving ODEs it often samples
times outside the sampling time specified in events.
When this happens, the time varying covariates are
interpolated. Currently this can be:
"linear" interpolation, which interpolates the covariate
by solving the line between the observed covariates and extrapolating the new
covariate value.
"constant" -- Last observation carried forward (the default).
"NOCB" -- Next Observation Carried Backward. This is the same method
that NONMEM uses.
"midpoint" Last observation carried forward to midpoint; Next observation
carried backward to midpoint.
A boolean indicating if covariates should be added to the output matrix or data frame. By default this is disabled.
A boolean indicating if a matrix should be returned instead of the RxODE's solved object.
Named sigma covariance or Cholesky decomposition of a covariance matrix. The names of the columns indicate parameters that are simulated. These are simulated for every observation in the solved system.
Degrees of freedom of the sigma t-distribution. By
default it is equivalent to Inf, or a normal distribution.
Lower bounds for simulated unexplained variability (by default -Inf)
Upper bounds for simulated unexplained variability (by default Inf)
Number of cores used for the simulation of the sigma variables. By default this is 1. To reproduce the results you need to run on the same platform with the same number of cores. This is the reason this is set to be one, regardless of what the number of cores are used in threaded ODE solving.
Boolean indicating if the sigma is in the Cholesky decomposition instead of a symmetric covariance
separation strategy for sigma;
Tells the type of separation strategy when
simulating covariance with parameter uncertainty with standard
deviations modeled in the thetaMat matrix.
"lkj" simulates the correlation matrix from the
rLKJ1 matrix with the distribution parameter eta
equal to the degrees of freedom nu by (nu-1)/2
"separation" simulates from the identity inverse Wishart
covariance matrix with nu degrees of freedom. This is then
converted to a covariance matrix and augmented with the modeled
standard deviations. While computationally more complex than the
"lkj" prior, it performs better when the covariance matrix
size is greater or equal to 10
"auto" chooses "lkj" when the dimension of the
matrix is less than 10 and "separation" when greater
than equal to 10.
When taking sigma values from the thetaMat
simulations (using the separation strategy for covariance
simulation), how should the thetaMat values be turned int
standard deviation values:
identity This is when standard deviation values are
directly modeled by the params and thetaMat matrix
variance This is when the params and thetaMat
simulates the variance that are directly modeled by the
thetaMat matrix
log This is when the params and thetaMat
simulates log(sd)
nlmixrSqrt This is when the params and
thetaMat simulates the inverse cholesky decomposed matrix
with the x^2 modeled along the diagonal. This only works
with a diagonal matrix.
nlmixrLog This is when the params and
thetaMat simulates the inverse cholesky decomposed matrix
with the exp(x^2) along the diagonal. This only works
with a diagonal matrix.
nlmixrIdentity This is when the params and
thetaMat simulates the inverse cholesky decomposed matrix.
This only works with a diagonal matrix.
An integer indicating the minimum number of c-based solves before a progress bar is shown. By default this is 10,000.
This supplies the dose units of a data frame supplied instead of an event table. This is for importing the data as an RxODE event table.
This supplies the time units of a data frame supplied instead of an event table. This is for importing the data as an RxODE event table.
a logical (TRUE by default) indicating whether
the ODE system is stiff or not.
For stiff ODE systems (`stiff = TRUE`), `RxODE` uses the
LSODA (Livermore Solver for Ordinary Differential Equations)
Fortran package, which implements an automatic method switching
for stiff and non-stiff problems along the integration
interval, authored by Hindmarsh and Petzold (2003).
For non-stiff systems (`stiff = FALSE`), `RxODE` uses
DOP853, an explicit Runge-Kutta method of order 8(5, 3) of
Dormand and Prince as implemented in C by Hairer and Wanner
(1993).
If stiff is not specified, the `method` argument is used instead.
A vector of parameters that will be named THETA\[#\] and
added to parameters
Lower bounds for simulated population parameter
variability (by default -Inf)
Upper bounds for simulated population unexplained
variability (by default Inf)
A vector of parameters that will be named ETA\[#\] and
added to parameters
Boolean indicating if the solve should add RxODE
EVID and related columns. This will also include dosing
information and estimates at the doses. Be default, RxODE
only includes estimates at the observations. (default
FALSE). When addDosing is NULL, only
include EVID=0 on solve and exclude any model-times or
EVID=2. If addDosing is NA the classic
RxODE EVID events are returned. When addDosing is TRUE
add the event information in NONMEM-style format; If
subsetNonmem=FALSE RxODE will also include extra event types
(EVID) for ending infusion and modeled times:
EVID=-1 when the modeled rate infusions are turned
off (matches rate=-1)
EVID=-2 When the modeled duration infusions are
turned off (matches rate=-2)
EVID=-10 When the specified rate infusions are
turned off (matches rate>0)
EVID=-20 When the specified dur infusions are
turned off (matches dur>0)
EVID=101,102,103,... Modeled time where 101 is the
first model time, 102 is the second etc.
When amounts/concentrations in one of the states
are above this value, trim them to be this value. By default
Inf. Also trims to -stateTrim for large negative
amounts/concentrations. If you want to trim between a range
say c(0, 2000000) you may specify 2 values with a lower and
upper range to make sure all state values are in the
reasonable range.
This is an internally used flag to update the
RxODE solved object (when supplying an RxODE solved object) as
well as returning a new object. You probably should not
modify it's FALSE default unless you are willing to
have unexpected results.
Estimate of Covariance matrix. When omega is a list, assume it is a block matrix and convert it to a full matrix for simulations.
The degrees of freedom of a t-distribution for
simulation. By default this is NULL which is
equivalent to Inf degrees, or to simulate from a normal
distribution instead of a t-distribution.
Indicates if the omega supplied is a
Cholesky decomposed matrix instead of the traditional
symmetric matrix.
Omega separation strategy
Tells the type of separation strategy when
simulating covariance with parameter uncertainty with standard
deviations modeled in the thetaMat matrix.
"lkj" simulates the correlation matrix from the
rLKJ1 matrix with the distribution parameter eta
equal to the degrees of freedom nu by (nu-1)/2
"separation" simulates from the identity inverse Wishart
covariance matrix with nu degrees of freedom. This is then
converted to a covariance matrix and augmented with the modeled
standard deviations. While computationally more complex than the
"lkj" prior, it performs better when the covariance matrix
size is greater or equal to 10
"auto" chooses "lkj" when the dimension of the
matrix is less than 10 and "separation" when greater
than equal to 10.
When taking omega values from the thetaMat
simulations (using the separation strategy for covariance
simulation), how should the thetaMat values be turned int
standard deviation values:
identity This is when standard deviation values are
directly modeled by the params and thetaMat matrix
variance This is when the params and thetaMat
simulates the variance that are directly modeled by the
thetaMat matrix
log This is when the params and thetaMat
simulates log(sd)
nlmixrSqrt This is when the params and
thetaMat simulates the inverse cholesky decomposed matrix
with the x^2 modeled along the diagonal. This only works
with a diagonal matrix.
nlmixrLog This is when the params and
thetaMat simulates the inverse cholesky decomposed matrix
with the exp(x^2) along the diagonal. This only works
with a diagonal matrix.
nlmixrIdentity This is when the params and
thetaMat simulates the inverse cholesky decomposed matrix.
This only works with a diagonal matrix.
Lower bounds for simulated ETAs (by default -Inf)
Upper bounds for simulated ETAs (by default Inf)
Number between subject variabilities (ETAs) simulated for every
realization of the parameters.
Named theta matrix.
The degrees of freedom of a t-distribution for
simulation. By default this is NULL which is
equivalent to Inf degrees, or to simulate from a normal
distribution instead of a t-distribution.
Indicates if the theta supplied is a
Cholesky decomposed matrix instead of the traditional
symmetric matrix.
Number virtual studies to characterize uncertainty in estimated parameters.
Degrees of freedom to sample the between subject variability matrix from the inverse Wishart distribution (scaled) or scaled inverse chi squared distribution.
Degrees of freedom to sample the unexplained variability matrix from the inverse Wishart distribution (scaled) or scaled inverse chi squared distribution.
This tells what type of object is returned. The currently supported types are:
"rxSolve" (default) will return a reactive data frame
that can change easily change different pieces of the solve and
update the data frame. This is the currently standard solving
method in RxODE, is used for rxSolve(object, ...), solve(object,...),
"data.frame" -- returns a plain, non-reactive data
frame; Currently very slightly faster than returnType="matrix"
"matrix" -- returns a plain matrix with column names attached
to the solved object. This is what is used object$run as well as object$solve
"data.table" -- returns a data.table; The data.table is
created by reference (ie setDt()), which should be fast.
"tbl" or "tibble" returns a tibble format.
an object specifying if and how the random number generator should be initialized
represents the number of simulations. For RxODE, if
you supply single subject event tables (created with
[eventTable()])
Minimum number of iterations for a steady-state dose
Maximum number of iterations for a steady-state dose
Step size for determining if a constant infusion has reached steady state. By default this is large value, 420.
Boolean indicating if a strict steady-state is
required. If a strict steady-state is (TRUE) required
then at least minSS doses are administered and the
total number of steady states doses will continue until
maxSS is reached, or atol and rtol for
every compartment have been reached. However, if ODE solving
problems occur after the minSS has been reached the
whole subject is considered an invalid solve. If
strictSS is FALSE then as long as minSS
has been reached the last good solve before ODE solving
problems occur is considered the steady state, even though
either atol, rtol or maxSS have not
been achieved.
When TRUE, reset the ISTATE variable to 1 for
lsoda and liblsoda with doses, like deSolve; When FALSE, do
not reset the ISTATE variable with doses.
subset to NONMEM compatible EVIDs only. By
default TRUE.
The maximum atol/rtol that
FOCEi and other routines may adjust to. By default 0.1
When there is no observations in the event table, start observations at this value. By default this is zero.
When there is no observations in the event table, end observations at this value. By default this is 24 + maximum dose time.
When there are no observations in the event table, this
is the amount to increment for the observations between from
and to.
The number of observations to create if there isn't any observations in the event table. By default this is 200.
A data frame of individual non-time varying covariates
to combine with the events dataset by merge.
Columns to keep from either the input dataset or the
iCov dataset. With the iCov dataset, the column
is kept once per line. For the input dataset, if any records
are added to the data LOCF (Last Observation Carried forward)
imputation is performed.
the requested accuracy tolerance on exponential matrix.
the maximum size for the Krylov basis
This is them matrix exponential type that is use for RxODE. Currently the following are supported:
Al-Mohy Uses the exponential matrix method of Al-Mohy Higham (2009)
arma Use the exponential matrix from RcppArmadillo
expokit Use the exponential matrix from Roger B. Sidje (1998)
an integer, the order of approximation to
be used, for the Al-Mohy and expokit values.
The best value for this depends on machine precision (and
slightly on the matrix). We use 6 as a default.
Columns to drop from the output
This boolean indicates if original ID values should be maintained. This changes the default sequentially ordered ID to a factor with the original ID values in the original dataset. By default this is enabled.
maximum number of messages printed (per problem)
warning that T + H = T on a step (H = step size). This must
be positive to result in a non-default value. The default
value is 0 (or infinite).
inverse of the maximum absolute value of H to are used.
hmxi = 0.0 is allowed and corresponds to an infinite hmax1 (default). hminandhmximay be changed at any time, but will not take effect until the next change ofHis considered. This option is only considered withmethod="liblsoda"`.
Warn if the ID is not present and RxODE assumes the order of the parameters/iCov are the same as the order of the parameters in the input dataset.
Warn if column(s) were supposed to be dropped, but were not present.
Steady state atol convergence factor. Can be a vector based on each state.
Steady state rtol convergence factor. Can be a vector based on each state.
Use safe zero divide and log routines. By default this is turned on but you may turn it off if you wish.
Sum type to use for sum() in
RxODE code blocks.
pairwise uses the pairwise sum (fast, default)
fsum uses Python's fsum function (most accurate)
kahan uses Kahan correction
neumaier uses Neumaier correction
c uses no correction: default/native summing
Product to use for prod() in RxODE blocks
long double converts to long double, performs the
multiplication and then converts back.
double uses the standard double scale for multiplication.
Sensitivity type for linCmt() model:
advan Use the direct advan solutions
autodiff Use the autodiff advan solutions
forward Use forward difference solutions
central Use central differences
This gives the linear difference amount for all the types of linear compartment model parameters where sensitivities are not calculated. The named components of this numeric vector are:
"lag" Central compartment lag
"f" Central compartment bioavailability
"rate" Central compartment modeled rate
"dur" Central compartment modeled duration
"lag2" Depot compartment lag
"f2" Depot compartment bioavailability
"rate2" Depot compartment modeled rate
"dur2" Depot compartment modeled duration
This gives the which parameters use central
differences for the linear compartment model parameters. The
are the same components as linDiff
A character vector of model variables to resample
from the input dataset; This sampling is done with replacement.
When NULL or FALSE no resampling is done. When
TRUE resampling is done on all covariates in the input
dataset
boolean representing if the resampling should be
done on an individual basis TRUE (ie. a whole patient is
selected) or each covariate is resampled independent of the
subject identifier FALSE. When resampleID=TRUE
correlations of parameters are retained, where as when
resampleID=FALSE ignores patient covariate correaltions.
Hence the default is resampleID=TRUE.
when using solve(), this is equivalent to the
object argument. If you specify object later in
the argument list it overwrites this parameter.
when using solve(), this is equivalent to the
params argument. If you specify params as a
named argument, this overwrites the output
A RxODE solved object