plot.ODEmorris: Plot of the Results of Morris Screening for Objects of Class `ODEmorris`

Description

plot.ODEmorris plots the results of Morris screening for objects of class ODEmorris.

Usage

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

Arguments

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

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.

kind

[character(1)] kind of the plot, choose between "sep" and "trajec" (see details).

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)] title for the plot. If kind = "sep", this is the overall title for the two separate plots. If NULL (the default), a standard 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.

Value

TRUE (invisible; for testing purposes).

Details

Morris sensitivity indices are plotted for one state variable (chosen by argument state_plot) and the parameters named in pars_plot. If no parameters are named in pars_plot, the sensitivity indices for all parameters are plotted. There are two kinds of plots:

kind = "sep": separate plots of the Morris sensitivity indices \(\mu^*\) and \(\sigma\) against time
kind = "trajec": plot of \(\mu^*\) against \(\sigma\)

Examples

Run this code

# NOT RUN {
##### 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(7292)
# Warning: The following code might take very long!
# }
# NOT RUN {
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)
my_cols <- c("firebrick", "orange2", "dodgerblue", 
             "forestgreen", "black")
plot(LVres_morris, kind = "sep", colors_pars = my_cols)
plot(LVres_morris, pars_plot = c("rGrow", "rMort"), state_plot = "Predator", 
     kind = "trajec", colors_pars = my_cols[2:3])
# }
# NOT RUN {
##### 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(283)
# Warning: The following code might take very long!
# }
# NOT RUN {
LFOres_morris <- ODEmorris(mod = lfonet,
                           pars = LFOpars,
                           times = LFOtimes,
                           binf = LFObinf,
                           bsup = LFObsup,
                           r = 500,
                           design = list(type = "oat", 
                                         levels = 10, grid.jump = 1),
                           scale = TRUE,
                           parallel_eval = TRUE,
                           parallel_eval_ncores = 2)
plot(LFOres_morris, pars_plot = paste0("k.", 1:4), state_plot = "x.2",
     kind = "sep", colors_pars = my_cols)
# }
# NOT RUN {
# }