Dear All,

please help with the following:

I have model:

require(deSolve)
require(FME)
pars <- list(cl = 0.348,v=14.3,EC50=4.1)

DOSE <-5
eventdat <- data.frame(var = rep("a", 5),time = seq(0,96,by=24) ,value = rep(DOSE,5),method = rep("add",5))
INRB <-0.89

derivs <- function(t, state, pars) {
     dy1 <- - 2 * state[1]
     dy2 <- 2 * state[1] - (pars$cl/pars$v) * state[2] 
     EDK50 <-pars$EC50*pars$cl
     DR <- (pars$cl/pars$v) * state[2]
     EFF <-(1*(DR^1.15))/((EDK50^1.15)+(DR^1.15))
     dy3 <-(1-EFF)*(3/28.6)-state[3]*(3/28.6)
     dy4 <-state[3]*(3/28.6)-state[4]*(3/28.6)
     dy5 <-state[4]*(3/28.6)-state[5]*(3/28.6)
     dy6 <-(1-EFF)*(3/118.3)-state[6]*(3/118.3)
     dy7 <-state[6]*(3/118.3)-state[7]*(3/118.3)
     dy8 <-state[7]*(3/118.3)-state[8]*(3/118.3)
     
     dc <-INRB + 20*(1-((state[5]+state[8])/2))   
     
     return(list(c(dy1,dy2,dy3,dy4,dy5,dy6,dy7,dy8),z=dc))

}

model <- function(pars, times=seq(0, 120, by = 0.01)) {
  state <- c(a = 0,b=0,c=1,d=1,e=1,f=1,g=1,h=1)
  return(ode(y = state, times = times, func = derivs, parms = pars,events = list(data = eventdat),method="lsoda",rtol = 1e-8,atol =1e-8,hmax=0.01))
}


plot(model(pars))
observed <- matrix (nc=2,byrow=2,data=c(0,1,30,1.1,55,1.2))
colnames(observed) <- c("time", "z.e")

tout <- unique(sort(c(eventdat$time,observed[,"time"])))

obj <- function(x, parset = names(x)) {
 pars[parset] <- x
 out <- model(pars, tout)
 return(modCost(obs = observed, model = out))
}

prior <- function(p)
  return( sum(((p-c(0.348,14.3,4.1))/c(0.348*0.298,14.3*0.232,4.1*0.332))^2 ))

final <- modMCMC(p = c(cl=pars$cl,v=pars$v,EC50=pars$EC50), f = obj, prior=prior,lower=c(0.000001,0.00000001,0.00000001),jump = NULL, niter = 200,updatecov=10, wvar0 = 1,var0=NULL,burninlength = 10)


summary(final)
plot(final)
pars2 <- as.list(final$bestpar)

out1 <- model(pars)
out2 <- model(pars2)

plot(out1, out2, obs=observed)
legend("topright", legend=c("start parameters", "fitted"), col=1:2, lty=1:2)

I have 2 questions:

1. for unknown reason to me the name of my global value in the output is changed to "z.e", instead of getting it back as z. For some reason the name of the 5th state variable (here it is "e") is automatically added to the name of my global value

2. I can not figure out why the modMCMC function is not running properly for me (ie: I have an error somewhere I think). it will initiate and run the chains but my bestpar at the end is the same as the starting point... 

All of your input is greatly appreciated,

thanks,

Andras

	[[alternative HTML version deleted]]

