Title: | A Flexible Modelling Environment for Inverse Modelling, Sensitivity, Identifiability and Monte Carlo Analysis |
---|---|
Description: | Provides functions to help in fitting models to data, to perform Monte Carlo, sensitivity and identifiability analysis. It is intended to work with models be written as a set of differential equations that are solved either by an integration routine from package 'deSolve', or a steady-state solver from package 'rootSolve'. However, the methods can also be used with other types of functions. |
Authors: | Karline Soetaert [aut, cre] , Thomas Petzoldt [aut] |
Maintainer: | Karline Soetaert <[email protected]> |
License: | GPL (>= 2) |
Version: | 1.3.6.3 |
Built: | 2024-11-30 06:34:57 UTC |
Source: | https://github.com/r-forge/fme |
R-package FME contains functions to run complex applications of models that produce output as a function of input parameters.
Although it was created to be used with models consisting of ordinary differential equations (ODE), partial differential equations (PDE) or differential algebraic equations (DAE), it can work with other models.
It contains:
Functions to allow fitting of the model to data.
Function modCost
estimates the (weighted) residuals between
model output and data, variable and model costs.
Function modFit
uses the output of modCost
to find
the best-fit parameters. It provides a wrapper around R's
built-in minimisation routines (optim
, nlm
,
nlminb
) and nls.lm
from package minpack.lm
.
Package FME also includes an implementation of the pseudo-random
search algorithm (function pseudoOptim
).
Function sensFun
estimates the sensitivity functions of
selected output variables as a function of model parameters. This
is the basis of uni-variate, bi-variate and multi-variate
sensitivity analysis.
Function collin
uses as input the sensitivity functions
and estimates the "collinearity" index for all possible parameter
sets. This multivariate sensitivity estimate measures approximate
linear dependence and is useful to derive which parameter sets are
identifiable given the data set.
Function sensRange
produces 'envelopes' around the
sensitivity variables, consisting of a time series or a
1-dimensional set, as a function of the sensitivity parameters.
It produces "envelopes" around the variables.
Function modCRL
calculates the values of single
variables as a function of the sensitivity parameters. This
function can be used to run simple "what-if" scenarios
Function modMCMC
runs a Markov chain Monte Carlo
(Bayesian analysis). It implements the delayed rejection -
adaptive Metropolis (DRAM) algorithm.
FME also contains functions to generate multiple parameter
values arranged according to a grid (Grid
) multinormal
(Norm
) or uniform (Unif
) design, and a latin
hypercube sampling (Latinhyper
) function
bug corrections:
version 1.3.6, sensFun: corrected calculation of L2 norm (now consistent with help page),
version 1.3, modCost: minlogp was not correctly estimated if more than one observed variable (used the wrong sd).
Karline Soetaert
Thomas Petzoldt
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. doi:10.18637/jss.v033.i03
## Not run: ## show examples (see respective help pages for details) example(modCost) example(sensFun) example(modMCMC) example(modCRL) ## open the directory with documents browseURL(paste(system.file(package = "FME"), "/doc", sep = "")) ## open the directory with examples browseURL(paste(system.file(package = "FME"), "/doc/examples", sep = "")) ## the vignettes vignette("FME") vignette("FMEdyna") vignette("FMEsteady") vignette("FMEother") vignette("FMEmcmc") edit(vignette("FME")) edit(vignette("FMEdyna")) edit(vignette("FMEsteady")) edit(vignette("FMEother")) edit(vignette("FMEmcmc")) ## End(Not run)
## Not run: ## show examples (see respective help pages for details) example(modCost) example(sensFun) example(modMCMC) example(modCRL) ## open the directory with documents browseURL(paste(system.file(package = "FME"), "/doc", sep = "")) ## open the directory with examples browseURL(paste(system.file(package = "FME"), "/doc/examples", sep = "")) ## the vignettes vignette("FME") vignette("FMEdyna") vignette("FMEsteady") vignette("FMEother") vignette("FMEmcmc") edit(vignette("FME")) edit(vignette("FMEdyna")) edit(vignette("FMEsteady")) edit(vignette("FMEother")) edit(vignette("FMEmcmc")) ## End(Not run)
Based on the sensitivity functions of model variables to a selection of parameters, calculates the "identifiability" of sets of parameter.
The sensitivity functions are a matrix whose (i,j)-th element contains
and where
is an output variable, at a certain (time) instance, i,
is the scaling of variable
,
is the scaling of
parameter
.
Function collin
estimates the collinearity, or identifiability of all
parameter sets or of one parameter set.
As a rule of thumb, a collinearity value less than about 20 is "identifiable".
collin(sensfun, parset = NULL, N = NULL, which = NULL, maxcomb = 5000) ## S3 method for class 'collin' print(x, ...) ## S3 method for class 'collin' plot(x, ...)
collin(sensfun, parset = NULL, N = NULL, which = NULL, maxcomb = 5000) ## S3 method for class 'collin' print(x, ...) ## S3 method for class 'collin' plot(x, ...)
sensfun |
model sensitivity functions as estimated by |
parset |
one selected parameter combination, a vector with their names or with the indices to the parameters. |
N |
the number of parameters in the set; if |
which |
the name or the index to the observed variables that should be used. Default = all observed variables. |
maxcomb |
the maximal number of combinations that can be tested.
If too large, this may produce a huge output. The number of combinations of
n parameters out of a total of p parameters is |
x |
an object of class |
... |
additional arguments passed to the methods. |
The collinearity is a measure of approximate linear dependence between sets of parameters. The higher its value, the more the parameters are related. With "related" is meant that several paraemter combinations may produce similar values of the output variables.
a data.frame of class collin
with one row for each parameter
combination (parameters as in sensfun
).
Each row contains:
... |
for each parameter whether it is present (1) or absent (0) in the set, |
N |
the number of parameters in the set, |
collinearity |
the collinearity value. |
The data.frame returned by collin
has methods for the generic
functions print
and plot
.
It is possible to use collin
for selecting parameter sets that
can be fine-tuned based on a data set. Thus it is a powerful
technique to make model calibration routines more robust, because
calibration routines often fail when parameters are strongly related.
In general, when the collinearity index exceeds 20, the linear dependence is assumed to be critical (i.e. it will not be possible or easy to estimate all the parameters in the combination together).
The procedure is explained in Omlin et al. (2001).
1. First the function collin
is used to test how far a dataset
can be used for estimating certain (combinations of) parameters.
After selection of an 'identifiable parameter set' (which has a low
"collinearity") they are fine-tuned by calibration.
2. As the sensitivity analysis is a local analysis (i.e. its outcome depends on the current values of the model parameters) and the fitting routine is used to estimate the best values of the parameters, this is an iterative procedure. This means that identifiable parameters are determined, fitted to the data, then a newly identifiable parameter set is determined, fitted, etcetera until convergenc is reached.
See the paper by Omlin et al. (2001) for more information.
Karline Soetaert <[email protected]>
Brun, R., Reichert, P. and Kunsch, H. R., 2001. Practical Identifiability Analysis of Large Environmental Simulation Models. Water Resour. Res. 37(4): 1015–1030.
Omlin, M., Brun, R. and Reichert, P., 2001. Biogeochemical Model of Lake Zurich: Sensitivity, Identifiability and Uncertainty Analysis. Ecol. Modell. 141: 105–123.
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. doi:10.18637/jss.v033.i03
## ======================================================================= ## Test collinearity values ## ======================================================================= ## linearly related set... => Infinity collin(cbind(1:5, 2*(1:5))) ## unrelated set => 1 MM <- matrix(nr = 4, nc = 2, byrow = TRUE, data = c(-0.400, -0.374, 0.255, 0.797, 0.690, -0.472, -0.546, 0.049)) collin(MM) ## ======================================================================= ## Bacterial 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) ## We wish to estimate parameters gmax and eff by fitting the model to ## these data: Data <- matrix(nc = 2, byrow = TRUE, data = c( 2, 0.14, 4, 0.2, 6, 0.38, 8, 0.42, 10, 0.6, 12, 0.107, 14, 1.3, 16, 2.0, 18, 3.0, 20, 4.5, 22, 6.15, 24, 11, 26, 13.8, 28, 20.0, 30, 31 , 35, 65, 40, 61) ) colnames(Data) <- c("time","Bact") head(Data) Data2 <- matrix(c(2, 100, 20, 93, 30, 55, 50, 0), ncol = 2, byrow = TRUE) colnames(Data2) <- c("time", "Sub") ## Objective function to minimise Objective <- function (x) { # Model cost pars[] <- x out <- solveBact(x) Cost <- modCost(obs = Data2, model = out) # observed data in 2 data.frames return(modCost(obs = Data, model = out, cost = Cost)) } ## 1. Estimate sensitivity functions - all parameters sF <- sensFun(func = Objective, parms = pars, varscale = 1) ## 2. Estimate the collinearity Coll <- collin(sF) ## The larger the collinearity, the less identifiable the data set Coll plot(Coll, log = "y") ## 20 = magical number above which there are identifiability problems abline(h = 20, col = "red") ## select "identifiable" sets with 4 parameters Coll [Coll[,"collinearity"] < 20 & Coll[,"N"]==4,] ## collinearity of one selected parameter set collin(sF, c(1, 3, 5)) collin(sF, 1:5) collin(sF, c("gmax", "eff")) ## collinearity of all combinations of 3 parameters collin(sF, N = 3) ## The collinearity depends on the value of the parameters: P <- pars P[1:2] <- 1 # was: 0.5 collin(sensFun(Objective, P, varscale = 1))
## ======================================================================= ## Test collinearity values ## ======================================================================= ## linearly related set... => Infinity collin(cbind(1:5, 2*(1:5))) ## unrelated set => 1 MM <- matrix(nr = 4, nc = 2, byrow = TRUE, data = c(-0.400, -0.374, 0.255, 0.797, 0.690, -0.472, -0.546, 0.049)) collin(MM) ## ======================================================================= ## Bacterial 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) ## We wish to estimate parameters gmax and eff by fitting the model to ## these data: Data <- matrix(nc = 2, byrow = TRUE, data = c( 2, 0.14, 4, 0.2, 6, 0.38, 8, 0.42, 10, 0.6, 12, 0.107, 14, 1.3, 16, 2.0, 18, 3.0, 20, 4.5, 22, 6.15, 24, 11, 26, 13.8, 28, 20.0, 30, 31 , 35, 65, 40, 61) ) colnames(Data) <- c("time","Bact") head(Data) Data2 <- matrix(c(2, 100, 20, 93, 30, 55, 50, 0), ncol = 2, byrow = TRUE) colnames(Data2) <- c("time", "Sub") ## Objective function to minimise Objective <- function (x) { # Model cost pars[] <- x out <- solveBact(x) Cost <- modCost(obs = Data2, model = out) # observed data in 2 data.frames return(modCost(obs = Data, model = out, cost = Cost)) } ## 1. Estimate sensitivity functions - all parameters sF <- sensFun(func = Objective, parms = pars, varscale = 1) ## 2. Estimate the collinearity Coll <- collin(sF) ## The larger the collinearity, the less identifiable the data set Coll plot(Coll, log = "y") ## 20 = magical number above which there are identifiability problems abline(h = 20, col = "red") ## select "identifiable" sets with 4 parameters Coll [Coll[,"collinearity"] < 20 & Coll[,"N"]==4,] ## collinearity of one selected parameter set collin(sF, c(1, 3, 5)) collin(sF, 1:5) collin(sF, c("gmax", "eff")) ## collinearity of all combinations of 3 parameters collin(sF, N = 3) ## The collinearity depends on the value of the parameters: P <- pars P[1:2] <- 1 # was: 0.5 collin(sensFun(Objective, P, varscale = 1))
Rearranges a data frame in cross tab format by putting all relevant columns below each other, replicating the independent variable and, if necessary, other specified columns. Optionally, an err column is added.
cross2long( data, x, select = NULL, replicate = NULL, error = FALSE, na.rm = FALSE)
cross2long( data, x, select = NULL, replicate = NULL, error = FALSE, na.rm = FALSE)
data |
a data frame (or matrix) with crosstab layout |
x |
name of the independent variable to be replicated |
select |
a vector of column names to be included (see details). All columns are included if not specified. |
replicate |
a vector of names of variables (apart from the independent variable that have to be replicated for every included column (e.g. experimental treatment specification). |
error |
boolean indicating whether the final dataset in long format should contain an extra column for error values (cf. modCost); here filled with 1's. |
na.rm |
whether or not to remove the |
The original data frame is converted from a wide (crosstab) layout (one variable per column) to a long (database) layout (all variable value in one column).
As an example of both formats consider the data, called Dat
consisting
of two observed variables, called "Obs1" and "Obs2", both containing two
observations, at time 1 and 2:
name | time | val | err |
Obs1 | 1 | 50 | 5 |
Obs1 | 2 | 150 | 15 |
Obs2 | 1 | 1 | 0.1 |
Obs2 | 2 | 2 | 0.2 |
for the long format and
time | Obs1 | Obs2 |
1 | 50 | 1 |
2 | 150 | 2 |
for the crosstab format.
The parameters x
, select
, and replicate
should
be disjoint. Although the independent variable always has to be replicated
it should not be given by the replicate
parameter.
A data frame with the following columns:
name |
Column containing the column names of the original crosstab data frame, |
x |
A replication of the independent variable |
y |
The actual data stacked upon each other in one column |
err |
Optional column, filled with NA values (necessary for some other functions) |
... |
all other columns from the original dataset that had to be replicated
(indicated by the parameter |
Tom Van Engeland <[email protected]>
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. doi:10.18637/jss.v033.i03
## ======================================================================= ## Suppose we have measured sediment oxygen concentration profiles ## ======================================================================= depth <- 0:7 O2mud <- c( 6, 1, 0.5, 0.1, 0.05,0, 0, 0) O2silt <- c( 6, 5, 3, 2, 1.5, 1, 0.5, 0) O2sand <- c( 6, 6, 5, 4, 3, 2, 1, 0) zones <- c("a", "b", "b", "c", "c", "d", "d", "e") oxygen <- data.frame(depth = depth, zone = zones, mud = O2mud, silt = O2silt, sand = O2sand ) cross2long(data = oxygen, x = depth, select = c(silt, mud), replicate = zone) cross2long(data = oxygen, x = depth, select = c(mud, -silt), replicate = zone) # twice the same column name: replicates colnames(oxygen)[4] <- "mud" cross2long(data=oxygen, x = depth, select = mud)
## ======================================================================= ## Suppose we have measured sediment oxygen concentration profiles ## ======================================================================= depth <- 0:7 O2mud <- c( 6, 1, 0.5, 0.1, 0.05,0, 0, 0) O2silt <- c( 6, 5, 3, 2, 1.5, 1, 0.5, 0) O2sand <- c( 6, 6, 5, 4, 3, 2, 1, 0) zones <- c("a", "b", "b", "c", "c", "d", "d", "e") oxygen <- data.frame(depth = depth, zone = zones, mud = O2mud, silt = O2silt, sand = O2sand ) cross2long(data = oxygen, x = depth, select = c(silt, mud), replicate = zone) cross2long(data = oxygen, x = depth, select = c(mud, -silt), replicate = zone) # twice the same column name: replicates colnames(oxygen)[4] <- "mud" cross2long(data=oxygen, x = depth, select = mud)
A calibration dataset in database format (cf. modCost for the database format) is extended in order to fit model output using a weighted least squares approach. To this end, the observations are replicated for a certain number of times, and weights are assigned to the replicates according to a Gaussian density function. This density has the relevant observation as mean value. The standard deviation, provided as a parameter, determines the number of inserted replicate observations (see Detail).
This weighted regression approach may be interesting when discontinuities exist in the observational data. Under these circumstances small changes in the timing (or more general the position along the axis of the independent variable) of the model output may have a disproportional impact on the overall goodness-of-fit (e.g. timing of nutrient depletion). Additionally, this approach may be used to model uncertainty in the independent variable (e.g. slices of sediment profiles, or the timing of a sampling).
gaussianWeights (obs, x = x, y = y, xmodel, spread, weight = "none", aggregation = x ,ordering)
gaussianWeights (obs, x = x, y = y, xmodel, spread, weight = "none", aggregation = x ,ordering)
obs |
dataset in long (database) format as is typically used by modCost |
x |
name of the independent variable (typically x, cf. modCost) in
|
y |
name of the dependent variable in |
xmodel |
an ordered vector of unique times at which model output is produced. If not given, the independent variable of the observational dataset is used. |
spread |
standard deviation used to calculate the weights from a
normal density function. This value also determines the number of
points from the model output that are compared to a specific observa-
tion in |
weight |
scaling factor of the modCost function ("sd", "mean", or "none"). The Gaussian weights are multiplied by this factor to account for differences in units. |
aggregation |
vector of column names from the dataset that are used to aggregate observations while calculating the scaling factor. Defaults to the variable name, "name". |
ordering |
Optional extra grouping and ordering of observations. Given as a vector of variable names. If none given, ordering will be done by variable name and independent variable. If both aggregation and ordering variables are given, ordering will be done as follows: x within ordering (in reverse order) within aggregation (in reverse order). Aggregation and ordering should be disjoint sets of variable names. |
Suppose: spread = 1/24 (days; = 1 hour) x = time in days, 1 per hour
Then: obs_i is replicated 7 times (spread = observational periodicity = 1 hour):
=> obs_i-3 = ... = obs_i-1 = obs_i = obs_i+1 = ... = obs_i+3
The weights (W_i+j, for j = -3 ...3) are calculated as follows: W'_i+j = 1/(spread * sqrt(2pi)) * exp(-1/2 * ((obs_i+j - obs_i)/spread)^2
W_i+j = W'_i+j/sum(W_i-3,...,W_i+3) (such that their sum equals 1)
A modified version of obs
is returned with the following extensions:
1. Each observation obs[i] is replicated n times were n represents the number
of modelx
values within the interval [obs_i - (3 * spread), obs_i +
3 * spread)].
2. These replicate observations get the same x
values as their
modeled counterparts (xmodel
).
3. Weights are given in column, called "err"
The returned data frame has the following columns:
"name" or another name specified by the first element of
aggregation
. Usually this column contains the names of the
observed variables.
"x" or another name specified by x
"y" or another name specified by y
"err" containing the calculated weights
The rest of the columns of the data frame given by obs
in
that order.
Tom Van Engeland <[email protected]>
## ======================================================================= ## A Sediment example ## ======================================================================= ## Sediment oxygen concentration is measured every ## centimeter in 3 sediment types depth <- 0:7 observations <- data.frame( profile = rep(c("mud","silt","sand"), each=8), depth = depth, O2 = c(c(6,1,0.5,0.1,0.05,0,0,0), c(6,5,3,2,1.5,1,0.5,0), c(6,6,5,4,3,2,1,0) ) ) ## A model generates profiles with a depth resolution of 1 millimeter modeldepths <- seq(0, 9, by = 0.05) ## All these model outputs are compared with weighed observations. gaussianWeights(obs = observations, x = depth, y = O2, xmodel = modeldepths, spread = 0.1, weight = "none", aggregation = profile) # Weights of one observation in silt at depth 2: Sub <- subset(observations, subset = (profile == "silt" & depth == 2)) plot(Sub[,-1]) SubWW <- gaussianWeights(obs = Sub, x = depth, y = O2, xmodel = modeldepths, spread = 0.5, weight="none", aggregation = profile) SubWW[,-1]
## ======================================================================= ## A Sediment example ## ======================================================================= ## Sediment oxygen concentration is measured every ## centimeter in 3 sediment types depth <- 0:7 observations <- data.frame( profile = rep(c("mud","silt","sand"), each=8), depth = depth, O2 = c(c(6,1,0.5,0.1,0.05,0,0,0), c(6,5,3,2,1.5,1,0.5,0), c(6,6,5,4,3,2,1,0) ) ) ## A model generates profiles with a depth resolution of 1 millimeter modeldepths <- seq(0, 9, by = 0.05) ## All these model outputs are compared with weighed observations. gaussianWeights(obs = observations, x = depth, y = O2, xmodel = modeldepths, spread = 0.1, weight = "none", aggregation = profile) # Weights of one observation in silt at depth 2: Sub <- subset(observations, subset = (profile == "silt" & depth == 2)) plot(Sub[,-1]) SubWW <- gaussianWeights(obs = Sub, x = depth, y = O2, xmodel = modeldepths, spread = 0.5, weight="none", aggregation = profile) SubWW[,-1]
Generates parameter sets arranged on a regular grid.
Grid(parRange, num)
Grid(parRange, num)
parRange |
the range (min, max) of the parameters, a matrix or a data.frame with one row for each parameter, and two columns with the minimum (1st) and maximum (2nd) column. |
num |
the number of random parameter sets to generate. |
The grid design produces the most regular parameter distribution;
there is no randomness involved.
The number of parameter sets generated with Grid
will be <=
num
.
a matrix with one row for each generated parameter set, and one column per parameter.
Karline Soetaert <[email protected]>
Norm
for (multi)normally distributed random parameter sets.
Latinhyper
to generates parameter sets using
latin hypercube sampling.
Unif
for uniformly distributed random parameter sets.
seq
the R-default for generating regular sequences of numbers.
## 4 parameters parRange <- data.frame(min = c(0, 1, 2, 3), max = c(10, 9, 8, 7)) rownames(parRange) <- c("par1", "par2", "par3", "par4") ## grid pairs(Grid(parRange, 500), main = "Grid")
## 4 parameters parRange <- data.frame(min = c(0, 1, 2, 3), max = c(10, 9, 8, 7)) rownames(parRange) <- c("par1", "par2", "par3", "par4") ## grid pairs(Grid(parRange, 500), main = "Grid")
Generates random parameter sets using a latin hypercube sampling algorithm.
Latinhyper(parRange, num)
Latinhyper(parRange, num)
parRange |
the range (min, max) of the parameters, a matrix or a data.frame with one row for each parameter, and two columns with the minimum (1st) and maximum (2nd) column. |
num |
the number of random parameter sets to generate. |
In the latin hypercube sampling, the space for each parameter is
subdivided into num
equally-sized segments and one parameter
value in each of the segments drawn randomly.
a matrix with one row for each generated parameter set, and one column per parameter.
The latin hypercube distributed parameter sets give better coverage in
parameter space than the uniform random design (Unif
).
It is a reasonable choice in case the number of parameter sets is
limited.
Karline Soetaert <[email protected]>
Press, W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery, B. P. (2007) Numerical Recipes in C. Cambridge University Press.
Norm
for (multi)normally distributed random parameter
sets.
Unif
for uniformly distributed random parameter sets.
Grid
to generate random parameter sets arranged on a
regular grid.
## 4 parameters parRange <- data.frame(min = c(0, 1, 2, 3), max = c(10, 9, 8, 7)) rownames(parRange) <- c("par1", "par2", "par3", "par4") ## Latin hypercube pairs(Latinhyper(parRange, 100), main = "Latin hypercube")
## 4 parameters parRange <- data.frame(min = c(0, 1, 2, 3), max = c(10, 9, 8, 7)) rownames(parRange) <- c("par1", "par2", "par3", "par4") ## Latin hypercube pairs(Latinhyper(parRange, 100), main = "Latin hypercube")
Given a solution of a model and observed data, estimates the residuals, and the variable and model costs (sum of squared residuals).
modCost(model, obs, x = "time", y = NULL, err = NULL, weight = "none", scaleVar = FALSE, cost = NULL, ...)
modCost(model, obs, x = "time", y = NULL, err = NULL, weight = "none", scaleVar = FALSE, cost = NULL, ...)
model |
model output, as generated by the integration routine or the steady-state solver, a matrix or a data.frame, with one column per dependent and independent variable. |
obs |
the observed data, either in long (database) format (name, x, y), a data.frame, or in wide (crosstable, or matrix) format - see details. |
x |
the name of the independent variable; it should be a
name occurring both in the |
y |
either |
err |
either |
cost |
if not |
weight |
only if |
scaleVar |
if |
... |
additional arguments passed to R-function |
This function compares model output with observed data.
It computes
the weighted residuals, one for each data point.
the variable costs, i.e. the sum of squared weight residuals per variable.
the model cost, the scaled sum of variable costs .
There are three steps:
1. For any observed data point, i, the weighted residuals are estimated as:
with and where
and
are the modeled, respectively observed value of
data point i.
The weights are equal to 1/error, where the latter can be inputted,
one for each data point by specifying err
as an extra column in
the observed data.
This can only be done when the data input is in long (database) format.
When err
is not inputted, then the weights are specified via argument
weight
which is either:
"none"
, which sets the weight equal to 1 (the default)
"std"
, which sets the weights equal to the reciprocal of the standard deviation
of the observed data (can only be used if there is more than 1 data point)
"mean"
, which uses 1/mean of the absolute value of the
observed data (can only be used if not 0).
2. Then for each observed variable, j, a variable cost is estimated as the sum of squared weighted residuals for this variable:
where is the number of observations for observed
variable j.
3. Finally, the model Cost is estimated as the scaled sum of variable costs:
and where allows to scale the variable
costs relative to the number of observations. This is set by
specifying argument
scaleVar
. If TRUE
, then the variable
costs are rescaled. The default is NOT to rescale
(i.e. =1).
The models typically consist of (a system of) differential equations, which are either solved by:
integration routines, e.g. the routines from package deSolve
,
steady-state estimators, as from package rootSolve
.
The data can be presented in two formats:
data table (long) format; this is a two to four column
data.frame that contains the name
of the observed variable (always
the FIRST column), the (optional) value of the independent variable
(default column name = "time"), the value of the observation
and
the (optional) value of the error
.
For data presented in this format, the names of the column(s) with the
independent variable (x
) and the name of the column that has the
value of the dependent variable y
must be passed to function
modCost
.
crosstable (wide) format; this is a matrix, where each
column denotes one dependent (or independent) variable; the column name
is the name of the observed variable.
When using this format, only the name of the column that contains the
dependent variable must be specified (x
).
As an example of both formats consider the data, called Dat
consisting
of two observed variables, called "Obs1" and "Obs2", both containing two
observations, at time 1 and 2:
name | time | val | err |
Obs1 | 1 | 50 | 5 |
Obs1 | 2 | 150 | 15 |
Obs2 | 1 | 1 | 0.1 |
Obs2 | 2 | 2 | 0.2 |
for the long format and
time | Obs1 | Obs2 |
1 | 50 | 1 |
2 | 150 | 2 |
for the crosstab format. Note, that in the latter case it is not possible to provide separate errors per data point.
By calling modCost several consecutive times (using the cost
argument),
it is possible to combine both types of data files.
a list of type modCost
containing:
model |
one value, the model cost, which equals the sum of scaled variable costs (see details). |
minlogp |
one value, -log(model probablity), where it is assumed
that the data are normally distributed, with standard deviation =
|
var |
the variable costs, a data.frame with, for each observed variable the following (see details):
|
residuals |
the data residual, a data.frame with several columns:
|
In the future, it should be possible to have more than one independent variable present. This is not yet implemented, but it should allow e.g. to fit time series of spatially dependent variables.
Karline Soetaert <[email protected]>
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. doi:10.18637/jss.v033.i03
## ======================================================================= ## Type 1 input: name, time, value ## ======================================================================= ## Create new data: two observed variables, "a", "b" Data <- data.frame(name = c(rep("a", 4), rep("b", 4)), time = c(1:4, 2:5), val = c(runif(4), 1:4)) ## "a nonsense model" Mod <- function (t, y, par) { da <- 0 db <- 1 return(list(c(da, db))) } out <- ode(y = c(a = 0.5, b = 0.5), times = 0:6, func = Mod, parms = NULL) Data # Show out ## The cost function modCost(model = out, obs = Data, y = "val") ## The cost function with a data error added Dat2 <- cbind(Data, Err = Data$val*0.1) # error = 10% of value modCost(model = out, obs = Dat2, y = "val", err = "Err") ## ======================================================================= ## Type 2 input: Matrix format; column names = variable names ## ======================================================================= ## logistic growth model TT <- seq(1, 100, 2.5) N0 <- 0.1 r <- 0.5 K <- 100 ## analytical solution Ana <- cbind(time = TT, N = K/(1 + (K/N0 - 1) * exp(-r*TT))) ## numeric solution logist <- function(t, x, parms) { with(as.list(parms), { dx <- r * x[1] * (1 - x[1]/K) list(dx) }) } time <- 0:100 parms <- c(r = r, K = K) x <- c(N = N0) ## Compare several numerical solutions Euler <- ode(x, time, logist, parms, hini = 2, method = "euler") Rk4 <- ode(x, time, logist, parms, hini = 2, method = "rk4") Lsoda <- ode(x, time, logist, parms) # lsoda is default method Ana2 <- cbind(time = time, N = K/(1 + (K/N0 - 1) * exp(-r * time))) ## the SSR and residuals with respect to the "data" cEuler <- modCost(Euler, Ana)$model cRk4 <- modCost(Rk4 , Ana)$model cLsoda <- modCost(Lsoda, Ana)$model cAna <- modCost(Ana2 , Ana)$model compare <- data.frame(method = c("euler", "rk4", "lsoda", "Ana"), cost = c(cEuler, cRk4, cLsoda, cAna)) ## Plot Euler, RK and analytic solution plot(Euler, Rk4, col = c("red", "blue"), obs = Ana, main = "logistic growth", xlab = "time", ylab = "N") legend("bottomright", c("exact", "euler", "rk4"), pch = c(1, NA, NA), col = c("black", "red", "blue"), lty = c(NA, 1, 2)) legend("right", ncol = 2, title = "SSR", legend = c(as.character(compare[,1]), format(compare[,2], digits = 2))) compare ## ======================================================================= ## Now suppose we do not know K and r and they are to be fitted... ## The "observations" are the analytical solution ## ======================================================================= ## Run the model with initial guess: K = 10, r = 2 parms["K"] <- 10 parms["r"] <- 2 init <- ode(x, time, logist, parms) ## FITTING algorithm uses modFit ## First define the objective function (model cost) to be minimised ## more general: using modFit Cost <- function(P) { parms["K"] <- P[1] parms["r"] <- P[2] out <- ode(x, time, logist, parms) return(modCost(out, Ana)) } (Fit<-modFit(p = c(K = 10, r = 2), f = Cost)) summary(Fit) ## run model with the optimized value: parms[c("K", "r")] <- Fit$par fitted <- ode(x, time, logist, parms) ## show results, compared with "observations" plot(init, fitted, col = c("green", "blue"), lwd = 2, lty = 1, obs = Ana, obspar = list(col = "red", pch = 16, cex = 2), main = "logistic growth", xlab = "time", ylab = "N") legend("right", c("initial", "fitted"), col = c("green", "blue"), lwd = 2) Cost(Fit$par)
## ======================================================================= ## Type 1 input: name, time, value ## ======================================================================= ## Create new data: two observed variables, "a", "b" Data <- data.frame(name = c(rep("a", 4), rep("b", 4)), time = c(1:4, 2:5), val = c(runif(4), 1:4)) ## "a nonsense model" Mod <- function (t, y, par) { da <- 0 db <- 1 return(list(c(da, db))) } out <- ode(y = c(a = 0.5, b = 0.5), times = 0:6, func = Mod, parms = NULL) Data # Show out ## The cost function modCost(model = out, obs = Data, y = "val") ## The cost function with a data error added Dat2 <- cbind(Data, Err = Data$val*0.1) # error = 10% of value modCost(model = out, obs = Dat2, y = "val", err = "Err") ## ======================================================================= ## Type 2 input: Matrix format; column names = variable names ## ======================================================================= ## logistic growth model TT <- seq(1, 100, 2.5) N0 <- 0.1 r <- 0.5 K <- 100 ## analytical solution Ana <- cbind(time = TT, N = K/(1 + (K/N0 - 1) * exp(-r*TT))) ## numeric solution logist <- function(t, x, parms) { with(as.list(parms), { dx <- r * x[1] * (1 - x[1]/K) list(dx) }) } time <- 0:100 parms <- c(r = r, K = K) x <- c(N = N0) ## Compare several numerical solutions Euler <- ode(x, time, logist, parms, hini = 2, method = "euler") Rk4 <- ode(x, time, logist, parms, hini = 2, method = "rk4") Lsoda <- ode(x, time, logist, parms) # lsoda is default method Ana2 <- cbind(time = time, N = K/(1 + (K/N0 - 1) * exp(-r * time))) ## the SSR and residuals with respect to the "data" cEuler <- modCost(Euler, Ana)$model cRk4 <- modCost(Rk4 , Ana)$model cLsoda <- modCost(Lsoda, Ana)$model cAna <- modCost(Ana2 , Ana)$model compare <- data.frame(method = c("euler", "rk4", "lsoda", "Ana"), cost = c(cEuler, cRk4, cLsoda, cAna)) ## Plot Euler, RK and analytic solution plot(Euler, Rk4, col = c("red", "blue"), obs = Ana, main = "logistic growth", xlab = "time", ylab = "N") legend("bottomright", c("exact", "euler", "rk4"), pch = c(1, NA, NA), col = c("black", "red", "blue"), lty = c(NA, 1, 2)) legend("right", ncol = 2, title = "SSR", legend = c(as.character(compare[,1]), format(compare[,2], digits = 2))) compare ## ======================================================================= ## Now suppose we do not know K and r and they are to be fitted... ## The "observations" are the analytical solution ## ======================================================================= ## Run the model with initial guess: K = 10, r = 2 parms["K"] <- 10 parms["r"] <- 2 init <- ode(x, time, logist, parms) ## FITTING algorithm uses modFit ## First define the objective function (model cost) to be minimised ## more general: using modFit Cost <- function(P) { parms["K"] <- P[1] parms["r"] <- P[2] out <- ode(x, time, logist, parms) return(modCost(out, Ana)) } (Fit<-modFit(p = c(K = 10, r = 2), f = Cost)) summary(Fit) ## run model with the optimized value: parms[c("K", "r")] <- Fit$par fitted <- ode(x, time, logist, parms) ## show results, compared with "observations" plot(init, fitted, col = c("green", "blue"), lwd = 2, lty = 1, obs = Ana, obspar = list(col = "red", pch = 16, cex = 2), main = "logistic growth", xlab = "time", ylab = "N") legend("right", c("initial", "fitted"), col = c("green", "blue"), lwd = 2) Cost(Fit$par)
Given a model consisting of differential equations, estimates the global effect of certain (sensitivity) parameters on 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 is useful for “what-if” scenarios.
If the output variables consist of a time-series or spatially dependent, use sensRange instead.
modCRL(func, parms = NULL, sensvar = NULL, dist = "unif", parInput = NULL, parRange = NULL, parMean = NULL, parCovar = NULL, num = 100, ...) ## S3 method for class 'modCRL' summary(object, ...) ## S3 method for class 'modCRL' plot(x, which = NULL, trace = FALSE, ask = NULL, ...) ## S3 method for class 'modCRL' pairs(x, which = 1:ncol(x), nsample = NULL, ...) ## S3 method for class 'modCRL' hist(x, which = 1:ncol(x), ask = NULL, ...)
modCRL(func, parms = NULL, sensvar = NULL, dist = "unif", parInput = NULL, parRange = NULL, parMean = NULL, parCovar = NULL, num = 100, ...) ## S3 method for class 'modCRL' summary(object, ...) ## S3 method for class 'modCRL' plot(x, which = NULL, trace = FALSE, ask = NULL, ...) ## S3 method for class 'modCRL' pairs(x, which = 1:ncol(x), nsample = NULL, ...) ## S3 method for class 'modCRL' hist(x, which = 1:ncol(x), ask = NULL, ...)
func |
an R-function that has as first argument |
parms |
parameters passed to |
sensvar |
the output variables for which the sensitivity needs to be
estimated. Either |
dist |
the distribution according to which the parameters should be
generated, one of The input parameters for the distribution are specified by
|
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 |
parInput |
a matrix with dimension (*, npar) with the values of the sensitivity parameters. |
parMean |
only when |
parCovar |
only when |
num |
the number of times the model has to be run. Set large enough.
If |
object |
an object of class |
x |
an object of class |
which |
the name or the index to the variables and parameters that should be plotted. Default = all variables and parameters. |
nsample |
the number of xy pairs to be plotted on the upper
panel in the pairs plot. When |
trace |
if |
ask |
logical; if |
... |
additional arguments passed to function |
a data.frame of type modCRL
containing the parameter(s) and
the corresponding values of the sensitivity output variables.
The list returned by modCRL
has a method for the generic functions
summary
and plot
– see note.
The following methods are included:
summary, estimates summary statistics for the
sensitivity variables, a table with as many rows as there are
variables (or elements in the vector 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 plot of the modCRL
output,
either one plot for each sensitivity variable and with the
parameter value on the x-axis. This only works when there is only
one parameter!
OR
one plot for each parameter value on the x-axis. This only works when there is only one variable!
hist, produces a histogram of the modCRL
output
parameters and variables.
pairs, produces a pairs plot of the modCRL
output.
The data.frame of type modCRL
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 - see
notes in help of sensRange
.
Karline Soetaert <[email protected]>.
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. doi:10.18637/jss.v033.i03
## ======================================================================= ## 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, main = "Bacteria", xlab = "time, hour", ylab = "molC/m3", type = "l", lwd = 2) ## Function that returns the last value of the simulation SF <- function (p) { pars["eff"] <- p out <- solveBact(pars) return(out[nrow(out), 2:3]) } parRange <- matrix(nr = 1, nc = 2, c(0.2, 0.8), dimnames = list("eff", c("min", "max"))) parRange CRL <- modCRL(func = SF, parRange = parRange) plot(CRL) # plots both variables plot(CRL, which = c("eff", "Bact"), trace = FALSE) #selects one
## ======================================================================= ## 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, main = "Bacteria", xlab = "time, hour", ylab = "molC/m3", type = "l", lwd = 2) ## Function that returns the last value of the simulation SF <- function (p) { pars["eff"] <- p out <- solveBact(pars) return(out[nrow(out), 2:3]) } parRange <- matrix(nr = 1, nc = 2, c(0.2, 0.8), dimnames = list("eff", c("min", "max"))) parRange CRL <- modCRL(func = SF, parRange = parRange) plot(CRL) # plots both variables plot(CRL, which = c("eff", "Bact"), trace = FALSE) #selects one
Fitting a model to data, with lower and/or upper bounds
modFit(f, p, ..., lower = -Inf, upper = Inf, method = c("Marq", "Port", "Newton", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Pseudo", "bobyqa"), jac = NULL, control = list(), hessian = TRUE) ## S3 method for class 'modFit' summary(object, cov=TRUE,...) ## S3 method for class 'modFit' deviance(object, ...) ## S3 method for class 'modFit' coef(object, ...) ## S3 method for class 'modFit' residuals(object, ...) ## S3 method for class 'modFit' df.residual(object, ...) ## S3 method for class 'modFit' plot(x, ask = NULL, ...) ## S3 method for class 'summary.modFit' print(x, digits = max(3, getOption("digits") - 3), ...)
modFit(f, p, ..., lower = -Inf, upper = Inf, method = c("Marq", "Port", "Newton", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Pseudo", "bobyqa"), jac = NULL, control = list(), hessian = TRUE) ## S3 method for class 'modFit' summary(object, cov=TRUE,...) ## S3 method for class 'modFit' deviance(object, ...) ## S3 method for class 'modFit' coef(object, ...) ## S3 method for class 'modFit' residuals(object, ...) ## S3 method for class 'modFit' df.residual(object, ...) ## S3 method for class 'modFit' plot(x, ask = NULL, ...) ## S3 method for class 'summary.modFit' print(x, digits = max(3, getOption("digits") - 3), ...)
f |
a function to be minimized, with first argument the vector
of parameters over which minimization is to take place. It should
return either a vector of residuals (of model versus data) or
an element of class modCost (as returned by a call to
|
p |
initial values for the parameters to be optimized over. |
... |
additional arguments passed to function |
lower |
lower bounds on the parameters; if unbounded set equal to
|
upper |
upper bounds on the parameters; if unbounded set equal to
|
method |
The method to be used, one of "Marq", "Port", "Newton", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Pseudo", "bobyqa" - see details. |
jac |
A function that calculates the Jacobian; it should be
called as |
hessian |
|
control |
additional control arguments passed to the optimisation
routine - see details of |
object |
an object of class |
x |
an object of class |
digits |
number of significant digits in printout. |
cov |
when |
ask |
logical; if |
Note that arguments after ...
must be matched exactly.
The method to be used is specified by argument method
which can
be one of the methods from function optim
:
"Nelder-Mead", the default from optim
,
"BFGS", a quasi-Newton method,
"CG", a conjugate-gradient method,
"L-BFGS-B", constrained quasi-Newton method,
"SANN", method of simulated annealing.
Or one of the following:
"Marq", the Levenberg-Marquardt algorithm
(nls.lm
from package minpack
) -
the default. Note that this method is the only least squares method.
"Newton", a Newton-type algorithm (see nlm
),
"Port", the Port algorithm (see nlminb
),
"Pseudo", a pseudorandom-search algorithm
(see pseudoOptim
),
"bobyqa", derivative-free optimization by quadratic
approximation from package minqa
.
For difficult problems it may be efficient to perform some iterations
with Pseudo
, which will bring the algorithm near the vicinity
of a (the) minimum, after which the default algorithm (Marq
) is
used to locate the minimum more precisely.
The implementation for the routines from optim
differs
from constrOptim
which implements an adaptive barrier
algorithm and which allows a more flexible implementation of linear
constraints.
For all methods except L-BFGS-B
, Port
, Pseudo
, and
bobyqa
that handle box constraints internally, bounds on parameters are
imposed by a transformation of the parameters to be fitted.
In case both lower and upper bounds are specified, this is achieved by a tangens and arc tangens transformation.
This is, parameter values, , generated by the optimisation
routine, and which are located in the range [-Inf, Inf] are
transformed, before they are passed to
f
as:
.
which maps them into the interval [lower, upper].
Before the optimisation routine is called, the original parameter values, as
given by argument p
are mapped from [lower,upper] to [-Inf, Inf] by:
In case only lower or upper bounds are specified, this is achieved by a log transformation and a corresponding exponential back transformation.
In case parameters are transformed (all methods) or in case the
method
Port
, Pseudo
, Marq
or bobyqa
is selected,
the Hessian is approximated as ,
where J is the Jacobian, estimated by finite differences.
This ignores the second derivative terms, but this is reasonable if the method has truly converged to the minimum.
Note that finite differences are not extremely precise.
In case the Levenberg-Marquard method (Marq
) is used, and parameters
are not transformed, 0.5 times the Hessian of the least squares problem is
returned by nls.lm
, the original Marquardt algorithm. To make it
compatible, this value is multiplied with 2 and the TRUE Hessian is thus
returned by modFit
.
a list of class modFit containing the results as returned from the called optimisation routines.
This includes the following:
par |
the best set of parameters found. |
ssr |
the sum of squared residuals, evaluated at the best set of parameters. |
Hessian |
A symmetric matrix giving an estimate of the Hessian at the solution found - see note. |
residuals |
the result of the last |
ms |
the mean squared residuals, i.e. |
var_ms |
the weighted and scaled variable mean squared residuals,
one value per observed variable; only when |
var_ms_unscaled |
the weighted, but not scaled variable mean squared residuals |
var_ms_unweighted |
the raw variable mean squared residuals, unscaled and unweighted. |
... |
any other arguments returned by the called optimisation routine. |
Note: this means that some return arguments of the original optimisation functions are renamed.
More specifically, "objective" and "counts" from routine
nlminb
(method = "Port") are renamed; "value" and
"counts"; "niter" and "minimum" from routine
nls.lm
(method=Marq) are renamed; "counts"
and "value"; "minimum" and "estimate" from routine nlm
(method = "Newton"
) are renamed.
The list returned by modFit
has a method for the
summary
, deviance
, coef
,
residuals
, df.residual
and
print.summary
– see note.
The summary
method is based on an estimate of the
parameter covariance matrix.
In computing the covariance matrix of the fitted parameters, the problem is
treated as if it were a linear least squares problem, linearizing around
the parameter values that minimize .
The covariance matrix is estimated as .
This computation relies on several things, i.e.:
the parameter values are located at the minimum (i.e. the fitting algorithm has converged).
the observations are subject to independent errors whose
variances are well estimated by
times the residual sum
of squares (where n = number of data points, p = number of parameters).
the model is not too nonlinear.
This means that the estimated covariance (correlation) matrix and the confidence intervals derived from it may be worthless if the assumptions behind the covariance computation are invalid.
If in doubt about the validity of the summary computations, use Monte Carlo
fitting instead, or run a modMCMC
.
Other methods included are:
deviance
, which returns the model deviance,
coef
, which extracts the values of the fitted parameters,
residuals
,which extracts the model residuals,
df.residual
which returns the residual degrees of freedom
print.summary
, producing a nice printout of the summary.
Specifying a function to estimate the Jacobian matrix via argument
jac
may increase speed. The Jacobian is used in the methods
"Marq", "BFGS", "CG", "L-BFGS", "Port", and is also used at the end,
to estimate the Hessian at the optimum.
Specification of the gradient
in routines "BFGS", "CG", "L-BFGS"
from optim
and "port" from nlminb
is not allowed here.
Within modFit
, the gradient is rather estimated from the Jacobian
jac
and the function f
.
Karline Soetaert <[email protected]>,
Thomas Petzoldt <[email protected]>
Bates, D., Mullen, K. D. Nash, J. C. and Varadhan, R. 2014. minqa: Derivative-free optimization algorithms by quadratic approximation. R package. https://cran.r-project.org/package=minqa
Gay, D. M., 1990. Usage Summary for Selected Optimization Routines. Computing Science Technical Report No. 153. AT&T Bell Laboratories, Murray Hill, NJ 07974.
Powell, M. J. D. (2009). The BOBYQA algorithm for bound constrained optimization without derivatives. Report No. DAMTP 2009/NA06, Centre for Mathematical Sciences, University of Cambridge, UK. https://www.damtp.cam.ac.uk/user/na/NA_papers/NA2009_06.pdf
Press, W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery, B. P., 2007. Numerical Recipes in C. Cambridge University Press.
Price, W.L., 1977. A Controlled Random Search Procedure for Global Optimisation. The Computer Journal, 20: 367-370. doi:10.1093/comjnl/20.4.367
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. doi:10.18637/jss.v033.i03
Please see also additional publications on the help pages of the individual algorithms.
constrOptim
for constrained optimization.
## ======================================================================= ## logistic growth model ## ======================================================================= TT <- seq(1, 60, 5) N0 <- 0.1 r <- 0.5 K <- 100 ## perturbed analytical solution Data <- data.frame( time = TT, N = K / (1+(K/N0-1) * exp(-r*TT)) * (1 + rnorm(length(TT), sd = 0.01)) ) plot(TT, Data[,"N"], ylim = c(0, 120), pch = 16, col = "red", main = "logistic growth", xlab = "time", ylab = "N") ##=================================== ## Fitted with analytical solution # ##=================================== ## initial "guess" parms <- c(r = 2, K = 10, N0 = 5) ## analytical solution model <- function(parms,time) with (as.list(parms), return(K/(1+(K/N0-1)*exp(-r*time)))) ## run the model with initial guess and plot results lines(TT, model(parms, TT), lwd = 2, col = "green") ## FITTING algorithm 1 ModelCost <- function(P) { out <- model(P, TT) return(Data$N-out) # residuals } (Fita <- modFit(f = ModelCost, p = parms)) times <- 0:60 lines(times, model(Fita$par, times), lwd = 2, col = "blue") summary(Fita) ##=================================== ## Fitted with numerical solution # ##=================================== ## numeric solution logist <- function(t, x, parms) { with(as.list(parms), { dx <- r * x[1] * (1 - x[1]/K) list(dx) }) } ## model cost, ModelCost2 <- function(P) { out <- ode(y = c(N = P[["N0"]]), func = logist, parms = P, times = c(0, TT)) return(modCost(out, Data)) # object of class modCost } Fit <- modFit(f = ModelCost2, p = parms, lower = rep(0, 3), upper = c(5, 150, 10)) out <- ode(y = c(N = Fit$par[["N0"]]), func = logist, parms = Fit$par, times = times) lines(out, col = "red", lty = 2) legend("right", c("data", "original", "fitted analytical", "fitted numerical"), lty = c(NA, 1, 1, 2), lwd = c(NA, 2, 2, 1), col = c("red", "green", "blue", "red"), pch = c(16, NA, NA, NA)) summary(Fit) plot(residuals(Fit)) ## ======================================================================= ## the use of the Jacobian ## ======================================================================= ## We use modFit to solve a set of linear equations A <- matrix(nr = 30, nc = 30, runif(900)) X <- runif(30) B <- A %*% X ## try to find vector "X"; the Jacobian is matrix A ## Function that returns the vector of residuals FUN <- function(x) as.vector(A %*% x - B) ## Function that returns the Jacobian JAC <- function(x) A ## The port algorithm print(system.time( mf <- modFit(f = FUN, p = runif(30), method = "Port") )) print(system.time( mf <- modFit(f = FUN, p = runif(30), method = "Port", jac = JAC) )) max(abs(mf$par - X)) # should be very small ## BFGS print(system.time( mf <- modFit(f = FUN, p = runif(30), method = "BFGS") )) print(system.time( mf <- modFit(f = FUN, p = runif(30), method = "BFGS", jac = JAC) )) max(abs(mf$par - X)) ## Levenberg-Marquardt print(system.time( mf <- modFit(f = FUN, p = runif(30), jac = JAC) )) max(abs(mf$par - X))
## ======================================================================= ## logistic growth model ## ======================================================================= TT <- seq(1, 60, 5) N0 <- 0.1 r <- 0.5 K <- 100 ## perturbed analytical solution Data <- data.frame( time = TT, N = K / (1+(K/N0-1) * exp(-r*TT)) * (1 + rnorm(length(TT), sd = 0.01)) ) plot(TT, Data[,"N"], ylim = c(0, 120), pch = 16, col = "red", main = "logistic growth", xlab = "time", ylab = "N") ##=================================== ## Fitted with analytical solution # ##=================================== ## initial "guess" parms <- c(r = 2, K = 10, N0 = 5) ## analytical solution model <- function(parms,time) with (as.list(parms), return(K/(1+(K/N0-1)*exp(-r*time)))) ## run the model with initial guess and plot results lines(TT, model(parms, TT), lwd = 2, col = "green") ## FITTING algorithm 1 ModelCost <- function(P) { out <- model(P, TT) return(Data$N-out) # residuals } (Fita <- modFit(f = ModelCost, p = parms)) times <- 0:60 lines(times, model(Fita$par, times), lwd = 2, col = "blue") summary(Fita) ##=================================== ## Fitted with numerical solution # ##=================================== ## numeric solution logist <- function(t, x, parms) { with(as.list(parms), { dx <- r * x[1] * (1 - x[1]/K) list(dx) }) } ## model cost, ModelCost2 <- function(P) { out <- ode(y = c(N = P[["N0"]]), func = logist, parms = P, times = c(0, TT)) return(modCost(out, Data)) # object of class modCost } Fit <- modFit(f = ModelCost2, p = parms, lower = rep(0, 3), upper = c(5, 150, 10)) out <- ode(y = c(N = Fit$par[["N0"]]), func = logist, parms = Fit$par, times = times) lines(out, col = "red", lty = 2) legend("right", c("data", "original", "fitted analytical", "fitted numerical"), lty = c(NA, 1, 1, 2), lwd = c(NA, 2, 2, 1), col = c("red", "green", "blue", "red"), pch = c(16, NA, NA, NA)) summary(Fit) plot(residuals(Fit)) ## ======================================================================= ## the use of the Jacobian ## ======================================================================= ## We use modFit to solve a set of linear equations A <- matrix(nr = 30, nc = 30, runif(900)) X <- runif(30) B <- A %*% X ## try to find vector "X"; the Jacobian is matrix A ## Function that returns the vector of residuals FUN <- function(x) as.vector(A %*% x - B) ## Function that returns the Jacobian JAC <- function(x) A ## The port algorithm print(system.time( mf <- modFit(f = FUN, p = runif(30), method = "Port") )) print(system.time( mf <- modFit(f = FUN, p = runif(30), method = "Port", jac = JAC) )) max(abs(mf$par - X)) # should be very small ## BFGS print(system.time( mf <- modFit(f = FUN, p = runif(30), method = "BFGS") )) print(system.time( mf <- modFit(f = FUN, p = runif(30), method = "BFGS", jac = JAC) )) max(abs(mf$par - X)) ## Levenberg-Marquardt print(system.time( mf <- modFit(f = FUN, p = runif(30), jac = JAC) )) max(abs(mf$par - X))
Performs a Markov Chain Monte Carlo simulation, using an adaptive Metropolis (AM) algorithm and including a delayed rejection (DR) procedure.
modMCMC(f, p, ..., jump = NULL, lower = -Inf, upper = +Inf, prior = NULL, var0 = NULL, wvar0 = NULL, n0 = NULL, niter = 1000, outputlength = niter, burninlength = 0, updatecov = niter, covscale = 2.4^2/length(p), ntrydr = 1, drscale = NULL, verbose = TRUE) ## S3 method for class 'modMCMC' summary(object, remove = NULL, ...) ## S3 method for class 'modMCMC' pairs(x, Full = FALSE, which = 1:ncol(x$pars), remove = NULL, nsample = NULL, ...) ## S3 method for class 'modMCMC' hist(x, Full = FALSE, which = 1:ncol(x$pars), remove = NULL, ask = NULL, ...) ## S3 method for class 'modMCMC' plot(x, Full = FALSE, which = 1:ncol(x$pars), trace = TRUE, remove = NULL, ask = NULL, ...)
modMCMC(f, p, ..., jump = NULL, lower = -Inf, upper = +Inf, prior = NULL, var0 = NULL, wvar0 = NULL, n0 = NULL, niter = 1000, outputlength = niter, burninlength = 0, updatecov = niter, covscale = 2.4^2/length(p), ntrydr = 1, drscale = NULL, verbose = TRUE) ## S3 method for class 'modMCMC' summary(object, remove = NULL, ...) ## S3 method for class 'modMCMC' pairs(x, Full = FALSE, which = 1:ncol(x$pars), remove = NULL, nsample = NULL, ...) ## S3 method for class 'modMCMC' hist(x, Full = FALSE, which = 1:ncol(x$pars), remove = NULL, ask = NULL, ...) ## S3 method for class 'modMCMC' plot(x, Full = FALSE, which = 1:ncol(x$pars), trace = TRUE, remove = NULL, ask = NULL, ...)
f |
the function to be evaluated, with first argument the vector
of parameters which should be varied. It should return either the
model residuals, an element of class modCost (as returned by
a call to |
p |
initial values for the parameters to be optimized over. |
... |
additional arguments passed to function |
jump |
jump length, either a number, a vector with length equal to the total number of parameters, a covariance matrix, or a function that takes as input the current values of the parameters and produces as output the perturbed parameters. See details. |
prior |
-2*log(parameter prior probability), either a function
that is called as |
var0 |
initial model variance; if |
wvar0 |
"weight" for the initial model variance – see details. |
n0 |
parameter used for weighing the initial model variance -
if |
lower |
lower bounds on the parameters; for unbounded parameters
set equal to |
upper |
upper bounds on the parameters; for unbounded parameters
set equal to |
niter |
number of iterations for the MCMC. |
outputlength |
number of iterations kept in the output; should
be smaller or equal to |
updatecov |
number of iterations after which the parameter covariance matrix is (re)evaluated based on the parameters kept thus far, and used to update the MCMC jumps. |
covscale |
scale factor for the parameter covariance matrix, used to perform the MCMC jumps. |
burninlength |
number of initial iterations to be removed from output. |
ntrydr |
maximal number of tries for the delayed rejection procedure. It is generally not a good idea to set this to a too large value. |
drscale |
for each try during delayed rejection, the cholesky
decomposition of the proposal matrix is scaled with this amount; if
|
verbose |
if |
object |
an object of class |
x |
an object of class |
Full |
If TRUE then not only the parameters will be plotted, but also the function value and (if appropriate) the model variance(s). |
which |
the name or the index to the parameters that should be
plotted. Default = all parameters. If |
trace |
if |
remove |
a list with indices of the runs that should be removed (e.g. to remove runs during burnin). |
nsample |
the number of xy pairs to be plotted on the upper
panel in the pairs plot. When |
ask |
logical; if |
Note that arguments after ... must be matched exactly.
R-function f
is called as f(p, ...)
. It should return
either -2 times the log likelihood of the model (one value), the
residuals between model and data or an item of class modFit
(as
created by function modFit
.
In the latter two cases, it is assumed that the prior distribution for
is either non-informative or gaussian. If gaussian, it
can be treated as a sum of squares (SS). If the measurement function
is defined as:
where is the measurement error, assumed normally
distribution, then the posterior for the parameters will be estimated
as:
and where is the error variance, SS is the sum of squares
function
.
If non-informative priors are used, then
.
The error variance is considered a nuisance parameter.
A prior distribution of it should be specified and a posterior distribution
is estimated.
If wvar0
or n0
is >0, then the variances are sampled as
conjugate priors from the inverse gamma distribution with parameters
var0
and n0=wvar0*n
. Larger values of wvar0
keep
these samples closer to var0
.
Thus, at each step, 1/ the error variance () is sampled
from a gamma distribution:
where n
is the number of data points and where
, and where the second argument to the gamma function
is the shape parameter.
The prior parameters (var0
and wvar0
) are the prior mean
for and the prior accuracy.
By setting wvar0
equal to 1, equal weight is given to the
prior and the current value.
If wvar0
is 0 then the prior is ignored.
If wvar0
is NULL
(the default) then the error variances are
assumed to be fixed.
var0
estimates the variance of the measured components. In case
independent estimates are not available, these variances can be
obtained from the mean squares of fitted residuals. (e.g. as reported
in modFit
). See the examples. (but note that this is not truly
independent information)
var0
is either one value, or a value for each observed
variable, or a value for each observed data point.
When var0
is not NULL
, then f
is assumed to
return the model residuals OR an instance of class modCost
.
When var0=NULL
, then f
should return either
-2*log(probability of the model), or an instance of class
modCost
.
modMCMC
implements the Metropolis-Hastings method. The proposal
distribution, which is used to generate new parameter values is the
(multidimensional) Gaussian density distribution, with standard
deviation given by jump
.
jump
can be either one value, a vector of length = number of
parameters or a parameter covariance matrix (nrow = ncol = number
parameters).
The jump parameter, jump
thus determines how much the new
parameter set will deviate from the old one.
If jump
is one value, or a vector, then the new parameter
values are generated by sampling a normal distribution with standard
deviation equal to jump
. A larger value will lead to larger
jumps in the parameter space, but acceptance of new points can get
very low. Smaller jump lengths increase the acceptance rate, but the
algorithm may move too slowly, and too many runs may be needed to scan
the parameter space.
If jump
is NULL
, then the jump length is taken as 10%
of the parameter value as given in p
.
jump
can also be a proposal covariance matrix. In this case,
the new parameter values are generated by sampling a multidimensional
normal distribution. It can be efficient to initialise jump
using the parameter covariance as resulting from fitting the model
(e.g. using modFit
) – see examples.
Finally, jump
can also be an R-function that takes as input the
current values of the parameters and returns the new parameter values.
Two methods are implemented to increase the number of accepted runs.
In the "adaptive Metropolis" method, new parameters are generated with a covariance matrix that is estimated from the parameters generated (and saved) thus far. The idea behind this is that the MCMC method is more efficient if the proposal covariance (to generate new parameter values) is somehow tuned to the shape and size of the target distribution.
Setting updatecov
smaller than niter
will trigger
this functionality. In this case, every updatecov
iterations, the jump covariance matrix will be estimated from the
covariance matrix of the saved parameter values. The covariance
matrix is scaled with where npar is the number
of parameters, unless
covscale
has been given a different
value. Thus, where the small number
is added on the diagonal of the covariance matrix to prevent it
from becoming singular.
Note that a problem of adapting the proposal distribution using the MCMC results so far is that standard convergence results do not apply. One solution is to use adaptation only for the burn-in period and discard the part of the chain where adaptation has been used.
Thus, when using updatecov
with a positive value of
burninlength
, the proposal distribution is only updated
during burnin. If burninlength
= 0 though, the updates
occur throughout the entire simulation.
When using the adaptive Metropolis method, it is best to start with a small value of the jump length.
In the "delayed rejection" method, new parameter
values are tried upon rejection. The process of delaying rejection
can be iterated for at most ntrydr
trials. Setting
ntrydr
equal to 1 (the default) toggles off delayed
rejection.
During the delayed rejection procedure, new parameters are
generated from the last accepted value by scaling the jump
covariance matrix with a factor as specified in
drscale
. The acceptance probability of this new set depends
on the candidates so far proposed and rejected, in such a way that
reversibility of the Markov chain is preserved. See Haario et
al. (2005, 2006) for more details.
Convergence of the MCMC chain can be checked via plot
, which
plots for each iteration the values of all parameters, and if
Full
is TRUE
, of the function value (SS) and (if
appropriate) the modeled variance. If converged, there should be no
visible drift.
In addition, the methods from package coda
become available by
making the object returned by modMCMC
of class mcmc
, as
used in the methods of coda
. For instance, if object
MCMCres
is returned by modMCMC
then
as.mcmc(MCMCres$pars)
will make an instance of class
mcmc
, usable by coda
.
The burninlength
is the number of initial steps that is not
included in the output. It can be useful if the initial value of the
parameters is far from the optimal value. Starting the MCMC with the
best fit parameter set will alleviate the need for using
burninlength
.
a list of class modMCMC containing the results as returned from the Markov chain.
This includes the following:
pars |
an array with dimension ( |
SS |
vector with the sum of squares function, one for each row
in |
naccepted |
the number of accepted runs. |
sig |
the sampled error variance |
bestpar |
the parameter set that gave the highest probability. |
bestfunp |
the function value corresponding to |
prior |
the parameter prior, one value for each row in
|
count |
information about the MCMC chain: number of delayed
rejection steps ( |
settings |
the settings for error covariance calculation,
i.e. arguments |
The list returned by modMCMC
has methods for the generic
functions summary
, plot
,
pairs
– see note.
The following S3 methods are provided:
summary, produces summary statistics of the MCMC results
plot, plots the MCMC results, for all parameters. Use it to check convergence.
pairs, produces a pairs plot of the MCMC results; overrides
the default gap = 0
, upper.panel = NA
, and
diag.panel
.
It is also possible to use the methods from the coda
package,
e.g. densplot
.
To do that, first the modMCMC
object has to be converted to an
mcmc object. See the examples for an application.
Karline Soetaert <[email protected]>
Marko Laine <[email protected]>
Laine, M., 2008. Adaptive MCMC Methods With Applications in Environmental and Geophysical Models. Finnish Meteorological Institute contributions 69, ISBN 978-951-697-662-7, Finnish Meteorological Institute, Helsinki.
Haario, H., Saksman, E. and Tamminen, J., 2001. An Adaptive Metropolis Algorithm. Bernoulli 7, pp. 223–242. doi:10.2307/3318737
Haario, H., Laine, M., Mira, A. and Saksman, E., 2006. DRAM: Efficient Adaptive MCMC. Statistics and Computing, 16(4), 339–354. doi:10.1007/s11222-006-9438-0
Haario, H., Saksman, E. and Tamminen, J., 2005. Componentwise Adaptation for High Dimensional MCMC. Computational Statistics 20(2), 265–274. doi:10.1007/BF02789703
Gelman, A. Varlin, J. B., Stern, H. S. and Rubin, D. B., 2004. Bayesian Data Analysis. Second edition. Chapman and Hall, Boca Raton.
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. doi:10.18637/jss.v033.i03
modFit
for constrained model fitting
## ======================================================================= ## Sampling a 3-dimensional normal distribution, ## ======================================================================= # mean = 1:3, sd = 0.1 # f returns -2*log(probability) of the parameter values NN <- function(p) { mu <- c(1,2,3) -2*sum(log(dnorm(p, mean = mu, sd = 0.1))) #-2*log(probability) } # simple Metropolis-Hastings MCMC <- modMCMC(f = NN, p = 0:2, niter = 5000, outputlength = 1000, jump = 0.5) # More accepted values by updating the jump covariance matrix... MCMC <- modMCMC(f = NN, p = 0:2, niter = 5000, updatecov = 100, outputlength = 1000, jump = 0.5) summary(MCMC) plot(MCMC) # check convergence pairs(MCMC) ## ======================================================================= ## test 2: sampling a 3-D normal distribution, larger standard deviation... ## noninformative priors, lower and upper bounds imposed on parameters ## ======================================================================= NN <- function(p) { mu <- c(1,2,2.5) -2*sum(log(dnorm(p, mean = mu, sd = 0.5))) #-2*log(probability) } MCMC2 <- modMCMC(f = NN, p = 1:3, niter = 2000, burninlength = 500, updatecov = 10, jump = 0.5, lower = c(0, 2, 1), upper = c(1, 3, 3)) plot(MCMC2) hist(MCMC2, breaks = 20) ## Compare output of p3 with theoretical distribution hist(MCMC2, which = "p3", breaks = 20) lines(seq(1, 3, 0.1), dnorm(seq(1, 3, 0.1), mean = 2.5, sd = 0.5)/pnorm(3, 2.5, 0.5)) summary(MCMC2) # functions from package coda... cumuplot(as.mcmc(MCMC2$pars)) summary(as.mcmc(MCMC2$pars)) raftery.diag(MCMC2$pars) ## ======================================================================= ## test 3: sampling a log-normal distribution, log mean=1:4, log sd = 1 ## ======================================================================= NL <- function(p) { mu <- 1:4 -2*sum(log(dlnorm(p, mean = mu, sd = 1))) #-2*log(probability) } MCMCl <- modMCMC(f = NL, p = log(1:4), niter = 3000, outputlength = 1000, jump = 5) plot(MCMCl) # bad convergence cumuplot(as.mcmc(MCMCl$pars)) MCMCl <- modMCMC (f = NL, p = log(1:4), niter = 3000, outputlength = 1000, jump = 2^(2:5)) plot(MCMCl) # better convergence but CHECK it! pairs(MCMCl) colMeans(log(MCMCl$pars)) apply(log(MCMCl$pars), 2, sd) MCMCl <- modMCMC (f = NL, p = rep(1, 4), niter = 3000, outputlength = 1000, jump = 5, updatecov = 100) plot(MCMCl) colMeans(log(MCMCl$pars)) apply(log(MCMCl$pars), 2, sd) ## ======================================================================= ## Fitting a Monod (Michaelis-Menten) function to data ## ======================================================================= # the observations #--------------------- Obs <- data.frame(x=c( 28, 55, 83, 110, 138, 225, 375), # mg COD/l y=c(0.053,0.06,0.112,0.105,0.099,0.122,0.125)) # 1/hour plot(Obs, pch = 16, cex = 2, xlim = c(0, 400), ylim = c(0, 0.15), xlab = "mg COD/l", ylab = "1/hr", main = "Monod") # the Monod model #--------------------- Model <- function(p, x) data.frame(x = x, N = p[1]*x/(x+p[2])) # Fitting the model to the data #--------------------- # define the residual function Residuals <- function(p) (Obs$y - Model(p, Obs$x)$N) # use modFit to find parameters P <- modFit(f = Residuals, p = c(0.1, 1)) # plot best-fit model x <-0:375 lines(Model(P$par, x)) # summary of fit sP <- summary(P) sP[] print(sP) # Running an MCMC #--------------------- # estimate parameter covariances # (to efficiently generate new parameter values) Covar <- sP$cov.scaled * 2.4^2/2 # the model variance s2prior <- sP$modVariance # set nprior = 0 to avoid updating model variance MCMC <- modMCMC(f = Residuals, p = P$par,jump = Covar, niter = 1000, var0 = s2prior, wvar0 = NULL, updatecov = 100) plot(MCMC, Full = TRUE) pairs(MCMC) # function from the coda package. raftery.diag(as.mcmc(MCMC$pars)) cor(MCMC$pars) cov(MCMC$pars) # covariances by MCMC sP$cov.scaled # covariances by Hessian of fit x <- 1:400 SR <- summary(sensRange(parInput = MCMC$pars, func = Model, x = x)) plot(SR, xlab="mg COD/l", ylab = "1/hr", main = "Monod") points(Obs, pch = 16, cex = 1.5)
## ======================================================================= ## Sampling a 3-dimensional normal distribution, ## ======================================================================= # mean = 1:3, sd = 0.1 # f returns -2*log(probability) of the parameter values NN <- function(p) { mu <- c(1,2,3) -2*sum(log(dnorm(p, mean = mu, sd = 0.1))) #-2*log(probability) } # simple Metropolis-Hastings MCMC <- modMCMC(f = NN, p = 0:2, niter = 5000, outputlength = 1000, jump = 0.5) # More accepted values by updating the jump covariance matrix... MCMC <- modMCMC(f = NN, p = 0:2, niter = 5000, updatecov = 100, outputlength = 1000, jump = 0.5) summary(MCMC) plot(MCMC) # check convergence pairs(MCMC) ## ======================================================================= ## test 2: sampling a 3-D normal distribution, larger standard deviation... ## noninformative priors, lower and upper bounds imposed on parameters ## ======================================================================= NN <- function(p) { mu <- c(1,2,2.5) -2*sum(log(dnorm(p, mean = mu, sd = 0.5))) #-2*log(probability) } MCMC2 <- modMCMC(f = NN, p = 1:3, niter = 2000, burninlength = 500, updatecov = 10, jump = 0.5, lower = c(0, 2, 1), upper = c(1, 3, 3)) plot(MCMC2) hist(MCMC2, breaks = 20) ## Compare output of p3 with theoretical distribution hist(MCMC2, which = "p3", breaks = 20) lines(seq(1, 3, 0.1), dnorm(seq(1, 3, 0.1), mean = 2.5, sd = 0.5)/pnorm(3, 2.5, 0.5)) summary(MCMC2) # functions from package coda... cumuplot(as.mcmc(MCMC2$pars)) summary(as.mcmc(MCMC2$pars)) raftery.diag(MCMC2$pars) ## ======================================================================= ## test 3: sampling a log-normal distribution, log mean=1:4, log sd = 1 ## ======================================================================= NL <- function(p) { mu <- 1:4 -2*sum(log(dlnorm(p, mean = mu, sd = 1))) #-2*log(probability) } MCMCl <- modMCMC(f = NL, p = log(1:4), niter = 3000, outputlength = 1000, jump = 5) plot(MCMCl) # bad convergence cumuplot(as.mcmc(MCMCl$pars)) MCMCl <- modMCMC (f = NL, p = log(1:4), niter = 3000, outputlength = 1000, jump = 2^(2:5)) plot(MCMCl) # better convergence but CHECK it! pairs(MCMCl) colMeans(log(MCMCl$pars)) apply(log(MCMCl$pars), 2, sd) MCMCl <- modMCMC (f = NL, p = rep(1, 4), niter = 3000, outputlength = 1000, jump = 5, updatecov = 100) plot(MCMCl) colMeans(log(MCMCl$pars)) apply(log(MCMCl$pars), 2, sd) ## ======================================================================= ## Fitting a Monod (Michaelis-Menten) function to data ## ======================================================================= # the observations #--------------------- Obs <- data.frame(x=c( 28, 55, 83, 110, 138, 225, 375), # mg COD/l y=c(0.053,0.06,0.112,0.105,0.099,0.122,0.125)) # 1/hour plot(Obs, pch = 16, cex = 2, xlim = c(0, 400), ylim = c(0, 0.15), xlab = "mg COD/l", ylab = "1/hr", main = "Monod") # the Monod model #--------------------- Model <- function(p, x) data.frame(x = x, N = p[1]*x/(x+p[2])) # Fitting the model to the data #--------------------- # define the residual function Residuals <- function(p) (Obs$y - Model(p, Obs$x)$N) # use modFit to find parameters P <- modFit(f = Residuals, p = c(0.1, 1)) # plot best-fit model x <-0:375 lines(Model(P$par, x)) # summary of fit sP <- summary(P) sP[] print(sP) # Running an MCMC #--------------------- # estimate parameter covariances # (to efficiently generate new parameter values) Covar <- sP$cov.scaled * 2.4^2/2 # the model variance s2prior <- sP$modVariance # set nprior = 0 to avoid updating model variance MCMC <- modMCMC(f = Residuals, p = P$par,jump = Covar, niter = 1000, var0 = s2prior, wvar0 = NULL, updatecov = 100) plot(MCMC, Full = TRUE) pairs(MCMC) # function from the coda package. raftery.diag(as.mcmc(MCMC$pars)) cor(MCMC$pars) cov(MCMC$pars) # covariances by MCMC sP$cov.scaled # covariances by Hessian of fit x <- 1:400 SR <- summary(sensRange(parInput = MCMC$pars, func = Model, x = x)) plot(SR, xlab="mg COD/l", ylab = "1/hr", main = "Monod") points(Obs, pch = 16, cex = 1.5)
Generates random parameter sets that are (multi)normally distributed.
Norm(parMean, parCovar, parRange = NULL, num)
Norm(parMean, parCovar, parRange = NULL, num)
parMean |
a vector, with the mean value of each parameter. |
parCovar |
the parameter variance-covariance matrix. |
parRange |
the range (min, max) of the parameters, a matrix or a data.frame with one row for each parameter, and two columns with the minimum (1st) and maximum (2nd) column. |
num |
the number of random parameter sets to generate. |
Function Norm
, draws parameter sets from a multivariate normal
distribution, as specified through the mean value and the
variance-covariance matrix of the parameters. In addition, it is
possible to impose a minimum and maximum of each parameter, via
parRange
. This will generate a truncated distribution. Use
this for instance if certain parameters cannot become negative.
a matrix with one row for each generated parameter set, and one column per parameter.
For function Norm
to work, parCovar
must be a valid
variance-covariance matrix. (i.e. positive definite). If this is not the
case, then the function will fail.
Karline Soetaert <[email protected]>
Unif
for uniformly distributed random parameter sets.
Latinhyper
to generates parameter sets using
latin hypercube sampling.
Grid
to generate random parameter sets arranged on a regular
grid
rnorm
the R-default for generating normally distributed random
numbers.
## multinormal parameters: variance-covariance matrix and parameter mean parCovar <- matrix(data = c(0.5, -0.2, 0.3, 0.4, -0.2, 1.0, 0.1, 0.3, 0.3, 0.1, 1.5, -0.7, 1.0, 0.3, -0.7, 4.5), nrow = 4) parCovar parMean <- 4:1 ## Generated sample Ndist <- Norm(parCovar = parCovar, parMean = parMean, num = 500) cov(Ndist) # check pairs(Ndist, main = "normal") ## truncated multinormal Ranges <- data.frame(min = rep(0, 4), max = rep(Inf, 4)) pairs(Norm(parCovar = parCovar, parMean = parMean, parRange = Ranges, num = 500), main = "truncated normal")
## multinormal parameters: variance-covariance matrix and parameter mean parCovar <- matrix(data = c(0.5, -0.2, 0.3, 0.4, -0.2, 1.0, 0.1, 0.3, 0.3, 0.1, 1.5, -0.7, 1.0, 0.3, -0.7, 4.5), nrow = 4) parCovar parMean <- 4:1 ## Generated sample Ndist <- Norm(parCovar = parCovar, parMean = parMean, num = 500) cov(Ndist) # check pairs(Ndist, main = "normal") ## truncated multinormal Ranges <- data.frame(min = rep(0, 4), max = rep(Inf, 4)) pairs(Norm(parCovar = parCovar, parMean = parMean, parRange = Ranges, num = 500), main = "truncated normal")
Plot all observed variables in matrix formalt
obsplot(x, ..., which = NULL, xyswap = FALSE, ask = NULL)
obsplot(x, ..., which = NULL, xyswap = FALSE, ask = NULL)
x |
a Several other matrices or data.frames can be passed in the If the first column of |
which |
the name(s) or the index to the variables that should be plotted. Default = all variables, except the first column. |
ask |
logical; if |
xyswap |
if |
... |
additional arguments. The graphical arguments are passed to The dots may contain other matrices and data.frames with observed data
to be plotted on the same graphs as The arguments after ... must be matched exactly. |
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.
## 'observed' data AIRquality <- cbind(DAY = 1:153, airquality[, 1:4]) head(AIRquality) obsplot(AIRquality, type="l", xlab="Day since May") ## second set of observed data AIR2 <- cbind( 1:100, Solar.R = 250 * runif(100), Temp = 90-30*cos(2*pi*1:100/365) ) obsplot(AIRquality, AIR2, type = "l", xlab = "Day since May" , lwd = 1:2) obsplot(AIRquality, AIR2, type = "l", xlab = "Day since May" , lwd = 1 : 2, which =c("Solar.R", "Temp"), xlim = list(c(0, 150), c(0, 100))) obsplot(AIRquality, AIR2, type = "l", xlab = "Day since May" , lwd = 1 : 2, which =c("Solar.R", "Temp"), log = c("y", "")) obsplot(AIRquality, AIR2, which = 1:3, xyswap = c(TRUE,FALSE,TRUE)) ## ' a data.frame, with 'treatments', presented in long database format Data <- ToothGrowth[,c(2,3,1)] head (Data) obsplot(Data, ylab = "len", xlab = "dose") # same, plotted as two observed data sets obsplot(subset(ToothGrowth, supp == "VC", select = c(dose, len)), subset(ToothGrowth, supp == "OJ", select = c(dose, len)))
## 'observed' data AIRquality <- cbind(DAY = 1:153, airquality[, 1:4]) head(AIRquality) obsplot(AIRquality, type="l", xlab="Day since May") ## second set of observed data AIR2 <- cbind( 1:100, Solar.R = 250 * runif(100), Temp = 90-30*cos(2*pi*1:100/365) ) obsplot(AIRquality, AIR2, type = "l", xlab = "Day since May" , lwd = 1:2) obsplot(AIRquality, AIR2, type = "l", xlab = "Day since May" , lwd = 1 : 2, which =c("Solar.R", "Temp"), xlim = list(c(0, 150), c(0, 100))) obsplot(AIRquality, AIR2, type = "l", xlab = "Day since May" , lwd = 1 : 2, which =c("Solar.R", "Temp"), log = c("y", "")) obsplot(AIRquality, AIR2, which = 1:3, xyswap = c(TRUE,FALSE,TRUE)) ## ' a data.frame, with 'treatments', presented in long database format Data <- ToothGrowth[,c(2,3,1)] head (Data) obsplot(Data, ylab = "len", xlab = "dose") # same, plotted as two observed data sets obsplot(subset(ToothGrowth, supp == "VC", select = c(dose, len)), subset(ToothGrowth, supp == "OJ", select = c(dose, len)))
Fits a model to data, using the pseudo-random search algorithm of Price (1977), a random-based fitting technique.
pseudoOptim(f, p,..., lower, upper, control = list())
pseudoOptim(f, p,..., lower, upper, control = list())
f |
function to be minimised, its first argument should be the vector of parameters over which minimization is to take place. It should return a scalar result, the model cost, e.g the sum of squared residuals. |
p |
initial values of the parameters to be optimised. |
... |
arguments passed to function |
lower |
minimal values of the parameters to be optimised; these must be specified; they cannot be -Inf. |
upper |
maximal values of the parameters to be optimised; these must be specified; they cannot be +Inf. |
control |
a list of control parameters - see details. |
The control
argument is a list that can supply any of the
following components:
npop, number of elements in the population. Defaults to max(5*length(p),50).
numiter, maximal number of iterations to be performed.
Defaults to 10000. The algorithm either stops when numiter
iterations has been performed or when the remaining variation is less
than varleft
.
centroid, number of elements from which to estimate a new parameter vector, defaults to 3.
varleft, relative variation remaining; if below this value the algorithm stops; defaults to 1e-8.
verbose, if TRUE, more verbose output will contain the parameters
in the final population, their respective population costs and the cost
at each succesful interation. Defaults to FALSE
.
see the book of Soetaert and Herman (2009) for a description of the algorithm AND for a line to line explanation of the function code.
a list containing:
par |
the optimised parameter values. |
cost |
the model cost, or function evaluation associated to the optimised parameter values, i.e. the minimal cost. |
iterations |
the number of iterations performed. |
and if control$verbose
is TRUE:
poppar |
all parameter vectors remaining in the population, matrix of dimension (npop,length(par)). |
popcost |
model costs associated with all population parameter vectors, vector of length npop. |
rsstrace |
a 2-columned matrix with the iteration number and the model cost at each succesful iteration. |
Karline Soetaert <[email protected]>
Soetaert, K. and Herman, P. M. J., 2009. A Practical Guide to Ecological Modelling. Using R as a Simulation Platform. Springer, 372 pp.
Price, W.L., 1977. A Controlled Random Search Procedure for Global Optimisation. The Computer Journal, 20: 367-370.
amp <- 6 period <- 5 phase <- 0.5 x <- runif(20)*13 y <- amp*sin(2*pi*x/period+phase) + rnorm(20, mean = 0, sd = 0.05) plot(x, y, pch = 16) cost <- function(par) sum((par[1] * sin(2*pi*x/par[2]+par[3])-y)^2) p1 <- optim(par = c(amplitude = 1, phase = 1, period = 1), fn = cost) p2 <- optim(par = c(amplitude = 1, phase = 1, period = 1), fn = cost, method = "SANN") p3 <- pseudoOptim(p = c(amplitude = 1, phase = 1, period = 1), lower = c(0, 1e-8, 0), upper = c(100, 2*pi, 100), f = cost, control = c(numiter = 3000, verbose = TRUE)) curve(p1$par[1]*sin(2*pi*x/p1$par[2]+p1$par[3]), lty = 2, add = TRUE) curve(p2$par[1]*sin(2*pi*x/p2$par[2]+p2$par[3]), lty = 3, add = TRUE) curve(p3$par[1]*sin(2*pi*x/p3$par[2]+p3$par[3]), lty = 1, add = TRUE) legend ("bottomright", lty = c(1, 2, 3), c("Price", "Mathematical", "Simulated annealing"))
amp <- 6 period <- 5 phase <- 0.5 x <- runif(20)*13 y <- amp*sin(2*pi*x/period+phase) + rnorm(20, mean = 0, sd = 0.05) plot(x, y, pch = 16) cost <- function(par) sum((par[1] * sin(2*pi*x/par[2]+par[3])-y)^2) p1 <- optim(par = c(amplitude = 1, phase = 1, period = 1), fn = cost) p2 <- optim(par = c(amplitude = 1, phase = 1, period = 1), fn = cost, method = "SANN") p3 <- pseudoOptim(p = c(amplitude = 1, phase = 1, period = 1), lower = c(0, 1e-8, 0), upper = c(100, 2*pi, 100), f = cost, control = c(numiter = 3000, verbose = TRUE)) curve(p1$par[1]*sin(2*pi*x/p1$par[2]+p1$par[3]), lty = 2, add = TRUE) curve(p2$par[1]*sin(2*pi*x/p2$par[2]+p2$par[3]), lty = 3, add = TRUE) curve(p3$par[1]*sin(2*pi*x/p3$par[2]+p3$par[3]), lty = 1, add = TRUE) legend ("bottomright", lty = c(1, 2, 3), c("Price", "Mathematical", "Simulated annealing"))
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
and where is an output variable (at a certain time instance),
is a parameter, and
is the scaling of variable
,
is the scaling of parameter
.
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), ...)
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), ...)
func |
an R-function 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 |
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 is perturbed as
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.
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 ,
L2: the L2-norm ,
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
.
Karline Soetaert <[email protected]>
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. 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. doi:10.18637/jss.v033.i03
## ======================================================================= ## 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)])
## ======================================================================= ## 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)])
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.
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, ...)
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, ...)
func |
an R-function that has as first argument |
parms |
parameters passed to |
sensvar |
the output variables for which the sensitivity needs
to be estimated. Either |
dist |
the distribution according to which the parameters should be
generated, one of |
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 |
parInput |
a matrix with dimension (*, npar) with the values of the sensitivity parameters. |
parMean |
only when |
parCovar |
only when |
num |
the number of times the model has to be run. Set large
enough. If |
map |
the column number with the (independent) mapping variable
in the output matrix returned by |
object |
an object of class |
x |
an object of class |
legpos |
position of the legend; set to |
xyswap |
if |
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 |
ask |
logical; if |
obs |
a |
obspar |
additional graphics arguments passed to |
... |
additional arguments passed to |
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
.
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.
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".
Karline Soetaert <[email protected]>
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. doi:10.18637/jss.v033.i03
## ======================================================================= ## 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))
## ======================================================================= ## 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))
Generates uniformly distributed random parameter sets.
Unif(parRange, num)
Unif(parRange, num)
parRange |
the range (min, max) of the parameters, a matrix or a data.frame with one row for each parameter, and two columns with the minimum (1st) and maximum (2nd) value. |
num |
the number of random parameter sets to generate. |
In the uniform sampling, each parameter is uniformly random distributed over its range.
a matrix with one row for each generated parameter set, and one column per parameter.
For small sample sizes, the latin hypercube distributed parameter sets
(Latinhyper
) may give better coverage in parameter space
than the uniform random design.
Karline Soetaert <[email protected]>
Norm
for (multi)normally distributed random parameter sets.
Latinhyper
to generates parameter sets using
latin hypercube sampling.
Grid
to generate random parameter sets arranged on a
regular grid
runif
the R-default for generating uniformally distributed
random numbers.
## 4 parameters parRange <- data.frame(min = c(0, 1, 2, 3), max = c(10, 9, 8, 7)) rownames(parRange) <- c("par1", "par2", "par3", "par4") ## uniform pairs(Unif(parRange, 100), main = "Uniformly random")
## 4 parameters parRange <- data.frame(min = c(0, 1, 2, 3), max = c(10, 9, 8, 7)) rownames(parRange) <- c("par1", "par2", "par3", "par4") ## uniform pairs(Unif(parRange, 100), main = "Uniformly random")