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
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
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 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.
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)
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 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.
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
## =======================================================================
## 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)])
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(FME)
Loading required package: deSolve
Attaching package: 'deSolve'
The following object is masked from 'package:graphics':
matplot
Loading required package: rootSolve
Loading required package: coda
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/FME/sensFun.Rd_%03d_medium.png", width=480, height=480)
> ### Name: sensFun
> ### Title: Local Sensitivity Analysis
> ### Aliases: sensFun summary.sensFun plot.sensFun plot.summary.sensFun
> ### pairs.sensFun
> ### Keywords: utilities
>
> ### ** 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 functions
> SnsBact <- sensFun(func = solveBact, parms = pars,
+ sensvar = "Bact", varscale = 1)
> head(SnsBact)
x var gmax eff ks rB dB
1 0.0 Bact 0.00000000 0.00000000 0.000000e+00 0.0000000000 0.0000000000
2 0.5 Bact 0.01394694 0.01394695 -7.024728e-05 -0.0005605677 -0.0005605676
3 1.0 Bact 0.03127202 0.03127206 -1.561543e-04 -0.0012570761 -0.0012570752
4 1.5 Bact 0.05259198 0.05259209 -2.623152e-04 -0.0021141409 -0.0021141387
5 2.0 Bact 0.07861772 0.07861795 -3.920489e-04 -0.0031603689 -0.0031603643
6 2.5 Bact 0.11017696 0.11017737 -5.494327e-04 -0.0044290424 -0.0044290341
> plot(SnsBact)
> plot(SnsBact, type = "b", pch = 15:19, col = 2:6,
+ main = "Sensitivity all vars")
>
> summary(SnsBact)
value scale L1 L2 Mean Min Max N
gmax 0.50 0.50 29.51 5.859 16.2 -17.1 266.360 101
eff 0.50 0.50 37.12 6.212 37.1 0.0 268.408 101
ks 0.50 0.50 0.17 0.037 -0.1 -1.8 0.097 101
rB 0.01 0.01 3.47 0.463 -3.5 -10.8 0.000 101
dB 0.01 0.01 2.06 0.297 -2.1 -10.8 0.000 101
> plot(summary(SnsBact))
>
> SF <- sensFun(func = solveBact, parms = pars,
+ sensvar = c("Bact", "Sub"), varscale = 1)
> head(SF)
x var gmax eff ks rB dB
1 0.0 Bact 0.00000000 0.00000000 0.000000e+00 0.0000000000 0.0000000000
2 0.5 Bact 0.01394694 0.01394695 -7.024728e-05 -0.0005605677 -0.0005605676
3 1.0 Bact 0.03127202 0.03127206 -1.561543e-04 -0.0012570761 -0.0012570752
4 1.5 Bact 0.05259198 0.05259209 -2.623152e-04 -0.0021141409 -0.0021141387
5 2.0 Bact 0.07861772 0.07861795 -3.920489e-04 -0.0031603689 -0.0031603643
6 2.5 Bact 0.11017696 0.11017737 -5.494327e-04 -0.0044290424 -0.0044290341
> tail(SF)
x var gmax eff ks rB dB
197 47.5 Sub -0.01041233 7.806256e-10 0.01020408 -1.561251e-11 0.01041233
198 48.0 Sub -0.01041233 6.938894e-10 0.01020408 -1.561251e-11 0.01041233
199 48.5 Sub -0.01041233 -9.540979e-10 0.01020408 2.428613e-11 0.01041233
200 49.0 Sub -0.01041233 3.469447e-10 0.01020408 0.000000e+00 0.01041233
201 49.5 Sub -0.01041233 8.673617e-11 0.01020408 -8.673617e-12 0.01041233
202 50.0 Sub -0.01041233 -3.469447e-10 0.01020408 3.469447e-12 0.01041233
>
> summary(SF, var = TRUE)
value scale L1 L2 Mean Min Max N var
gmax1 0.50 0.50 29.51 58.88 16.25 -1.7e+01 2.7e+02 101 Bact
gmax2 0.50 0.50 48.40 122.95 -48.40 -5.6e+02 0.0e+00 101 Sub
eff1 0.50 0.50 37.12 62.43 37.12 0.0e+00 2.7e+02 101 Bact
eff2 0.50 0.50 39.64 102.50 -39.64 -4.8e+02 6.9e-06 101 Sub
ks1 0.50 0.50 0.17 0.37 -0.10 -1.8e+00 9.7e-02 101 Bact
ks2 0.50 0.50 0.29 0.77 0.29 0.0e+00 3.8e+00 101 Sub
rB1 0.01 0.01 3.47 4.65 -3.47 -1.1e+01 0.0e+00 101 Bact
rB2 0.01 0.01 1.59 4.12 1.59 -2.8e-07 1.9e+01 101 Sub
dB1 0.01 0.01 2.06 2.98 -2.06 -1.1e+01 0.0e+00 101 Bact
dB2 0.01 0.01 1.78 4.54 1.78 0.0e+00 2.1e+01 101 Sub
>
> 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)])
gmax eff ks rB dB
gmax 1.0000000 0.9184218 -0.9879345 -0.2602264 -0.7165957
eff 0.9184218 1.0000000 -0.9265083 -0.5575108 -0.8883637
ks -0.9879345 -0.9265083 1.0000000 0.2878314 0.7302562
rB -0.2602264 -0.5575108 0.2878314 1.0000000 0.8599353
dB -0.7165957 -0.8883637 0.7302562 0.8599353 1.0000000
>
>
>
>
>
>
> dev.off()
null device
1
>