sensFun: Local Sensitivity Analysis

In FME: A Flexible Modelling Environment for Inverse Modelling, Sensitivity, Identifiability and Monte Carlo Analysis

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 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 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 \Theta_j is perturbed as \max{(tiny,\Theta_j \cdot (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 \frac{1}{n} \cdot \sum{|S_{ij}|},

  • L2: the L2-norm \cdot \sqrt{\frac{1}{n}\sum{S_{ij} \cdot S_{ij}}},

  • 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 <[e-mail]>

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. \Sexpr[results=rd]{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. \Sexpr[results=rd]{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 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)])

FME documentation built on July 9, 2023, 5:59 p.m.