Giter Site home page Giter Site logo

Comments (2)

Vincent-AC avatar Vincent-AC commented on June 3, 2024

It is especially surprising because when you remove the 10^-5 added to every time-point and set ourzero to 0 the OFVs make sense again.

poped.db$design$xt

       obs_1 obs_2 obs_3 obs_4
grp_1      0    10    20  22.5
grp_2      0    10    20  22.5
grp_3      0    10    20  22.5
...

poped.db2$design$xt

       obs_1 obs_2 obs_3
grp_1     10    20  22.5
grp_2     10    20  22.5
grp_3     10    20  22.5
...

evaluate_design(poped.db)

$ofv
[1] 44.11343

$fim
                    KNET       BMAX       INOC ESLOPE_PMB_NR  GAMMA_PMB SIGMA[1,1]
KNET           53703.462  6729.2954  2715.7124    -54735.118  71330.641      0.000
BMAX            6729.295   996.8936   376.4853     -6851.225  10780.754      0.000
INOC            2715.712   376.4853   812.1461     -2788.706   3319.721      0.000
ESLOPE_PMB_NR -54735.118 -6851.2252 -2788.7059     55823.626 -70947.690      0.000
GAMMA_PMB      71330.641 10780.7541  3319.7208    -70947.690 203222.874      0.000
SIGMA[1,1]         0.000     0.0000     0.0000         0.000      0.000   6840.097

$rse
         KNET          BMAX          INOC ESLOPE_PMB_NR     GAMMA_PMB    SIGMA[1,1] 
   43.5387124     1.6707773     0.5786627    38.8936676    36.8925865     6.4549722 

evaluate_design(poped.db2)

$ofv
[1] 38.55947

$fim
                    KNET       BMAX       INOC ESLOPE_PMB_NR  GAMMA_PMB SIGMA[1,1]
KNET           53703.462  6729.2954  2715.7124    -54735.118  71330.641      0.000
BMAX            6729.295   996.8936   376.4853     -6851.225  10780.754      0.000
INOC            2715.712   376.4853   171.5168     -2788.706   3319.721      0.000
ESLOPE_PMB_NR -54735.118 -6851.2252 -2788.7059     55823.626 -70947.690      0.000
GAMMA_PMB      71330.641 10780.7541  3319.7208    -70947.690 203222.874      0.000
SIGMA[1,1]         0.000     0.0000     0.0000         0.000      0.000   5130.072

$rse
         KNET          BMAX          INOC ESLOPE_PMB_NR     GAMMA_PMB    SIGMA[1,1] 
   103.626144      1.672481      8.053468     93.355070     69.042877      7.453560 

Full code:

library(PopED)
library(ggplot2)
library(deSolve)
set.seed(1998)

sfg <- function(x,a,bpop,b,bocc){
  parameters=c(KA = bpop[1],  #no IIV
               Vma= bpop[2],
               Ama= bpop[3],
               CL = bpop[4],
               Vme= bpop[5],
               Cm = bpop[6],
               V1 = bpop[7],
               V2 = bpop[8],
               Q  = bpop[9],
               KNET = bpop[10],
               BMAX = bpop[11],
               INOC = bpop[12],
               ESLOPE_PMB_NR = bpop[13],
               GAMMA_PMB = bpop[14],
               DOSE=a[1],
               TAU=a[2])
  return(parameters)
}

PKPD.resistant.ip.ode <- function(Time, State, Pars){
  with(as.list(c(State, Pars)), {
    CFREE = A3/V1*0.086               #free concentration, 0.086 as fraction unbound
    PLATEAU = 1 - (A5/(10**BMAX))
    KILL_PMB_NR = 0
    if(CFREE > 0) KILL_PMB_NR = (ESLOPE_PMB_NR * (CFREE**GAMMA_PMB))

    dA1 <- -Vma*A1/(A1+Ama)
    dA2 <- -KA*A2
    dA3 <- Vma*A1/(A1 + Ama) + KA*A2 -
      Q/V1*A3 + Q/V2*A4 -
      Vme*A3/(A3+(Cm*V1)) - CL/V1*A3
    dA4 <- Q/V1*A3 - Q/V2*A4
    dA5 <- (KNET*PLATEAU - KILL_PMB_NR) * A5
    return(list(c(dA1, dA2, dA3, dA4, dA5)))
  })
}

