FME: modCRL – R documentation

Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!

modCRL

Monte Carlo Analysis

Description

Given a model consisting of differential equations, estimates the global effect of certain (sensitivity) parameters on selected sensitivity variables.

This is done by drawing parameter values according to some predefined distribution, running the model with each of these parameter combinations, and calculating the values of the selected output variables at each output interval.

This function is useful for “what-if” scenarios.

If the output variables consist of a time-series or spatially dependent, use sensRange instead.

Usage

modCRL(func, parms = NULL, sensvar = NULL, dist = "unif",
       parInput = NULL, parRange = NULL, parMean = NULL, parCovar = NULL,
       num = 100, ...)

## S3 method for class 'modCRL'
summary(object,  ...)

## S3 method for class 'modCRL'
plot(x, which = NULL, trace = FALSE, ask = NULL, ...)

## S3 method for class 'modCRL'
pairs(x, which = 1:ncol(x), nsample = NULL, ...)

## S3 method for class 'modCRL'
hist(x, which = 1:ncol(x), ask = NULL, ...)

Arguments

`func`	an R-function that has as first argument `parms` and that returns a vector with variables whose sensitivity should be estimated.
`parms`	parameters passed to `func`; should be either a vector, or a list with named elements. If `NULL`, then the first element of `parInput` is taken.
`sensvar`	the output variables for which the sensitivity needs to be estimated. Either `NULL`, the default=all output variables, or a vector with output variable `names` (which should be present in the vector returned by `func`), or a vector with `indices` to output variables as present in the output vector returned by `func`.
`dist`	the distribution according to which the parameters should be generated, one of `"unif"` (uniformly random samples), `"norm"`, (normally distributed random samples), `"latin"` (latin hypercube distribution), `"grid"` (parameters arranged on a grid). The input parameters for the distribution are specified by `parRange` (min,max), except for the normally distributed parameters, in which case the distribution is specified by the parameter means `parMean` and the variance-covariance matrix, `parCovar`. Note that, if the distribution is `"norm"` and `parRange` is given, then a truncated distribution will be generated. (This is useful to prevent for instance that certain parameters become negative). Ignored if `parInput` is specified.
`parRange`	the range (min, max) of the sensitivity parameters, a matrix or (preferred) a data.frame with one row for each parameter, and two columns with the minimum (1st) and maximum (2nd) value. The rownames of `parRange` should be parameter names that are known in argument `parms`. Ignored if `parInput` is specified.
`parInput`	a matrix with dimension (*, npar) with the values of the sensitivity parameters.
`parMean`	only when `dist` is `"norm"`: the mean value of each parameter. Ignored if `parInput` is specified.
`parCovar`	only when `dist` is `"norm"`: the parameter's variance-covariance matrix.
`num`	the number of times the model has to be run. Set large enough. If `parInput` is specified, then `num` parameters are selected randomly (from the rows of `parInput`.
`object`	an object of class `modCRL`.
`x`	an object of class `modCRL`.
`which`	the name or the index to the variables and parameters that should be plotted. Default = all variables and parameters.
`nsample`	the number of xy pairs to be plotted on the upper panel in the pairs plot. When `NULL` all xy pairs plotted. Set to a lower number in case the graph becomes too dense (and the exported picture too large). This does not affect the histograms on the diagonal plot (which are estimated using all values).
`trace`	if `TRUE`, adds smoothed line to the plot.
`ask`	logical; if `TRUE`, the user is asked before each plot, if `NULL` the user is only asked if more than one page of plots is necessary and the current graphics device is set interactive, see `par(ask=.)` and `dev.interactive`.
`...`	additional arguments passed to function `func` or to the methods.

Value

a data.frame of type modCRL containing the parameter(s) and the corresponding values of the sensitivity output variables.

The list returned by modCRL has a method for the generic functions summary and plot – see note.

Note

The following methods are included:

summary, estimates summary statistics for the sensitivity variables, a table with as many rows as there are variables (or elements in the vector returned by func) and the following columns: x, the mapping value, Mean, the mean, sd, the standard deviation, Min, the minimal value, Max, the maximal value, q25, q50, q75, the 25th, 50 and 75% quantile.
plot, produces a plot of the modCRL output, either one plot for each sensitivity variable and with the parameter value on the x-axis. This only works when there is only one parameter!

OR

one plot for each parameter value on the x-axis. This only works when there is only one variable!
hist, produces a histogram of the modCRL output parameters and variables.
pairs, produces a pairs plot of the modCRL output.

The data.frame of type modCRL has several attributes, which remain hidden, and which are generally not of practical use (they are needed for the S3 methods). There is one exception - see notes in help of sensRange.

Author(s)

Karline Soetaert <karline.soetaert@nioz.nl>.

References

Soetaert, K. and Petzoldt, T., 2010. Inverse Modelling, Sensitivity and Monte Carlo Analysis in R Using Package FME. Journal of Statistical Software 33(3) 1–28. http://www.jstatsoft.org/v33/i03

Examples

## =======================================================================
## Bacterial growth model as in Soetaert and Herman, 2009
## =======================================================================

pars <- list(gmax = 0.5,eff = 0.5,
              ks = 0.5, rB = 0.01, dB = 0.01)

solveBact <- function(pars) {
  derivs <- function(t,state,pars) {    # returns rate of change
    with (as.list(c(state,pars)), {
      dBact <-  gmax*eff * Sub/(Sub + ks)*Bact - dB*Bact - rB*Bact
      dSub  <- -gmax     * Sub/(Sub + ks)*Bact + dB*Bact
      return(list(c(dBact, dSub)))
    })
  }

 state <- c(Bact = 0.1, Sub = 100)
 tout  <- seq(0, 50, by = 0.5)
 ## ode solves the model by integration...
 return(as.data.frame(ode(y = state, times = tout, func = derivs,
                          parms = pars)))
}

out <- solveBact(pars)

plot(out$time, out$Bact, main = "Bacteria",
  xlab = "time, hour", ylab = "molC/m3", type = "l", lwd = 2)

## Function that returns the last value of the simulation
SF <- function (p) {
  pars["eff"] <- p
  out <- solveBact(pars)
  return(out[nrow(out), 2:3])
}

parRange <- matrix(nr = 1, nc = 2, c(0.2, 0.8),
  dimnames = list("eff", c("min", "max")))
parRange

CRL <- modCRL(func = SF, parRange = parRange)

plot(CRL)  # plots both variables
plot(CRL, which = c("eff", "Bact"), trace = FALSE) #selects one

FME

A Flexible Modelling Environment for Inverse Modelling, Sensitivity, Identifiability and Monte Carlo Analysis

v1.3.6.1

GPL (>= 2)

Authors

Karline Soetaert [aut, cre] (<https://orcid.org/0000-0003-4603-7100>), Thomas Petzoldt [aut] (<https://orcid.org/0000-0002-4951-6468>)

Initial release