Message-ID: <4CFFAE1B.6090402@tu-dresden.de>
Date: 2010-12-08T16:11:07Z
From: Thomas Petzoldt
Subject: [R-sig-dyn-mod] Turbidostat
In-Reply-To: <545390A0-D865-4FF0-9A26-107366FE61BB@tu-dresden.de>
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)