ff.PKPD.resistant.ip.ode <- function(model_switch, xt, parameters, poped.db){
  with(as.list(parameters),{
    # initial conditions of ODE system
    A_ini <- c(A1=0, A2=0, A3=0, A4=0, A5=(10**INOC))
    #Set up time points to get ODE solutions
    times_xt <- drop(xt) # sample times
    times_start <- c(0) # add extra time for start of study
    times_dose = seq(from=0,to=max(times_xt),by=TAU) # dose times
    times <- unique(sort(c(times_start,times_xt,times_dose))) # combine it all


    # Dosing
    dose_set <- c()
    for (i in DOSE)
    {
      if (i == 0)
      {
        dose_set <- append(dose_set, 0)
        dose_set <- append(dose_set, 0)
      } else

        if (i > 0 & i < 4.8) {
          dose_set <- append(dose_set, i)
          dose_set <- append(dose_set, 0)
        } else {
          dose_set <- append(dose_set, 4.8)
          dose_set <- append(dose_set, i - 4.8)
        }
    }
    dose_dat <- data.frame(
      var = c("A1","A2"),
      time = rep(times_dose, each = 2),  #dose administration in two compartment (A1 and A2) at the same time
      value = c(dose_set),
      method = c("add")
    )
    out <- ode(A_ini, times, PKPD.resistant.ip.ode, parameters,
               events = list(data = dose_dat))#atol=1e-13,rtol=1e-13)
    pd <- out[,"A5"]
    y = log10(pd)
    # extract the time points of the observations
    y = y[match(times_xt,out[,"time"])]
    y = cbind(y)
    return(list(y=y,poped.db=poped.db))
  })
}

feps_ODE_compiled <- function(model_switch,xt,parameters,epsi,poped.db){

  returnArgs <- ff.PKPD.resistant.ip.ode(model_switch,xt,parameters,poped.db)
  y <- returnArgs[[1]]
  poped.db <- returnArgs[[2]]
  y = y + epsi[,1]

  return(list(y=y,poped.db=poped.db))
}

poped.db <- create.poped.database(ff_file="ff.PKPD.resistant.ip.ode",
                                  fg_file="sfg",
                                  fError_file="feps_ODE_compiled",
                                  bpop=c(KA = 1.83,
                                         Vma= 11.3,
                                         Ama= 11.5,
                                         CL = 0.178,
                                         Vme= 0.255,
                                         Cm = 0.823,
                                         V1 = 0.325,
                                         V2 = 0.808,
                                         Q  = 0.198,
                                         KNET = 1.11,
                                         BMAX = 7.453,
                                         INOC = 6.81,
                                         ESLOPE_PMB_NR = 1.216,
                                         GAMMA_PMB = 0.02626),
                                  notfixed_bpop=c(rep(0, 9),1,1,1,1,1),
                                  sigma=0.4328^2,
                                  #design
                                  groupsize = 4,
                                  m = 30,
                                  a =list(c(DOSE=0, TAU=24),
                                          c(DOSE=0.5, TAU=24), c(DOSE=10, TAU=24),
                                          c(DOSE=20, TAU=24), c(DOSE=30, TAU=24), c(DOSE=45, TAU=24),
                                          c(DOSE=5, TAU=12), c(DOSE=10, TAU=12), c(DOSE=15, TAU=12),
                                          c(DOSE=22.5, TAU=12), c(DOSE=30, TAU=12), c(DOSE=45, TAU=12),
                                          c(DOSE=0.83, TAU=8), c(DOSE=1.67, TAU=8), c(DOSE=3.33, TAU=8),
                                          c(DOSE=5, TAU=8), c(DOSE=6.67, TAU=8), c(DOSE=7.5, TAU=8),
                                          c(DOSE=10, TAU=8), c(DOSE=15, TAU=8), c(DOSE=20, TAU=8),
                                          c(DOSE=25, TAU=8), c(DOSE=30, TAU=8), c(DOSE=40, TAU=8),
                                          c(DOSE=1.67, TAU=4), c(DOSE=3.33, TAU=4), c(DOSE=5, TAU=4),
                                          c(DOSE=7.5, TAU=4), c(DOSE=10, TAU=4), c(DOSE=15, TAU=4)),
                                  mina = c(DOSE=0,TAU=4),
                                  maxa = c(DOSE=45, TAU=24),
                                  ourzero = 0,
                                  #Design Space
                                  xt = c(0, 10, 20, 22.5), discrete_xt = list(0, seq(0, 24, 0.5),seq(0, 24, 0.5),22.5),
                                  bUseGrouped_xt=TRUE)

