ODEsobol.default {ODEsensitivity}

Sobol' Sensitivity Analysis for General ODE Models

Description

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

Usage

## 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,
  ...)

Arguments

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.

Details

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).

Value

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").

Note

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.

Author(s)

Stefan Theers, Frank Weber

References

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.

See Also

soboljansen, sobolmartinez, plot.ODEsobol

Examples

##### 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)