metrumresearchgroup / pbpk-qsp-mrgsolve

http://metrumrg.com/events/
24 stars 14 forks source link

[lsoda] error #8

Open FanWn opened 2 years ago

FanWn commented 2 years ago

@kylebaron Hi Kyle, I'm trying to optimize two parameters of my PBPK model against some observed data (11 mean concentration data points). The codes are as follow:

data<-read.csv("data.csv")

yobs <- filter(data,EVID==0) %>% dplyr::select(DV) %>% unlist %>% unname
yobs

wss <- function(dv, pred, par=NULL) {
  sum(((dv-pred)/dv)^2)
}

pred <- function(.par, .data, yobs=NULL, pred=FALSE) {

  .par <- .par %>% setNames(names(theta))

  .mod <- param(mod,.par)

  if(pred) {
    out <- mrgsim(.mod,data=.data,carry.out="type", maxsteps = 10000)
    return(as_data_frame(out))
  }

  out <- mrgsim(.mod, data=.data, obsonly=TRUE, Req="Cvenous", maxsteps = 10000)

  return(wss(yobs,out$Cvenous))

  #return(-1*sum(dnorm(log(yobs),log(out$CP),.par$sigma,log=TRUE)))

}

theta<-c(CL_r=9.8, Kp=0.001)
fit1b <- optim(theta, pred,.data=data, yobs=yobs)

Then I got the lsoda error:

[lsoda] warning..internal t = 1.25203 and h = 6.11389e-17 are such that in the machine, t + h = t on the next step solver will continue anyway. [lsoda] warning..internal t = 1.25203 and h = 6.11389e-17 are such that in the machine, t + h = t on the next step solver will continue anyway. [lsoda] above warning has been issued 2 times it will not be issued again for this problem. [lsoda] 10000 steps taken before reaching tout; consider increasing maxsteps.

[mrgsolve] lsoda returned with negative istate: -1 excess work done on this call; check the model or increase maxsteps. current value of maxsteps: 10000

Error: simulation terminated.

I tried to increase the maxsteps value but that didn't work and the same error occurred. Can you tell what is wrong? Is it because there is not enough observed data?

kylebaron commented 2 years ago

Hi @FanWn -

warning..internal t = 1.25203 and h_ = 6.11389e-17

It looks like the ode solver is having an issue at time = 1.25 so I don't think it should take that many steps to get there.

The first thing I'd try doing is log-transforming the parameter list

data<-read.csv("data.csv")

yobs <- filter(data,EVID==0) %>% dplyr::select(DV) %>% unlist %>% unname
yobs

wss <- function(dv, pred, par=NULL) {
  sum(((dv-pred)/dv)^2)
}

pred <- function(.par, .data, yobs=NULL, pred=FALSE) {

  .par <- lapply(.par, exp)

  .par <- .par %>% setNames(names(theta))

  .mod <- param(mod,.par)

  if(pred) {
    out <- mrgsim(.mod,data=.data,carry.out="type", maxsteps = 10000)
    return(as_data_frame(out))
  }

  out <- mrgsim(.mod, data=.data, obsonly=TRUE, Req="Cvenous", maxsteps = 10000)

  return(wss(yobs,out$Cvenous))

  #return(-1*sum(dnorm(log(yobs),log(out$CP),.par$sigma,log=TRUE)))

}

theta<- log(c(CL_r=9.8, Kp=0.001))
fit1b <- optim(theta, pred,.data=data, yobs=yobs)

Can you try that and let's go from there.

Kyle

FanWn commented 2 years ago

Hi Kyle @kylebaron , thanks for troubleshooting! The codes worked with what you suggested (only the function for pred should be pred <- function(par, .data, yobs=NULL, pred=FALSE) {). However, the simulation is still off after I plugged in the optimized parameter values. I tried a couple of different initials but didn't improve.

kylebaron commented 2 years ago

Ok; thanks for reporting back. Next, I'd look to see what might be going on in the data or model; this will need to be shared so we can take a look.