search

andrewhooker / PopED

Population Experimental Design (PopED) in R
https://andrewhooker.github.io/PopED/
GNU Lesser General Public License v3.0
31 stars 21 forks source link

Illogical OFV when xts have a lot of significant digits #57

Open Vincent-AC opened 3 years ago

Vincent-AC commented 3 years ago

Hello,

I am having an issue with the following code. I have two poped.db that differ only by 1 sampling time (xt). poped.db one having an additional sample at t=10^-5. poped.db$design$xt

       obs_1    obs_2    obs_3    obs_4
grp_1  1e-05 10.00001 20.00001 22.50001
grp_2  1e-05 10.00001 20.00001 22.50001
grp_3  1e-05 10.00001 20.00001 22.50001
...

poped.db2$design$xt

          obs_1    obs_2    obs_3
grp_1  10.00001 20.00001 22.50001
grp_2  10.00001 20.00001 22.50001
grp_3  10.00001 20.00001 22.50001
...

So in theory, the ofv of poped.db should be equal or better than the ofv of poped.db2 but in practice : evaluate_design(poped.db)

$ofv
[1] 43.91572

$fim
                    KNET       BMAX       INOC ESLOPE_PMB_NR  GAMMA_PMB SIGMA[1,1]
KNET           53703.510  6729.2907  2715.7226    -54735.168  71329.543      0.000
BMAX            6729.291   984.2588   390.4884     -6851.221  10780.391      0.000
INOC            2715.723   390.4884   828.8236     -2788.715   3319.783      0.000
ESLOPE_PMB_NR -54735.168 -6851.2210 -2788.7153     55823.678 -70946.687      0.000
GAMMA_PMB      71329.543 10780.3908  3319.7827    -70946.687 203204.863      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] 
   45.5228679     1.8679539     0.5727989    40.7431924    39.1848724     6.4549722

evaluate_design(poped.db2)

$ofv
[1] 67.28574

$fim
                     KNET       BMAX         INOC ESLOPE_PMB_NR  GAMMA_PMB SIGMA[1,1]
KNET           363225.581  9204.1376 5.843678e+06     -27047.92  71330.730      0.000
BMAX             9204.138   995.3279 4.708744e+04      -2952.73  10780.763      0.000
INOC          5843678.039 47087.4356 1.102244e+08     521117.71   3319.722      0.000
ESLOPE_PMB_NR  -27047.917 -2952.7295 5.211177e+05     691485.22 -70947.776      0.000
GAMMA_PMB       71330.730 10780.7627 3.319722e+03     -70947.78 203223.119      0.000
SIGMA[1,1]          0.000     0.0000 0.000000e+00          0.00      0.000   5130.072

$rse
         KNET          BMAX          INOC ESLOPE_PMB_NR     GAMMA_PMB    SIGMA[1,1] 
   1.22385846    1.49656707    0.01010109    0.11616444   13.61248583    7.45355992

Could there be a problem with comptation of the fim when xt have a lot of significant digits ?

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),

                                  #Design Space
                                  xt = c(10^-5, 10+10^-5, 20+10^-5, 22.5+10^-5), discrete_xt = list(10^-5, seq(0.5+10^-5, 24+10^-5, 0.5), seq(0.5+10^-5, 24+10^-5, 0.5),22.5+10^-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+10^-5, 20+10^-5, 22.5+10^-5), discrete_xt = list(seq(10^-5, 24+10^-5, 0.5), seq(10^-5, 24+10^-5, 0.5), seq(10^-5, 24+10^-5, 0.5)),
                                   bUseGrouped_xt=TRUE)

eval_1 <- evaluate_design(poped.db)
poped.db$design$xt
eval_2 <- evaluate_design(poped.db2)
poped.db2$design$xt
Vincent-AC commented 3 years ago

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)
andrewhooker commented 3 years ago

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