Calculate variations

Calculates the variational matrix of the given model, that is, the first and optionally second order derivatives of the model outputs with respect to the model parameters, at given times. The user can choose whether the computation uses numerical or symbolic derivatives.

calcVariations(
  model,
  p,
  init,
  output,
  times,
  theta,
  eta,
  eps,
  vartheta = names(theta),
  vareta = names(eta),
  vareps = names(eps),
  secOrd = FALSE,
  symbolic = TRUE,
  chkModel = TRUE,
  showWarn = TRUE,
  ...
)

Arguments

model

Function(t, y, p) of time, state variables and parameters, specifying the differential equations. This function should return the numeric vector dy/dt. See the notes for requirements in case of symbolic derivations.

p

Function(theta, eta) of population and individual parameters specifying the model parameters. This function should return a named numeric vector. See the notes for requirements in case of symbolic derivations.

init

Function(p) of parameters specifying the initial state (at time 0) of the model. See the notes for requirements in case of symbolic derivations.

output

Function(y, p, eps) of state variables, parameters and residual errors specifying the model outputs. This function should return a numeric vector. See the notes for requirements in case of symbolic derivations.

times

Numeric vector of times where variations are to be evaluated. Should all be >= 0, but do not need to include the dose time 0.

theta

Population parameter values, as named numeric vector. May be NULL if none are needed.

eta

Individual parameter values, as named numeric vector. May be NULL if none are needed.

eps

Residual parameter values, as named numeric vector. May be NULL if none are needed.

vartheta

Vector of names of population parameters for which variations are to be calculated. Should be a subset of names(theta), or NULL for none. By default equal to names(theta).

vareta

Vector of names of individual parameters for which variations are to be calculated. Should be a subset of names(eta), or NULL for none. By default equal to names(eta).

vareps

Vector of names of residual parameters for which variations are to be calculated. Should be a subset of names(eps), or NULL for none. By default equal to names(eps).

secOrd

TRUE if second order derivatives have to be computed, FALSE (default) if not.

symbolic

TRUE (default) if derivatives are to be computed symbolically, FALSE if numerically.

chkModel

TRUE (default) if it has to be checked whether model components (model, p, init, output) are formatted correctly for symbolic manipulation, FALSE if not.

showWarn

TRUE (default) if warnings about missing random parameters are to be shown, FALSE if not. The latter is of use in case only the structural model is of interest.

...

Named arguments to be passed to lsoda. Can be used for example to pass events or error tolerances.

Value

A data frame with columns 't' for time, 'i' for the index of the output element, 'y' for the outputs, 'dy/d<v1>' for the first derivatives and (if secOrd == TRUE) 'd2y/d<v1>_<v2>' for the second derivatives. In the column names, variables v1 and v2 are replaced by the names in vartheta, vareta and vareps. The data frame has attributes 'theta', 'eta' and 'eps' listing the variables theta, eta and eps, as named vectors. The 'type' attribute is set to 'VariationalMatrix', and the 'secOrd' and 'symbolic' ones to the values of secOrd and symbolic, respectively. The variational matrix is not normalized, as reflected in the attributes 'parameterNormalized' and 'outputNormalized', that record input and output normalization, respectively. The attribute 'parameterNormalized' is set to a named boolean vector of FALSE, with the theta, eta and eps parameter names as names. The attribute 'outputNormalized' is set to FALSE. The function displays an error and returns NULL if vartheta contains elements not in names(theta), and likewise for eta and eps, or if times includes a negative value.

Details

If symbolic==TRUE, then the variational matrix is calculated by first taking symbolic derivatives of the differential equations and associated functions (p, init and output), and solving the resulting system of differential equations numerically. If symbolic==FALSE, then the different equations are solved numerically, and subsequently their derivatives are computed numerically (using Richardson's method, https://en.wikipedia.org/wiki/Richardson_extrapolation). The first method tends to give more stable results, and hence is set as default. It is not clear how the complexities of these methods compare. The first requires solving a system of differential equations of higher dimension, while the second requires multiple solutions to the original system (of lower dimension).

Note

It is recommended to use symbolic derivation (symbolic==TRUE) rather than numerical derivation, as it is more accurate, especially for derivatives that are formally zero. In numerical calculation, their value may be estimated as close to zero rather than exactly zero, and this may give the wrong answers, especially for the sensitivity and FIM methods. Downsides of symbolic derivation are that it can be considerably slower than the numerical option, and it places more constraints on the definition of the functions model, p, init and output. Therefore it is advisable to compare symbolic and numerical versions, in order to spot any mistakes in these definitions. The constraints are necessary so that derivatives are properly computed. They are checked (to an extent) in case chkModel = TRUE. They are:

  • The function arguments of these functions should be named exactly as specified in the arguments above. I.e., it is wrong to define model <- function(t, x, q) ....

  • Do not use with in the function definitions.

  • Do not use a return statement in the function definitions. It is allowed to use composite function bodies (i.e., containing multiple statements). The final statement is the return value.

  • Refer to function arguments only component-wise, i.e., as y[["a"]], y[[1]] etc, and not as vectors. A function definition of the form p <- function(theta, eta) theta should not be used. Rather use p <- function(theta, eta) c(theta[["CL"]], theta[["V"]]). Components may be referred to by [[]], but not by [] or $.

  • Components of the state vector should be referred to either by index or by name, but not both, also not between functions. Likewise for components of the parameter vectors.

  • Do not define any functions with names starting with "FU_". These are used internally to represent state variables and therefore may not be used for other purposes.

The functions model, p, init and output should satisfy the rules imposed by lsoda, where output corresponds to the additional output of the lsoda function func.

If this function is used for sensitivity or aliasing calculations, then individual and random parameters are not used, so the function p and output may ignore these inputs, and eta and eps may be set to NULL. Furthermore, in this case secOrd = FALSE. See calcFirstVariations.

If this function is used to calculate the Fisher Information matrix, then derivatives are evaluated at eta = 0 and eps = 0, so the values of these inputs should be set to 0. Furthermore, in this case secOrd = TRUE. See calcVariationsFim.

Author

Martijn van Noort