ODEmorris.default {ODEsensitivity}

Morris Screening for General ODE Models

Description

ODEmorris.default is the default method of ODEmorris. It performs a sensitivity analysis for general ODE models using the Morris screening method.

Usage

## Default S3 method:
ODEmorris(mod, pars, state_init, times, binf = 0,
  bsup = 1, r = 500, design = list(type = "oat", levels = 10, grid.jump =
  1), scale = TRUE, 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 Morris screening.

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.

binf

[character(1 or k)]
vector of lower borders of possible input parameter values. If they are all equal, a single value can be set.

bsup

[character(1 or k)]
vector of upper borders of possible input parameter values. If they are all equal, a single value can be set.

r

[integer(1 or 2)]
if of length 1, the number of repetitions of the design. If of length 2, a space-filling optimization of the sampling design is used, see morris. However, this space-filling optimization might lead to long runtimes, so length 1 is recommended for r. Defaults to 500.

design

[list]
a list specifying the design type and its parameters, cf. morris.

scale

[logical(1)]
if TRUE, scaling is done for the input design of experiments after building the design and before calculating the elementary effects, cf. morris. Defaults to TRUE, which is highly recommended if the factors have different orders of magnitude, see morris.

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 using morris from the package sensitivity.

For non-ODE models, values for r are typically between 10 and 50. However, much higher values are recommended for ODE models (the default is r = 500).

Value

List of class ODEmorris of length length(state_init) containing in each element a matrix for one state variable. The matrices themselves contain the Morris screening results for all timepoints (rows: mu, mu.star and sigma for every parameter; columns: timepoints).

Note

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

If morris throws a warning message stating "In ... keeping ... repetitions out of ...", try using a bigger number of levels in the design argument (only possible for OAT design).

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

morris, plot.ODEmorris

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))
# Morris screening:
set.seed(7292)
# Warning: The following code might take very long!

LVres_morris <- ODEmorris(mod = LVmod,
                          pars = LVpars,
                          state_init = LVinit,
                          times = LVtimes,
                          binf = LVbinf,
                          bsup = LVbsup,
                          r = 500,
                          design = list(type = "oat", 
                                        levels = 10, grid.jump = 1),
                          scale = TRUE,
                          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_morris <- ODEmorris(mod = FHNmod,
                           pars = c("a", "b", "s"),
                           state_init = c(Voltage = -1, Current = 1),
                           times = seq(0.1, 50, by = 5),
                           binf = c(0.18, 0.18, 2.8),
                           bsup = c(0.22, 0.22, 3.2),
                           r = 500,
                           design = list(type = "oat", 
                                         levels = 50, grid.jump = 1),
                           scale = TRUE,
                           ode_method = "adams",
                           parallel_eval = TRUE,
                           parallel_eval_ncores = 2)