ODEsobol.ODEnetwork {ODEsensitivity}

Sobol' Sensitivity Analysis for Objects of Class ODEnetwork

Description

ODEsobol.ODEnetwork performs the variance-based Sobol' sensitivity analysis for objects of class ODEnetwork. Package ODEnetwork is required for this function to work.

Usage

## S3 method for class 'ODEnetwork'
ODEsobol(mod, pars, 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

[ODEnetwork]
list of class ODEnetwork.

pars

[character(k)]
names of the parameters to be included as input variables in the Sobol' sensitivity analysis. All parameters must be contained in names(ODEnetwork::createParamVec(mod)) and must not be derivable from other parameters supplied (e.g. "k.2.1" can be derived from "k.1.2", so supplying "k.1.2" suffices).

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 ODEs in situations where the solution has to be determined numerically, see ode for details. 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

If the object of class ODEnetwork supplied for mod doesn't include any events, the solution of the ODE network is determined analytically using simuNetwork. In the presence of events, simuNetwork uses ode to solve the ODE network numerically.

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 2 * nrow(mod$state) and of class ODEsobol containing in each element a list of the Sobol' sensitivity analysis results for the corresponding state 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

In situations where the solution of the ODE model has to be determined numerically, it might be helpful to try a different type of ODE-solver (argument ode_method) if the simulation of the model takes too long. 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)

Frank Weber

See Also

soboljansen, sobolmartinez, plot.ODEsobol

Examples

##### A network of 4 mechanical oscillators connected in a circle #####
# Definition of the network using the package "ODEnetwork":
M_mat <- rep(2, 4)
K_mat <- diag(rep(2 * (2*pi*0.17)^2, 4))
K_mat[1, 2] <- K_mat[2, 3] <- 
  K_mat[3, 4] <- K_mat[1, 4] <- 2 * (2*pi*0.17)^2 / 10
D_mat <- diag(rep(0.05, 4))
library("ODEnetwork")
lfonet <- ODEnetwork(masses = M_mat, dampers = D_mat, springs = K_mat)
# The parameters to be included in the sensitivity analysis and their lower
# and upper boundaries:
LFOpars <- c("k.1", "k.2", "k.3", "k.4",
             "d.1", "d.2", "d.3", "d.4")
LFObinf <- c(rep(0.2, 4), rep(0.01, 4))
LFObsup <- c(rep(20, 4), rep(0.1, 4))
# Setting of the initial values of the state variables:
lfonet <- setState(lfonet, state1 = rep(2, 4), state2 = rep(0, 4))
# The timepoints of interest:
LFOtimes <- seq(25, 150, by = 2.5)
# Sobol' sensitivity analysis (here only with n = 500, but n = 1000 is
# recommended):
set.seed(1739)
# Warning: The following code might take very long! There are warnings
# occurring which might be due to "n" being too low.

suppressWarnings(
  LFOres_sobol <- ODEsobol(mod = lfonet,
                           pars = LFOpars,
                           times = LFOtimes,
                           n = 500,
                           rfuncs = "runif",
                           rargs = paste0("min = ", LFObinf,
                                          ", max = ", LFObsup),
                           sobol_method = "Martinez",
                           parallel_eval = TRUE,
                           parallel_eval_ncores = 2)
)