sensRange: Sensitivity Ranges of a Timeseries or 1-D Variables

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

Sensitivity Ranges of a Timeseries or 1-D Variables

Description

Given a model consisting of differential equations, estimates the global effect of certain (sensitivity) parameters on a time series or on 1-D spatial series of selected sensitivity variables.

This is done by drawing parameter values according to some predefined distribution, running the model with each of these parameter combinations, and calculating the values of the selected output variables at each output interval.

This function thus produces 'envelopes' around the sensitivity variables.

Usage

sensRange(func, parms = NULL, sensvar = NULL, dist = "unif",
          parInput = NULL, parRange = NULL, parMean = NULL, 
          parCovar = NULL, map = 1, num = 100,  ...)
  
## S3 method for class 'sensRange'
summary(object, ...)

## S3 method for class 'summary.sensRange'
plot(x, xyswap = FALSE,
              which = NULL, legpos = "topleft",
              col = c(grey(0.8), grey(0.7)),
              quant = FALSE, ask = NULL, obs = NULL, 
              obspar = list(), ...)

## S3 method for class 'sensRange'
plot(x, xyswap = FALSE,
              which = NULL, ask = NULL, ...)

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).
dist	the distribution according to which the parameters should be generated, one of "unif" (uniformly random samples), "norm", (normally distributed random samples), "latin" (latin hypercube distribution), "grid" (parameters arranged on a grid). The input parameters for the distribution are specified by parRange (min,max), except for the normally distributed parameters, in which case the distribution is specified by the parameter means parMean and the variance-covariance matrix, parCovar. Note that, if the distribution is "norm" and parRange is given, then a truncated distribution will be generated. (This is useful to prevent for instance that certain parameters become negative). Ignored if parInput is specified.
parRange	the range (min, max) of the sensitivity parameters, a matrix or (preferred) a data.frame with one row for each parameter, and two columns with the minimum (1st) and maximum (2nd) value. The rownames of parRange should be parameter names that are known in argument parms. Ignored if parInput is specified.
parInput	a matrix with dimension (*, npar) with the values of the sensitivity parameters.
parMean	only when dist is "norm": the mean value of each parameter. Ignored if parInput is specified.
parCovar	only when dist is "norm": the parameter's variance-covariance matrix.
num	the number of times the model has to be run. Set large enough. If parInput is specified, then num parameters are selected randomly (from the rows of parInput.
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 ranges; they are used for plotting.
object	an object of class sensRange.
x	an object of class sensRange.
legpos	position of the legend; set to NULL to avoid plotting a legend.
xyswap	if TRUE, then x-and y-values are swapped and the y-axis is from top to bottom. Useful for drawing vertical profiles.
which	the name or the index to the variables that should be plotted. Default = all variables.
col	the two colors of the polygons that should be plotted.
quant	if TRUE, then the median surrounded by the quantiles q25-q75 and q95-q95 are plotted, else the min-max and mean +- sd are plotted.
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.
obs	a data.frame or matrix with "observed data" that will be added as points to the plots. obs can also be a list with multiple data.frames and/or matrices containing observed data. The first column of obs should contain the time or space-variable. If obs is not NULL and which is NULL, then the variables, common to both obs and x will be plotted.
obspar	additional graphics arguments passed to points, for plotting the observed data. If obs is a list containing multiple observed data sets, then the graphics arguments can be a vector or a list (e.g. for xlim, ylim), specifying each data set separately.
...	additional arguments passed to func or to the methods.

Details

Models solved by integration (i.e. by using one of 'ode', 'ode.1D', 'ode.band', 'ode.2D'), have the output already in a form usable by sensRange.

Value

a data.frame of type sensRange containing the parameter set and the corresponding values of the sensitivity output variables.

The list returned by sensRange has a method for the generic functions summary,plot and plot.summary – see note.

Note

The following methods are included:

• summary, estimates summary statistics for the sensitivity variables, a data.frame with as many rows as there are mapping variables (or rows in the matrix returned by func) and the following columns: x, the mapping value, Mean, the mean, sd, the standard deviation, Min, the minimal value, Max, the maximal value, q25, q50, q75, the 25th, 50 and 75% quantile

• plot, produces a "matplot" of the sensRange output, one plot for each sensitivity variable and with the mapping variable on the x-axis.

  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.

• summary.plot, produces a plot of the summary of the sensRange output, one plot for each sensitivity variable and with the ranges and mean +- standard deviation or the quantiles as coloured polygons.

  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.

The output for models solved by a steady-state solver (i.e. one of 'steady', 'steady.1D', 'steady.band', 'steady.2D', needs to be rearranged – see examples.

For plot.summary.sensRange and plot.sensRange, the number of panels per page is automatically determined up to 3 x 3 (par(mfrow = c(3, 3))). This default can be overwritten by specifying user-defined settings for mfrow or mfcol. Set mfrow equal to NULL to avoid the plotting function to change user-defined mfrow or mfcol settings.

Other graphical parameters can be passed as well. Parameters are vectorized, either according to the number of plots (xlab, ylab, main, sub, xlim, ylim, log, asp, ann, axes, frame.plot,panel.first,panel.last, cex.lab,cex.axis,cex.main) or according to the number of lines within one plot (other parameters e.g. col, lty, lwd etc.) so it is possible to assign specific axis labels to individual plots, resp. different plotting style. Plotting parameter ylim, or xlim can also be a list to assign different axis limits to individual plots.

Similarly, the graphical parameters for observed data, as passed by obspar can be vectorized, according to the number of observed data sets (when obs is a list).

The data.frame of type sensRange has several attributes, which remain hidden, and which are generally not of practical use (they are needed for the S3 methods).

There is one exception, i.e. if parameter values are imposed via argument parInput, and these parameters are generated by a Markov chain (modMCMC). If the number of draws, num, is less than the number of rows in parInput, then num random draws will be taken. Attribute, "pset" then contains the index to the parameters that have been selected.

The sensRange method only represents the distribution of the model response variables as a function of the parameter values. But an additional source of noise is due to the model error, as represented by the sampled values of sigma in the Markov chain. In order to represent also this source of error, gaussian noise should be added to each sensitivity output variables, with a standard deviation that corresponds to the original parameter draw – see vignette "FMEother".

Author(s)

Karline Soetaert <[e-mail]>

References

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 from 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)

mf  <-par(mfrow = c(2,2))

plot(out$time, out$Bact, main = "Bacteria",
     xlab = "time, hour", ylab = "molC/m3", type = "l", lwd = 2)

## the sensitivity parameters
parRanges <- data.frame(min = c(0.4, 0.4, 0.0), max = c(0.6, 0.6, 0.02))
rownames(parRanges)<- c("gmax", "eff", "rB")
parRanges

tout <- 0:50
## sensitivity to rB; equally-spaced parameters ("grid")
SensR <- sensRange(func = solveBact, parms = pars, dist = "grid",
                   sensvar = "Bact", parRange = parRanges[3,], num = 50)

Sens  <-summary(SensR)
plot(Sens, legpos = "topleft", xlab = "time, hour", ylab = "molC/m3",
     main = "Sensitivity to rB", mfrow = NULL)

## sensitivity to all; latin hypercube
Sens2 <- summary(sensRange(func = solveBact, parms = pars, dist = "latin",
           sensvar = c("Bact", "Sub"), parRange = parRanges, num = 50))

## Plot all variables; plot mean +- sd, min max
plot(Sens2, xlab = "time, hour", ylab = "molC/m3",
     main = "Sensitivity to gmax,eff,rB", mfrow = NULL)

par(mfrow = mf)

## Select one variable for plotting; plot the quantiles
plot(Sens2, xlab = "time, hour", ylab = "molC/m3", which = "Bact", quant = TRUE)

## Add data
data <- cbind(time = c(0,10,20,30), Bact = c(0,1,10,45))
plot(Sens2, xlab = "time, hour", ylab = "molC/m3", quant = TRUE, 
  obs = data, obspar = list(col = "darkblue", pch = 16, cex = 2))

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