ODEsobol.default: Sobol' Sensitivity Analysis for General ODE Models
In ODEsensitivity: Sensitivity Analysis of Ordinary Differential Equations

Description Usage Arguments Details Value Note Author(s) References See Also Examples

ODEsobol.default is the default method of ODEsobol. It performs the variance-based Sobol' sensitivity analysis for general ODE models.

## Default S3 method:
ODEsobol(mod, pars, state_init, times, n = 1000,
  rfuncs = "runif", rargs = "min = 0, max = 1", sobol_method = "Martinez",
  ode_method = "lsoda", parallel_eval = FALSE, parallel_eval_ncores = NA,
  ...)

`mod`	[`function(Time, State, Pars)`] model to examine, supplied in the manner as needed for `ode` (see example below).
`pars`	[`character(k)`] names of the parameters to be included as input variables in the Sobol' sensitivity analysis.
`state_init`	[`numeric(z)`] vector of `z` initial values. Must be named (with unique names).
`times`	[`numeric`] points of time at which the sensitivity analysis should be executed (vector of arbitrary length). The first point of time must be greater than zero.
`n`	[`integer(1)`] number of random parameter values used to estimate the Sobol' sensitivity indices by Monte Carlo simulation. Defaults to 1000.
`rfuncs`	[`character(1` or `k)`] names of the functions used to generate the `n` random values for the `k` parameters. Can be of length 1 or `k`. If of length 1, the same function is used for all parameters. Defaults to `"runif"`, so a uniform distribution is assumed for all parameters.
`rargs`	[`character(1` or `k)`] arguments to be passed to the functions in `rfuncs`. Can be of length 1 or `k`. If of length 1, the same arguments are used for all parameters. Each element of `rargs` has to be a string of the form `"tag1 = value1, tag2 = value2, ..."`, see example below. Default is `"min = 0, max = 1"`, so (together with the default value of `rfuncs`) a uniform distribution on [0, 1] is assumed for all parameters.
`sobol_method`	[`character(1)`] either `"Jansen"` or `"Martinez"`, specifying which modification of the variance-based Sobol' method shall be used. Defaults to `"Martinez"`.
`ode_method`	[`character(1)`] method to be used for solving the differential equations, see `ode`. Defaults to `"lsoda"`.
`parallel_eval`	[`logical(1)`] logical indicating if the evaluation of the ODE model shall be performed parallelized.
`parallel_eval_ncores`	[`integer(1)`] number of processor cores to be used for parallelization. Only applies if `parallel_eval = TRUE`. If set to `NA` (as per default) and `parallel_eval = TRUE`, 1 processor core is used.
`...`	further arguments passed to or from other methods.

Function ode from deSolve is used to solve the ODE system.

The sensitivity analysis is done for all state variables and all timepoints simultaneously. If sobol_method = "Jansen", soboljansen from the package sensitivity is used to estimate the Sobol' sensitivity indices and if sobol_method = "Martinez", sobolmartinez is used (also from the package sensitivity).

List of length length(state_init) and of class ODEsobol containing in each element a list of the Sobol' sensitivity analysis results for the corresponding state_init-variable (i.e. first order sensitivity indices S and total sensitivity indices T) for every point of time in the times vector. This list has an extra attribute "sobol_method" where the value of argument sobol_method is stored (either "Jansen" or "Martinez").

If the evaluation of the model function takes too long, it might be helpful to try a different type of ODE-solver (argument ode_method). The ode_methods "vode", "bdf", "bdf_d", "adams", "impAdams" and "impAdams_d" might be faster than the standard ode_method "lsoda".

If n is too low, the Monte Carlo estimation of the sensitivity indices might be very bad and even produce first order indices < 0 or total indices > 1. First order indices in the interval [-0.05, 0) and total indices in (1, 1.05] are considered as minor deviations and set to 0 resp. 1 without a warning. First order indices < -0.05 or total indices > 1.05 are considered as major deviations. They remain unchanged and a warning is thrown. Up to now, first order indices > 1 or total indices < 0 haven't occured yet. If this should be the case, please contact the package author.

Stefan Theers, Frank Weber

J. O. Ramsay, G. Hooker, D. Campbell and J. Cao, 2007, Parameter estimation for differential equations: a generalized smoothing approach, Journal of the Royal Statistical Society, Series B, 69, Part 5, 741–796.

soboljansen, sobolmartinez, plot.ODEsobol

##### Lotka-Volterra equations #####
# The model function:
LVmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    Ingestion    <- rIng  * Prey * Predator
    GrowthPrey   <- rGrow * Prey * (1 - Prey/K)
    MortPredator <- rMort * Predator
    
    dPrey        <- GrowthPrey - Ingestion
    dPredator    <- Ingestion * assEff - MortPredator
    
    return(list(c(dPrey, dPredator)))
  })
}
# The parameters to be included in the sensitivity analysis and their lower 
# and upper boundaries:
LVpars  <- c("rIng", "rGrow", "rMort", "assEff", "K")
LVbinf <- c(0.05, 0.05, 0.05, 0.05, 1)
LVbsup <- c(1.00, 3.00, 0.95, 0.95, 20)
# The initial values of the state variables:
LVinit  <- c(Prey = 1, Predator = 2)
# The timepoints of interest:
LVtimes <- c(0.01, seq(1, 50, by = 1))
set.seed(59281)
# Sobol' sensitivity analysis (here only with n = 500, but n = 1000 is
# recommended):
# Warning: The following code might take very long!

LVres_sobol <- ODEsobol(mod = LVmod,
                        pars = LVpars,
                        state_init = LVinit,
                        times = LVtimes,
                        n = 500,
                        rfuncs = "runif",
                        rargs = paste0("min = ", LVbinf,
                                       ", max = ", LVbsup),
                        sobol_method = "Martinez",
                        ode_method = "lsoda",
                        parallel_eval = TRUE,
                        parallel_eval_ncores = 2)


##### FitzHugh-Nagumo equations (Ramsay et al., 2007) #####
FHNmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    
    dVoltage <- s * (Voltage - Voltage^3 / 3 + Current)
    dCurrent <- - 1 / s *(Voltage - a + b * Current)
    
    return(list(c(dVoltage, dCurrent)))
  })
}
# Warning: The following code might take very long!

FHNres_sobol <- ODEsobol(mod = FHNmod,
                         pars = c("a", "b", "s"),
                         state_init = c(Voltage = -1, Current = 1),
                         times = seq(0.1, 50, by = 5),
                         n = 500,
                         rfuncs = "runif",
                         rargs = c(rep("min = 0.18, max = 0.22", 2),
                                   "min = 2.8, max = 3.2"),
                         sobol_method = "Martinez",
                         ode_method = "adams",
                         parallel_eval = TRUE,
                         parallel_eval_ncores = 2)

# Just for demonstration purposes: The use of different distributions for the 
# parameters (here, the distributions and their arguments are chosen 
# completely arbitrarily):
# Warning: The following code might take very long!

demo_dists <- ODEsobol(mod = FHNmod,
                       pars = c("a", "b", "s"),
                       state_init = c(Voltage = -1, Current = 1),
                       times = seq(0.1, 50, by = 5),
                       n = 500,
                       rfuncs = c("runif", "rnorm", "rexp"),
                       rargs = c("min = 0.18, max = 0.22",
                                 "mean = 0.2, sd = 0.2 / 3",
                                 "rate = 1 / 3"),
                       sobol_method = "Martinez",
                       ode_method = "adams",
                       parallel_eval = TRUE,
                       parallel_eval_ncores = 2)