Description
Given a model consisting of differential equations, estimates the local effect of certain parameters on selected sensitivity variables by calculating a matrix of so-called sensitivity functions. In this matrix the (i,j)-th element contains
$$\frac{\partial y_i}{\partial \Theta _j}\cdot \frac{\Delta \Theta _j} {\Delta y_i}$$
and where \(y_i\) is an output variable (at a certain time instance), \(\Theta_j\) is a parameter, and \(\Delta y_i\) is the scaling of variable \(y_i\), \(\Delta \Theta_j\) is the scaling of parameter \(\Theta_j\).
Usage
sensFun(func, parms, sensvar = NULL, senspar = names(parms), varscale = NULL, parscale = NULL, tiny = 1e-8, map = 1, ...)# S3 method for sensFunsummary(object, vars = FALSE, ...)
# S3 method for sensFunpairs(x, which = NULL, ...)
# S3 method for sensFunplot(x, which = NULL, legpos="topleft", ask = NULL, ...)
# S3 method for summary.sensFunplot(x, which = 1:nrow(x), ...)
Value
a data.frame of class sensFun
containing the sensitivity functions this is one row for each sensitivity variable at each independent (time or position) value and the following columns:
x
, the value of the independent (mapping) variable, usually time (solver= "ode.."), or distance (solver= "steady.1D")
var
, the name of the observed variable,
...
, a number of columns, one for each sensitivity parameter
The data.frame returned by sensFun
has methods for the generic functions summary
, plot
,
pairs
-- see note.
Arguments
an R-function that has as first argument parameters passed to the output variables for which the sensitivity needs to be estimated. Either the parameters whose sensitivity needs to be estimated, the default=all parameters. Either a vector with parameter names, or a vector with indices to positions of parameters in the scaling (weighing) factor for sensitivity variables, the scaling (weighing) factor for sensitivity parameters, the perturbation, or numerical difference, factor, see details. the column number with the (independent) mapping variable in the output matrix returned by additional arguments passed to an object of class an object of class if FALSE: summaries per parameter are returned; if the name or the index to the variables that should be plotted. Default = all variables. position of the legend; set to logical; if parms
and that returns a matrix or data.frame with the values of the output variables (columns) at certain output intervals (rows), and -- optionally -- a mapping variable (by default the first column).func
; should be either a vector, or a list with named elements. If NULL
, then the first element of parInput
is taken.NULL
, the default, which selects all variables, or a vector with variable names
(which should be present in the matrix returned by func
), or a vector with indices
to variables as present in the output matrix (note that the column of this matrix with the mapping variable should not be selected).parms
.NULL
indicates that the variable value is used.NULL
indicates that the parameter value is used.func
. For dynamic models solved by integration, this will be the (first) column with time
. For 1-D spatial output, this column will be some distance variable. Set to NULL if there is no mapping variable. Mapping variables should not be selected for estimating sensitivity functions; they are used for plotting.func
or to the methods.sensFun
.sensFun
.TRUE
, summaries per parameter and per variable are returned.NULL
to avoid plotting a legend.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
.
Author
Karline Soetaert <karline.soetaert@nioz.nl>
Details
There are essentially two ways in which to use function sensFun
.
When
func
returns a matrix or data frame with output values,sensFun
can be used for sensitivity analysis, estimating the impact of parameters on output variables.When
func
returns an instance of classmodCost
(as returned by a call to functionmodCost
), thensensFun
can be used for parameter identifiability. In this case the results fromsensFun
are used as input to function collin. See the help file forcollin
.
For each sensitivity parameter, the number of sensitivity functions estimated is: length(sensvar) * length(mapping variable), i.e. one for each element returned by func
(except the mapping variable).
The sensitivity functions are estimated numerically. This means that each parameter value \(\Theta_j\) is perturbed as \(\max{(tiny,\Theta_j \cdot (1+tiny))}\)
References
Soetaert, K. and Herman, P. M. J., 2009. A Practical Guide to Ecological Modelling -- Using R as a Simulation Platform. Springer, 390 pp.
Brun, R., Reichert, P. and Kunsch, H.R., 2001. Practical Identificability Analysis of Large Environmental Simulation Models. Water Resour. Res. 37(4): 1015--1030. tools:::Rd_expr_doi("10.1029/2000WR900350")
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. 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, ylim = range(c(out$Bact, out$Sub)), xlab = "time, hour", ylab = "molC/m3", type = "l", lwd = 2)lines(out$time, out$Sub, lty = 2, lwd = 2)lines(out$time, out$Sub + out$Bact)legend("topright", c("Bacteria", "Glucose", "TOC"), lty = c(1, 2, 1), lwd = c(2, 2, 1))## sensitivity functionsSnsBact <- sensFun(func = solveBact, parms = pars, sensvar = "Bact", varscale = 1)head(SnsBact)plot(SnsBact)plot(SnsBact, type = "b", pch = 15:19, col = 2:6, main = "Sensitivity all vars")summary(SnsBact)plot(summary(SnsBact))SF <- sensFun(func = solveBact, parms = pars, sensvar = c("Bact", "Sub"), varscale = 1)head(SF)tail(SF)summary(SF, var = TRUE)plot(SF)plot(SF, which = c("Sub","Bact"))pm <- par(mfrow = c(1,3))plot(SF, which = c("Sub", "Bact"), mfrow = NULL)plot(SF, mfrow = NULL)par(mfrow = pm)## Bivariate sensitivitypairs(SF) # same colorpairs(SF, which = "Bact", col = "green", pch = 15)pairs(SF, which = c("Bact", "Sub"), col = c("green", "blue"))mtext(outer = TRUE, side = 3, line = -2, "Sensitivity functions", cex = 1.5)## pairwise correlationcor(SnsBact[,-(1:2)])
Run the code above in your browser using DataLab