plot.ODEsobol {ODEsensitivity}

Plot of the Results of Sobol' Sensitivity Analysis for Objects of Class ODEsobol

Description

plot.ODEsobol plots the results of Sobol' SA for objects of class ODEsobol.

Usage

## S3 method for class 'ODEsobol'
plot(x, pars_plot = NULL, state_plot = names(x)[1],
  colors_pars = NULL, main_title = NULL, legendPos = "outside",
  type = "l", ...)

Arguments

x

[ODEsobol]
output of ODEsobol (of class ODEsobol).

pars_plot

[character(k)]
names of the k parameters to be plotted. If NULL (the default), all parameters are plotted.

state_plot

[character(1)]
name of the state variable to be plotted. Defaults to the name of the first state variable.

colors_pars

[character(>= k)]
vector of the colors to be used for the k different parameters. Must be at least of length k (only the first k elements will be used, though). If NULL (the default), rainbow(k) is used.

main_title

[character(1)]
common title for the two graphics. Default is NULL, which means an automatic title is generated.

legendPos

[character(1)]
keyword for the legend position, either one of those specified in legend or "outside" (the default), which means the legend is placed under the plot (useful, if there are many parameters in the model).

type

[character(1)]
plot type, i.e. "p", "l", "b", "c", "o", "s", "h" or "n". Defaults to "l".

...

additional arguments passed to plot.default.

Details

First order and total Sobol' sensitivity indices are plotted for one state variable (chosen by argument state_plot) and the parameters named in pars_plot against time. If no parameters are named in pars_plot, the sensitivity indices for all parameters are plotted.

Value

TRUE (invisible; for testing purposes).

Note

Not all arguments of plot.default can be passed by ..., for example xlab and ylab are fixed.

Author(s)

Stefan Theers, Frank Weber

See Also

ODEsobol, soboljansen, sobolmartinez

Examples

##### Lotka-Volterra equations #####
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)))
  })
}
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)
LVinit  <- c(Prey = 1, Predator = 2)
LVtimes <- c(0.01, seq(1, 50, by = 1))
set.seed(59281)
# 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)
my_cols <- c("firebrick", "orange2", "dodgerblue", 
             "forestgreen", "black")
plot(LVres_sobol, colors_pars = my_cols)
plot(LVres_sobol, pars_plot = c("rGrow", "rMort"), state_plot = "Predator", 
     colors_pars = my_cols[2:3])


##### A network of 4 mechanical oscillators connected in a circle #####
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)
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))
lfonet <- setState(lfonet, state1 = rep(2, 4), state2 = rep(0, 4))
LFOtimes <- seq(25, 150, by = 2.5)
set.seed(1739)
# Warning: The following code might take very long!

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)
)
plot(LFOres_sobol, pars_plot = paste0("k.", 1:4), state_plot = "x.2",
     colors_pars = my_cols)