Local Sensitivity Analysis
Description
Given a model consisting of differential equations, estimates the local effect of certain parameters on selected sensitivity variables by calculating a matrix of socalled sensitivity functions. In this matrix the (i,j)th element contains
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.
Usage
1 2 3 4 5 6 7 8 9 10 11 12 13 14  sensFun(func, parms, sensvar = NULL, senspar = names(parms),
varscale = NULL, parscale = NULL, tiny = 1e8, 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 Rfunction that has as first argument 
parms 
parameters passed to 
sensvar 
the output variables for which the sensitivity needs
to be estimated. Either 
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 
varscale 
the scaling (weighing) factor for sensitivity
variables, 
parscale 
the scaling (weighing) factor for sensitivity
parameters, 
tiny 
the perturbation, or numerical difference, factor, see details. 
map 
the column number with the (independent) mapping variable
in the output matrix returned by 
... 
additional arguments passed to 
object 
an object of class 
x 
an object of class 
vars 
if FALSE: summaries per parameter are returned; if

which 
the name or the index to the variables that should be plotted. Default = all variables. 
legpos 
position of the legend; set to 
ask 
logical; if 
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 par_j is perturbed as max(tiny,par_j)*(1+tiny)
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.
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 L1norm sum(abs(Sij))/n,
L2: the L2norm 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 parmfrow
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
, anddiag.panel
.
Author(s)
Karline Soetaert <karline.soetaert@nioz.nl>
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.
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 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65  ## =======================================================================
## 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)])
