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 “whatif” scenarios.
If the output variables consist of a timeseries or spatially dependent, use sensRange instead.
Usage
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15  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 Rfunction that has as first argument 
parms 
parameters passed to 
sensvar 
the output variables for which the sensitivity needs to be
estimated. Either 
dist 
the distribution according to which the parameters should be
generated, one of The input parameters for the distribution are specified by

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 
parInput 
a matrix with dimension (*, npar) with the values of the sensitivity parameters. 
parMean 
only when 
parCovar 
only when 
num 
the number of times the model has to be run. Set large enough.
If 
object 
an object of class 
x 
an object of class 
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 
trace 
if 
ask 
logical; if 
... 
additional arguments passed to function 
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 xaxis. This only works when there is only one parameter!OR
one plot for each parameter value on the xaxis. 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
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43  ## =======================================================================
## 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