poped.db2 <- create.poped.database(ff_file="ff.PKPD.resistant.ip.ode",
                                   fg_file="sfg",
                                   fError_file="feps_ODE_compiled",
                                   bpop=c(KA = 1.83,
                                          Vma= 11.3,
                                          Ama= 11.5,
                                          CL = 0.178,
                                          Vme= 0.255,
                                          Cm = 0.823,
                                          V1 = 0.325,
                                          V2 = 0.808,
                                          Q  = 0.198,
                                          KNET = 1.11,
                                          BMAX = 7.453,
                                          INOC = 6.81,
                                          ESLOPE_PMB_NR = 1.216,
                                          GAMMA_PMB = 0.02626),
                                   notfixed_bpop=c(rep(0, 9),1,1,1,1,1),
                                   sigma=0.4328^2,
                                   #design
                                   groupsize = 4,
                                   m = 30,
                                   a =list(c(DOSE=0, TAU=24),
                                           c(DOSE=0.5, TAU=24), c(DOSE=10, TAU=24),
                                           c(DOSE=20, TAU=24), c(DOSE=30, TAU=24), c(DOSE=45, TAU=24),
                                           c(DOSE=5, TAU=12), c(DOSE=10, TAU=12), c(DOSE=15, TAU=12),
                                           c(DOSE=22.5, TAU=12), c(DOSE=30, TAU=12), c(DOSE=45, TAU=12),
                                           c(DOSE=0.83, TAU=8), c(DOSE=1.67, TAU=8), c(DOSE=3.33, TAU=8),
                                           c(DOSE=5, TAU=8), c(DOSE=6.67, TAU=8), c(DOSE=7.5, TAU=8),
                                           c(DOSE=10, TAU=8), c(DOSE=15, TAU=8), c(DOSE=20, TAU=8),
                                           c(DOSE=25, TAU=8), c(DOSE=30, TAU=8), c(DOSE=40, TAU=8),
                                           c(DOSE=1.67, TAU=4), c(DOSE=3.33, TAU=4), c(DOSE=5, TAU=4),
                                           c(DOSE=7.5, TAU=4), c(DOSE=10, TAU=4), c(DOSE=15, TAU=4)),
                                   mina = c(DOSE=0,TAU=4),
                                   maxa = c(DOSE=45, TAU=24),
                                   #Design Space
                                   xt = c(10, 20, 22.5), discrete_xt = list(seq(0, 24, 0.5), seq(0, 24, 0.5), seq(0, 24, 0.5)),
                                   bUseGrouped_xt=TRUE)

eval_1 <- evaluate_design(poped.db)

eval_2 <- evaluate_design(poped.db2)

from poped.

andrewhooker avatar andrewhooker commented on June 3, 2024

Hi Vincent,

I agree this is a surprising result and seems to have something to do with the "ourzero" functionality. I'll take a look and get back to you.

Andy

from poped.

Related Issues (20)

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.