fitOdeModel {simecol}

Parameter Fitting for odeModel Objects

Description

Fit parameters of odeModel objects to measured data.

Usage

fitOdeModel(simObj, whichpar = names(parms(simObj)), obstime, yobs, 
  sd.yobs = as.numeric(lapply(yobs, sd)), initialize = TRUE, 
  weights = NULL, debuglevel = 0, fn = ssqOdeModel, 
  method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "PORT",
   "newuoa", "bobyqa"),
  lower = -Inf, upper = Inf, scale.par = 1,
  control = list(), ...)

Arguments

simObj

a valid object of class odeModel,

whichpar

character vector with names of parameters which are to be optimized (subset of parameter names of the simObj),

obstime

vector with time steps for which observational data are available,

yobs

data frame with observational data for all or a subset of state variables. Their names must correspond exacly with existing names of state variables in the odeModel,

sd.yobs

vector of given standard deviations (or scale) for all observational variables given in yobs. If no standard deviations (resp. scales) are given, these are estimated from yobs,

initialize

optional boolean value whether the simObj should be re-initialized after the assignment of new parameter values. This can be necessary in certain models to assign consistent values to initial state variables if they depend on parameters.

weights

optional weights to be used in the fitting process. See cost function (currently only ssqOdeModel) for details.

debuglevel

a positive number that specifies the amount of debugging information printed,

fn

objective function, i.e. function that returns the quality criterium that is minimized, defaults to ssqOdeModel,

method

optimization method, see nlminb for the PORT algorithm, newuoa resp. bobyqa for the newuoa and bobyqa algorithms, and optim for all other methods,

lower, upper

bounds of the parameters for method L-BFGS-B, see optim, PORT see nlminb and bobyqa bobyqa. The bounds are also respected by other optimizers by means of an internal transformation of the parameter space (see p.constrain). In this case, named vectors are required.

scale.par

scaling of parameters for method PORT see nlminb. In many cases, automatic scaling (scale.par = 1) does well, but sometimes (e.g. if parameter ranges differ several orders of magnitude) manual adjustment is required. Often you get a reasonable choice if you set scale.par = 1/upper. The parameter is ignored by all other methods. For "Nelder-Mead", "BFGS", "CG" and "SANN" parameter scaling occurs as a side effect of parameter transformation with p.constrain.

control

a list of control parameters for optim resp. nlminb,

...

additional parameters passed to the solver method (e.g. to lsoda).

Details

This function works currently only with odeModel objects where parms is a vector, not a list.

Note also that the control parameters of the PORT algorithm are different from the control parameters of the other optimizers.

Value

A list with the optimized parameters and other information, see optim resp. nlminb for details.

References

Gay, D. M. (1990) Usage Summary for Selected Optimization Routines. Computing Science Technical Report No. 153. AT&T Bell Laboratories, Murray Hill, NJ.

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

See Also

ssqOdeModel, optim, nlminb, bobyqa

Note also that package FME function modFit has even more flexible means to fit model parameters.

Examples are given in the package vignettes.

Examples

## ======== load example model =========
data(chemostat)

#source("chemostat.R")

## derive scenarios
cs1 <- cs2 <- chemostat

## generate some noisy data
parms(cs1)[c("vm", "km")] <- c(2, 10)
times(cs1) <- c(from=0, to=20, by=2)
yobs <- out(sim(cs1))
obstime <- yobs$time
yobs$time <- NULL
yobs$S <- yobs$S + rnorm(yobs$S, sd= 0.1 * sd(yobs$S))*2
yobs$X <- yobs$X + rnorm(yobs$X, sd= 0.1 * sd(yobs$X))

## ======== optimize it! =========

## time steps for simulation, either small for rk4 fixed step
# times(cs2)["by"] <- 0.1
# solver(cs2) <- "rk4"

## or, faster: use lsoda and and return only required steps that are in the data
times(cs2) <- obstime
solver(cs2) <- "lsoda"

## Nelder-Mead (default)
whichpar  <- c("vm", "km")

res <- fitOdeModel(cs2, whichpar=whichpar, obstime, yobs,
  debuglevel=0,
  control=list(trace=TRUE))

coef(res)

## assign fitted parameters to the model, i.e. as start values for next step
parms(cs2)[whichpar] <- coef(res)

## alternatively, L-BFGS-B (allows lower and upper bounds for parameters)
res <- fitOdeModel(cs2, whichpar=c("vm", "km"), obstime, yobs,
  debuglevel=0, fn = ssqOdeModel,
  method = "L-BFGS-B", lower = 0,
  control=list(trace=TRUE),
  atol=1e-4, rtol=1e-4)

coef(res)

## alternative 2, transform parameters to constrain unconstrained method
## Note: lower and upper are *named* vectors
res <- fitOdeModel(cs2, whichpar=c("vm", "km"), obstime, yobs,
  debuglevel=0, fn = ssqOdeModel,
  method = "BFGS", lower = c(vm=0, km=0), upper=c(vm=4, km=20),
  control=list(trace=TRUE),
  atol=1e-4, rtol=1e-4)

coef(res)


## alternative 3a, use PORT algorithm
parms(cs2)[whichpar] <- c(vm=1, km=2)

lower <- c(vm=0, km=0)
upper <- c(vm=4, km=20)

res <- fitOdeModel(cs2, whichpar=c("vm", "km"), obstime, yobs,
  debuglevel=0, fn = ssqOdeModel,
  method = "PORT", lower = lower, upper = upper,
  control=list(trace=TRUE),
  atol=1e-4, rtol=1e-4)

coef(res)

## alternative 3b, PORT algorithm with manual parameter scaling
res <- fitOdeModel(cs2, whichpar=c("vm", "km"), obstime, yobs,
  debuglevel=0, fn = ssqOdeModel,
  method = "PORT", lower = lower, upper = upper, scale.par = 1/upper,
  control=list(trace=TRUE),
  atol=1e-4, rtol=1e-4)

coef(res)

## set model parameters to  fitted values and simulate again
parms(cs2)[whichpar] <- coef(res)
times(cs2) <- c(from=0, to=20, by=1)
ysim <- out(sim(cs2))

## plot results
par(mfrow=c(2,1))
plot(obstime, yobs$X, ylim = range(yobs$X, ysim$X))
lines(ysim$time, ysim$X, col="red")
plot(obstime, yobs$S, ylim= range(yobs$S, ysim$S))
lines(ysim$time, ysim$S, col="red")