Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

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

Open
Vincent-AC opened this issue May 11, 2021 · 2 comments
Open

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

Vincent-AC opened this issue May 11, 2021 · 2 comments
Assignees
@andrewhooker

Comments

@Vincent-AC
Copy link
Contributor

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
Copy link
Contributor Author

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
Copy link
Owner

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

@andrewhooker andrewhooker self-assigned this May 11, 2021
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees

@andrewhooker andrewhooker

Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@andrewhooker @Vincent-AC