# Local Sensitivity Analysis 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 [REMOVE_ME]$$ \frac{\partial y_i}{\partial \Theta _j}\cdot \frac{\Delta \Theta _j}{\Delta y_i}dy_i/dpar_j*parscale_j/varscale_i [REMOVE_ME_2]$$ and where $y_i$ is an output variable (at a certain time instance), $par_j$ is a parameter, and $varscale_i$ is the scaling of variable $y_i$, $parscale_j$ is the scaling of parameter $par_j$. ## 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}dy_i/dpar_j*parscale_j/varscale_i $$ and where $y_i$ is an output variable (at a certain time instance), $par_j$ is a parameter, and $varscale_i$ is the scaling of variable $y_i$, $parscale_j$ is the scaling of parameter $par_j$. ```r sensFun(func, parms, sensvar = NULL, senspar = names(parms), varscale = NULL, parscale = NULL, tiny = 1e-8, map = 1, ...) ## S3 method for class 'sensFun' summary(object, vars = FALSE, ...) ## S3 method for class 'sensFun' pairs(x, which = NULL, ...) ## S3 method for class 'sensFun' plot(x, which = NULL, legpos="topleft", ask = NULL, ...) ## S3 method for class 'summary.sensFun' plot(x, which = 1:nrow(x), ...) ``` ## Arguments - `func`: an R-function that has as first argument `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). - `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, 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). - `senspar`: 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 `parms`. - `varscale`: the scaling (weighing) factor for sensitivity variables, `NULL` indicates that the variable value is used. - `parscale`: the scaling (weighing) factor for sensitivity parameters, `NULL` indicates that the parameter value is used. - `tiny`: the perturbation, or numerical difference, factor, see details. - `map`: the column number with the (independent) mapping variable in the output matrix returned by `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. - `...`: additional arguments passed to `func` or to the methods. - `object`: an object of class `sensFun`. - `x`: an object of class `sensFun`. - `vars`: if FALSE: summaries per parameter are returned; if `TRUE`, summaries per parameter and per variable are returned. - `which`: the name or the index to the variables that should be plotted. Default = all variables. - `legpos`: position of the legend; set to `NULL` to avoid plotting a legend. - `ask`: logical; if `TRUE`, the user is **ask**ed 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`. ## Returns 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. ## 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 class `modCost` (as returned by a call to function `modCost`), then `sensFun` can be used for **parameter identifiability**. In this case the results from `sensFun` are used as input to function collin . See the help file for `collin`. 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 $par_j$ is perturbed as $max(tiny,par_j)*(1+tiny)$ ## Note Sensitivity functions are generated by perturbing one by one the parameters with a very small amount, and quantifying the differences in the output. It is important that the output is generated with high precision, else it is possible, that the sensitivity functions are just noise. For instance, when used with a dynamic model (using solver from `deSolve`) set the tolerances `atol` and `rtol` to a lower value, to see if the sensitivity results make sense. The following methods are provided: * **summary**. Produces summary statistics of the sensitivity functions, a `data.frame` with: one row for each parameter and the following columns: * L1: the L1-norm $sum(abs(Sij))/n$, * L2: the L2-norm $sqrt(sum(Sij^2)/n)$, * Mean: the mean of the sensitivity functions, * Min: the minimal value of the sensitivity functions, * Max: the maximal value of the sensitivity functions. * **var**the summary of the variables sensitivity functions, a data.frame with the same columns as `model` and one row for each parameter + variable combination. This is only outputted if the variable names are effectively known * **plot** plots the sensitivity functions for each parameter; each parameter has its own color. By default, the sensitivity functions for all variables are plotted in one figure, unless `which` gives a selection of variables; in that case, each variable will be plotted in a separate figure, and the figures aligned in a rectangular grid, unless par `mfrow` is passed as an argument. * **pairs** produces a pairs plot of the sensitivity results; per parameter. By default, the sensitivity functions for all variables are plotted in one figure, unless `which` gives a selection of variables. Overrides the default `gap = 0`, `upper.panel = NA`, and `diag.panel`. ## Examples ```r ## ======================================================================= ## 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 functions SnsBact <- 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 sensitivity pairs(SF) # same color pairs(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 correlation cor(SnsBact[,-(1:2)]) ``` ## Author(s) Karline Soetaert ## 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")

sensFun function