[R-sig-dyn-mod] modMCMC question

[McGrath, Justin M]

z, as far as I can tell, is not a global variable. It does not exist in the global environment.
It has "e" appended because it is keeping the name of state[5] in the following line 
dc <-INRB + 20*(1-((state[5]+state[8])/2))
You can stop that from happening by using this return statement
 return(list(c(dy1,dy2,dy3,dy4,dy5,dy6,dy7,dy8),z=as.numeric(dc)))  # as.numeric() returns "dc" without attributes.
or by adding "names(dc) <- NULL" before your current return statement. I don't know which is faster, but I'd guess as.numeric() is.

I don't have any help for modMCM()

Justin
________________________________________
From: r-sig-dynamic-models-bounces at r-project.org [r-sig-dynamic-models-bounces at r-project.org] on behalf of Andras Farkas [motyocska at yahoo.com]
Sent: Saturday, February 22, 2014 8:12 PM
To: r-sig-dynamic-models at r-project.org
Subject: [R-sig-dyn-mod] modMCMC question

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]]