[R-sig-dyn-mod] Turbidostat

Uta,

thanks for the nice example. I would do it like follows: declare the 
correction factors (i.e. time dependent parameters) as parameters for 
the modCost function, so that they can be fitted like ordinary model 
parameters. Note also that in the example below we now pass the "signal" 
as an optional parameter explicitly to the turbidostat function.

A few additional notes:

1) The parameters may be not simultaneously identifiable because of 
strong correlation. Ypou may consider to use other start values, e.g.

parms <- c(r=1, korr_s=1)

and/or to adapt the interpolation method (e.g. "constant").

2) Why does the simulation start at time step 1 and not at zero like the 
data?

3) Is there any evidence why this inhibition occurs at the specified 
times only? Is there any chance to describe this a little bit more 
process oriented?

Thomas P.



###########################################################
library(deSolve)
library(FME)

maxTime <- 15
korr_s <- 0.3  #not fitted but badly assumed - just as an example
signal <- approxfun(1:maxTime, c(rep(korr_s,3),rep(1,10),rep(korr_s,2)))

turbidostat <- function(time, E0, parms, signal){
     korr <- signal(time)
     with(as.list(c(E0, parms)), {
       dE <- r*korr*E
     list(c(dE))
   })
}

E0 <- c(E=0.1)
parms <- c(r=0.5)
times <- seq(1,15,0.1)

eventdat <- data.frame(var = c("E","E","E","E","E"), time = c(7,9,11,13,15),
   value=c(rep(0.5,5)),method=c("rep","rep","rep","rep","rep"))


#####observed data
Obs <- (data.frame(time = c(0,1,2,3,7,9,11,13,15),
                    E = c(0.1,0.11,0.2,0.3,1.8,1.4,1.45,1.46,1.0)))

out <- ode(E0,times,turbidostat,parms, events=list(data=eventdat), 
signal=signal)

plot(Obs$time,Obs$E,pch=16,cex=2.0)
lines(out, type = "l", lwd = 2)

Cost <- function(parms){
   pks <- parms["korr_s"]
   signal <- approxfun(1:maxTime, c(rep(pks, 3),rep(1, 10),rep(pks, 2)))
   model <- ode(E0, times, turbidostat, parms,
     events=list(data=eventdat), signal=signal)
   modCost(model = model, obs = Obs, x = "time")
}

parms <- c(r=0.5, korr_s=0.3)

Cost(parms)

Fit <- modFit(f=Cost,p=parms)
deviance(Fit)
Fit$par

out <- ode(E0,times,turbidostat,Fit$par, events=list(data=eventdat), 
signal=signal)

plot(Obs$time,Obs$E,pch=16,cex=2.0)
lines(out, type = "l", lwd = 2)

## But note that there is *very strong* correlation
## between these two parameters:
summary(Fit)