[R-sig-dyn-mod] ] Bayesian model fitting and deSolve correction

[Andras Farkas]

Correction on the code below, all the rlnorm should be rnorm. Sorry about that,

Andras



------------------------------
On Wed, Apr 24, 2013 3:54 PM EDT Andras Farkas wrote:

>Dear List,
> 
>If you can, please help with writing the code for the following problem:
>Currently I have (and thanks to Dr Petzoldt):
> 
>require(deSolve)
>pkmod <- function(t, y, p) {
>  inp <- forc(t)
>  dy1 <- - (p["k12"]+p["k"])* y[1] + p["k21"]*y[2]   + inp
>  dy2 <- p["k12"]* y[1] - p["k21"] *y[2]
>  list(c(dy1, dy2), inp=inp)
>}
>intimes <- c(0,0.5, 4, 4.5, 10,10.5, 24)
>input   <- c(800,0, 200,0, 100, 0,0)
>forc <- approxfun(intimes, input, method="constant")
>times <- seq(0, 24, 0.1)
>yini <- c(Central = 0, Peripheral = 0)
>p    <- c(k=0.05, k12=0.197118, k21=0.022665)
>result <- lsoda(yini, times, pkmod, p, rtol = 1e-08, atol = 1e-08)
>result <- data.frame(result)
>names(result) <- c("time","Amount1", "Amount2")
>plot(result$Amount1 ~ result$time, type="b", main="Central compartment", xlab="time", ylab="Amount")
> 
>this will produce the time vs amount 'Central compartment" graph...What I am trying to figure out is to fit the "pkmod" model to a few values extracted from time vs amount in the 'Central compartment" plot to find parameters p, basically doing the whole thing "backwards" with a bayesian model fitting approach to find the model paramaters "p" as best as I can. I am rather new to deSolve, and I used OpenBugs via R2openBugs before, but would like to convert to R based computations fully. So a few data points from the plot "Central compartment" extracted are:
> 
>time <-c(1,2,3,4,6,10,12)
>amount <-c(332,261.7,207.1,165.2,173.7,82.05,94.6)
> 
>here we also know that all model parameters follow normal distribution and can be simulated with the following relationship (this is the way I do it in OpenBugs to avoid from negative values being generated)
> 
>kpre=0.05
>kvar <-rlnorm(5000,0,0.3)
>k <-kpre*exp(kvar) # parameter 1
> 
>k12pre=0.197118
>k12var <-rlnorm(5000,0,0.3)
>k12 <-k12pre*exp(k12var)  # parameter 2
> 
>k21pre=0.022665
>k21var <-rlnorm(5000,0,0.3)
>k21 <-k21pre*exp(k21var) # parameter 3
> 
>Well, for OpenBugs this information would be enough to run, and I would just run it through my code and get the stats for the fitted parameters. Any input on how to approach the problem would be appreciated. I looked at laplacesDemon, which I found to be a bit complicated for my level of understanding, looked also at FME, but not sure if that does bayesian model fitting to find model parameters, it seems like the fitting function are likely different. One another point I would like to mention is that at times I have only 1 time vs amount data point to fit, which often has "OpenBugs hung up". The same thing does not happen when I use analytical solutions, so perhaps it is related to the diff equations?
> 
>Appreciate all the help I can get,
> 
>Sincerely,
> 
>Andras
>	[[alternative HTML version deleted]]
>