Dear All,
 
Currently I am running the following code:
 
library(stats4)
library(odesolve)
library(rgenoud)
Input<-data.frame(SUB=c(1),time=c(0.5,3,10,15),lev=c(2.05,12.08,9.02,8))
XD<-500
IT<-3
diffeqfun<-function(time, y, parms) { 
  if(time<=IT)  
    dCpdt <- (XD/IT)/parms["Vol"] - (parms["Vm"]/parms["Vol"])*y[1]/(parms["Km"]/parms["Vol"]+y[1]) 
  else
    dCpdt <- -(parms["Vm"]/parms["Vol"])*y[1]/(parms["Km"]/parms["Vol"]+y[1])
  list(dCpdt)
}
modelfunc<-function(time,Vm,Km,Vol) { 
  out <- lsoda(0,c(0,time),diffeqfun,parms=c(Vm=Vm,Km=Km,Vol=Vol),rtol=1e-5,atol=1e-5)
  out[-1,2] 
} 
objectfunc <- function(par) {
  out <- modelfunc(Input$time, par[1], par[2], par[3])
  gift <- which( Input$lev != 0 )
  sum((Input$lev[gift]-out[gift])^2)
}        
gen <- genoud(objectfunc,nvars=3,max=FALSE,pop.size=10,max.generations=100,wait.generations=100, starting.value=c(40,8,12),BFGS=FALSE, print.level=0,boundary.enforcement=0,Domains=matrix(c(0,0,0,100,100,100),3,2),MemoryMatrix=TRUE)
outgen<-c(gen$par[1],gen$par[2],gen$par[3])
opt<-optim(c(gen$par[1],gen$par[2],gen$par[3]),objectfunc,method="Nelder-Mead")  
outopt<-c(opt$par[1],opt$par[2],opt$par[3])     
fm<-nls(lev~modelfunc(time,Vm,Km,Vol),data=Input,start=list(Vm=opt$par[1],Km=opt$par[2],Vol=opt$par[3]),trace=TRUE,nls.control(tol=1))
x<-Input$time
y<-Input$lev
cal<-predict(fm,list(time=x))
plot(x,y,type="l",col="green")
lines(x,cal,col="red")
 
as you can see, the parameter optimization using this code is done via a non-linear approach using the Nelder-Mead optimization. I was wondering if you could give me some ideas on how to rewrite this code so that it could be run using OpenBUGS via R2OpenBUGS? The issue I am hung up on is how to write the differential equation into a code that could run through R2OpenBUGS and OpenBUGS. Below is an example for an analytical solution for a different model I currently use that is run through R2OpenBUGS and OpenBUGS that generates excelent results (to be able to run this code OpenBUGS installed, also please note this is only an example, I really need to reproduce the ideas from above using the bayesian approach):
 
library(R2OpenBUGS)
D <-400
tin <-1
VDPRE <-0.5
LWEIGHT <-68
CRCL <-60
tau <-12
T <-2
ts <-c(2,8)
c <-c(15,6)
mfsssc <- function() {
           for (i in 1:1) { 
           for (j in 1:T) {            
           c[j]~dnorm(mu[j],100)      
           mu[j]<-((D/tin)*(1/cl))*(((1-exp(-(cl/v)*tin))*exp(-(cl/v)*(ts[j])))/(1-exp(-(cl/v)*tau)))  
           }    
           cla~dnorm(0, 6.25)
           va~dnorm(0, 6.25)
           theta1<-0.615*(CRCL*0.06)+0.7194      
           theta2<-VDPRE*LWEIGHT 
           cl<-exp(cla)*theta1
           v<-exp(va)*theta2              
       }
}
res <- bugs(list(T = T, c = c, tau = tau, ts = ts, tin = tin, D = D, CRCL = CRCL, VDPRE = VDPRE, LWEIGHT = LWEIGHT), model.file = mfsssc, inits = NULL, parameters.to.save=c("v", "cl"), n.chains=2, n.iter=10000)
finalsssc <- res$summary
 
I hope my question made enough sense, but if not please ask and I will try to further explain my goals here. All your input would be greatly apreciated,
 
Sincerely,
 
Andras
	[[alternative HTML version deleted]]


