fitOdeModel {simecol} | R Documentation

## 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"),
      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 and `optim` for all other methods,
`lower, upper` | bounds of the parameters for method L-BFGS-B, see `optim` and for method PORT see `nlminb`. 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.

### Note

The interface to the PORT algorithm will change for future R versions.

### References

Gay, D. M. (1990) Usage Summary for Selected Optimization Routines. Computing Science Technical Report No. 153. AT&T Bell Laboratories, Murray Hill, NJ. http://netlib.bell-labs.com/cm/cs/cstr/153.pdf

### See Also

`ssqOdeModel`, `optim`

### Examples

    library(simecol)

    ## ======== 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))

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

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

    res$par

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

    res$par

    ## alternative 3a, use PORT algorithm
    ## Note: lower and upper are *named* vectors
    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)

    res$par

    ## alternative 3b, PORT algorithm with manual parameter scaling
    ## Note: lower and upper are *named* vectors

    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)

    res$par

    ## set model parameters to  fitted values and simulate again
    parms(cs2)[whichpar] <- res$par
    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")

[Package simecol version 0.6-9 Index]