modCRL: Monte Carlo Analysis

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

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 <[e-mail]>.

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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.18637/jss.v033.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 documentation built on July 9, 2023, 5:59 p.m.