[R-sig-dyn-mod] how to fit parameters of simecol models

[Thomas Petzoldt]

Hello,

I do not feel very comfortable with the posted example. Nevertheless, 
parameter fitting for simecol objects seems of general interest, so here 
are a few notes and a pharmacokinetic model.

Package simecol already had a preliminary model fitting function 
(fitOdeModel) but the recent package FME (Flexible Modeling Environment) 
has much more flexible fitting procedures. So it is possible to define 
an own cost function, that may use either the full or only parts of the 
data set. In addition, one gets better information about the model fit 
(using summary).

Details are described in our JSS article:
http://cran.r-project.org/web/packages/FME/

and in the package vignettes:
http://cran.r-project.org/web/packages/FME/


An example for how to use FME with simecol is attached below.

Thomas Petzoldt


#####################################################

library(simecol)

## A two compartment pharmacokinetic model
twocomp <- new("odeModel",
   main = function (time, init, parms, ...) {
     with(as.list(c(parms, init)), {
       dCl <- - (Ke + Klf) * Cl + Kfl * Cf  #Liver
       dCf <-    Klf * Cl  - Kfl * Cf       #Fat
       list(c(dCl, dCf))
     })
   },
   parms = c(Ke = 0.1, Kfl = 0.2, Klf = 0.02),
   times = 0:40,
   init =  c(Cl = 1, Cf = 0),
   solver = "lsoda"
)

## simulate and plot the model
twocomp <- sim(twocomp)
plot(twocomp)

## use data from simulation and add some noise
dat <- out(twocomp)[seq(1, 30, 4), ]
N <- nrow(dat)
dat$Cl <- dat$Cl + rnorm(N, sd=0.1)
dat$Cf <- dat$Cf + rnorm(N, sd=0.002)

## === Now let's fit the model =========================
library(FME) ## see http://www.jstatsoft.org/v33/i03

## define a cost function
ModelCost <- function(p) {
  parms(twocomp) <- p
  out   <- out(sim(twocomp))
  cost <- modCost(out, dat, weight="std")

  ## OR: stepwise build up of cost function
  #cost1 <- modCost(out, dat[c("time", "Cl")], weight="std")
  #cost <- modCost(out, dat[c("time", "Cf")], weight="std", cost=cost1)
  cost
}


## set times to the time steps of data
times(twocomp) <- dat$time

## disturb parameters, so that we see something
parms(twocomp) <- c(Ke = 0.2,	Kfl = 0.1,	Klf = 0.05)

## Fit the model
Fit <- modFit(f = ModelCost, p = parms(twocomp))
summary(Fit)

## Now simulate and plot model and data
times(twocomp) <- 0:40    # restore original time steps

twocomp.fitted <- twocomp # make clone of model
parms(twocomp.fitted) <- coef(Fit) # assign fitted parms

## simulate original and fitted model
twocomp <- sim(twocomp)
twocomp.fitted <- sim(twocomp.fitted)

oo <- out(twocomp)
of <- out(twocomp.fitted)

## plot results
par(mfrow = c(2,1))
plot(oo$time, oo$Cl, type="l")
lines(of$time, of$Cl, col="blue")
points(dat$time, dat$Cl, col="red", pch="+")

plot(oo$time, oo$Cf, type="l")
lines(of$time, of$Cf, col="blue")
points(dat$time, dat$Cf, col="red", pch="+")