individual model with a discreate covariate

dzianismr commented 3 years ago

I get to test new nlmixr (nlmixr_2.0.4, RxODE_1.0.8 ), it is great that it is now on CRAN! The goal is to fit a simple individual PK model with SPECIE as a discreate covariate for Cl. Following discussion in 474 (https://github.com/nlmixrdevelopment/nlmixr/issues/474) I set up it the following way: ini({ tClrat <- log(3) tCldog <- log(1) ... }) model({ if ( SPECIE == 'Rat'){ ## SPECIE is a column in the data table Cl = exp(tClrat) } if ( SPECIE == 'Dog'){ Cl = exp(tCldog) } ...
}) No random effects were defined and focei, as suggest in discusion 409 (https://github.com/nlmixrdevelopment/nlmixr/issues/409) was chosen as algorithm. This crashed R w/o any error message after. The last output in console: |-----+---------------+-----------+-----------+-----------+-----------+-----------| | #| Objective Fun | tVc | tVp | tCl | tClD | prop.err |

Than I tried to introduce a dummy random effect ( eta.Vc ~ fix(0.0001) ) and use SAEM algorithm. In this case I got the following error message: warning: 'Cl' may be used uninitialized in this function [-Wmaybe-uninitialized] DDtStateVar[0] = ((double)(_ON[0]))(_IR[0] -Cl/safe_zero(Vc)center-ClD/safe_zero(Vc)center+ClD/safe_zero(Vp)periph); Error in configsaem(model = model, data = dat, inits = inits, mcmc = .mcmc, : nphi and covstruct dim mismatch

Am I using if / else in an incorrect way?

billdenney commented 3 years ago

I'm guessing that this is a duplicate of nlmixrdevelopment/RxODE#410. Can you try changing SPECIE to a numeric covariate and see if it still causes issues?

dzianismr commented 3 years ago

Thanks! changing SPECIE to numeric fixes the problem and indeed it seems to be duplicate of nlmixrdevelopment/RxODE#410.

mattfidler commented 3 years ago

Thanks @billdenney for pointing this out; Thank you @dzianismr for trying out nlmixr

dzianismr commented 3 years ago

@mattfidler I really like that nlmixr is fully R integrated @billdenney I have noticed another two issues with covariate and if/else statements: a) it seem that if/else cannot preceed ODE definition e.g. snippet below runs just fine: if ( SPECIEN == 1){ Cl = exp(tClrat) } ClD = exp(tClD) d/dt(center) = ... where is this ClD = exp(tClD) if ( SPECIEN == 1){ Cl = exp(tClrat) } d/dt(center) = ... crushes with Error in parse(text = paste(c("function() {", rx.txt[1:w], "}"), collapse = "\n")) :

:17:0: unexpected end of input 15: Cl = exp(tClrat) 16: } ^ b) model output does not include covariate column SPECIEN: m1 <- nlmixr( two.cpt.cov, dat, "focei",list(print=10)) if i try to force it with the tableControl() m1 <- nlmixr( two.cpt.cov dat, "focei",list(print=10), table = tableControl( covariates = TRUE ) ) it does not seems to have any effect How do I get covariates in the output object m1?

billdenney commented 3 years ago

@dzianismr , that one is out of my depth. Can you please provide a full reproducible example? (with data that can be blinded and a minimal model that causes the problem)

mattfidler commented 3 years ago

In this case a simple model may be sufficient; Let me know when that happens so I can help you with your problem.

dzianismr commented 3 years ago

I have produced a simple example

generate test data

m2c_rx <- function() {

Vc  = exp(tVc )  
Vp  = exp(tVp)

if ( SPECIEN == 2){
  Cl  = exp(tCl2) 
}
if ( SPECIEN == 1){
  Cl  = exp(tCl1) 
} 
ClD = exp(tClD) 

d/dt(center) = - Cl / Vc * center - ClD / Vc * center + ClD / Vp * periph 
d/dt(periph) =                      ClD / Vc * center - ClD / Vp * periph
cp = center / Vc

}

ev <- eventTable() %>% et(amt = 1, dur = 1, evid = 1, cmt = 1 ) %>% et(0.25, 3, by=0.25) %>% et(id = 1:2) ev$SPECIEN <- ev$id

m2rx <- RxODE(model = m2c_rx ) theta <- c( tVc = log(0.2), tVp = log(0.5), tClD = log(1), tCl1 = log(0.5), tCl2 = log(2) )

df_sim <- rxSolve( m2rx, theta, ev )

dat <- data.frame( time = df_sim$time, dv = df_sim$cp * ( 1 + rnorm( length(df_sim$cp), 0, sd = 0.1) ), id = df_sim$id, SPECIEN = df_sim$id, amt = 0, evid = 0, dur = 0, cmt = 1 ) tmp1 <- data.frame( time = 0, dv = NA, id = 1, SPECIEN = 1, amt = 1, evid = 1, dur = 1, cmt = 1 ) dat <- rbind( dat, tmp1 ) tmp1[,c("SPECIEN", "id") ] <-2 dat <- rbind( dat, tmp1 ) dat <- dat[ order(dat$id, dat$time), ]

fit

m2c_cov <- function() { ini({ tVc <- log(0.3) tVp <- log(0.3) tCl1 <- log(0.75) tCl2 <- log(1.5) tClD <- log(1) prop.err <- 0.05 }) model({ Vc = exp(tVc )
Vp = exp(tVp)

if ( SPECIEN == 2){
  Cl  = exp(tCl2) 
}
if ( SPECIEN == 1){
  Cl  = exp(tCl1) 
}
ClD = exp(tClD) # crash if you move this line above if / else statement 

d/dt(center) = - Cl / Vc * center - ClD / Vc * center + ClD / Vp * periph 
d/dt(periph) =                      ClD / Vc * center - ClD / Vp * periph
cp = center / Vc
 cp ~ prop(prop.err)

}) }

mfit <- nlmixr( m2c_cov, dat, "focei", list(print=10) )

comment

the code above runs just fine, however if you move ClD = exp(tClD) above if / else statemnt in nlmixr model definition an error will be thrown: Error in parse(text = paste(c("function() {", rx.txt[1:w], "}"), collapse = "\n")) : :11:0: unexpected end of input 9: Cl = exp(tCl1) 10: } ^

Note that in RxODE model sequence of parameter definition plays no role, so it seem sto be nlmixr problem.

mattfidler commented 3 years ago

I can reproduce the issue:

library(nlmixr)

rxClean()

##### generate test data
m2c_rx <- function() {

  Vc  = exp(tVc )
  Vp  = exp(tVp)

  if ( SPECIEN == 2){
    Cl  = exp(tCl2)
  }
  if ( SPECIEN == 1){
    Cl  = exp(tCl1)
  }
  ClD = exp(tClD)

  d/dt(center) = - Cl / Vc * center - ClD / Vc * center + ClD / Vp * periph
  d/dt(periph) =                      ClD / Vc * center - ClD / Vp * periph
  cp = center / Vc
}

ev <-  eventTable() %>%
  et(amt = 1,    dur = 1,  evid = 1, cmt = 1 ) %>%
  et(0.25, 3, by=0.25) %>%
  et(id = 1:2)
ev$SPECIEN <- ev$id

m2rx    <- RxODE(model = m2c_rx  )
#> → creating RxODE include directory
#> → getting R compile options
#> → precompiling headers
#> ✔ done
theta   <- c( tVc  = log(0.2),   tVp  = log(0.5),  tClD = log(1),
              tCl1 = log(0.5),   tCl2 = log(2)     )

df_sim  <- rxSolve( m2rx, theta, ev  )

dat     <- data.frame(
                   time = df_sim$time,
                   dv   = df_sim$cp * ( 1 + rnorm( length(df_sim$cp), 0, sd = 0.1) ),
                   id   = df_sim$id,
                   SPECIEN = df_sim$id,
                   amt  = 0,
                   evid = 0,
                   dur  = 0,
                   cmt  = 1
)
tmp1    <- data.frame(
                   time = 0,
                   dv   = NA,
                   id   = 1,
                   SPECIEN = 1,
                   amt  = 1,
                   evid = 1,
                   dur  = 1,
                   cmt  = 1
)
dat <- rbind( dat, tmp1 )
tmp1[,c("SPECIEN", "id") ] <-2
dat <- rbind( dat, tmp1 )
dat <- dat[ order(dat$id, dat$time), ]

### fit
m2c_cov <- function() {
  ini({
    tVc    <- log(0.3)
    tVp    <- log(0.3)
    tCl1   <- log(0.75)
    tCl2   <- log(1.5)
    tClD   <- log(1)
    prop.err <- 0.05
  })
  model({
    Vc  = exp(tVc )
    Vp  = exp(tVp)

    if ( SPECIEN == 2){
      Cl  = exp(tCl2)
    }
    if ( SPECIEN == 1){
      Cl  = exp(tCl1)
    }
    ClD = exp(tClD) # crash if you move this line above if / else statement

    d/dt(center) = - Cl / Vc * center - ClD / Vc * center + ClD / Vp * periph
    d/dt(periph) =                      ClD / Vc * center - ClD / Vp * periph
    cp = center / Vc
     cp ~ prop(prop.err)
  })
}

mfit      <- nlmixr( m2c_cov, dat, "focei", list(print=10) )
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling EBE model...
#> ✔ done
#> Needed Covariates:
#> [1] "SPECIEN"
#> RxODE 1.0.8.1 using 4 threads (see ?getRxThreads)
#> Key: U: Unscaled Parameters; X: Back-transformed parameters; G: Gill difference gradient approximation
#> F: Forward difference gradient approximation
#> C: Central difference gradient approximation
#> M: Mixed forward and central difference gradient approximation
#> Unscaled parameters for Omegas=chol(solve(omega));
#> Diagonals are transformed, as specified by foceiControl(diagXform=)
#> |-----+---------------+-----------+-----------+-----------+-----------|
#> |    #| Objective Fun |       tVc |       tVp |      tCl1 |      tCl2 |
#> |.....................|      tClD |  prop.err |...........|...........|
#> |-----+---------------+-----------+-----------+-----------+-----------|
#> |   10|     36.181183 |   -0.8691 |   -0.6028 |   -0.4601 |     1.281 |
#> |.....................|    0.3727 |    0.8735 |...........|...........|
#> |    U|     36.181183 |    -1.073 |   -0.8067 |   -0.8865 |    0.6865 |
#> |.....................|   -0.1235 |   0.05788 |...........|...........|
#> |    X|     36.181183 |    0.3420 |    0.4463 |    0.4121 |     1.987 |
#> |.....................|    0.8838 |   0.05788 |...........|...........|
#> |    F| Forward Diff. |     20.22 |    -498.3 |    -639.4 |     934.1 |
#> |.....................|     149.8 |    -144.7 |...........|...........|
#> |-----+---------------+-----------+-----------+-----------+-----------|
#> |   20|    -139.48405 |    -1.161 |   -0.5655 |   -0.2716 |     1.279 |
#> |.....................|    0.3758 |     1.429 |...........|...........|
#> |    U|    -139.48405 |    -1.365 |   -0.7694 |   -0.6979 |    0.6842 |
#> |.....................|   -0.1204 |   0.07178 |...........|...........|
#> |    X|    -139.48405 |    0.2554 |    0.4633 |    0.4976 |     1.982 |
#> |.....................|    0.8866 |   0.07178 |...........|...........|
#> |    F| Forward Diff. |     3.729 |    -36.55 |    -4.969 |     45.08 |
#> |.....................|     6.414 |    -25.21 |...........|...........|
#> |-----+---------------+-----------+-----------+-----------+-----------|
#> |   30|    -150.92514 |    -1.140 |   -0.5429 |   -0.2578 |     1.307 |
#> |.....................|    0.3994 |     1.885 |...........|...........|
#> |    U|    -150.92514 |    -1.344 |   -0.7469 |   -0.6841 |    0.7125 |
#> |.....................|  -0.09678 |   0.08318 |...........|...........|
#> |    X|    -150.92514 |    0.2608 |    0.4738 |    0.5045 |     2.039 |
#> |.....................|    0.9078 |   0.08318 |...........|...........|
#> |    F| Forward Diff. |    -9.676 |    -12.27 |    -36.00 |     1.358 |
#> |.....................|    -7.647 |    -6.350 |...........|...........|
#> |-----+---------------+-----------+-----------+-----------+-----------|
#> |   40|    -158.54368 |    -1.196 |   -0.5036 |   -0.2494 |     1.315 |
#> |.....................|    0.4660 |     3.195 |...........|...........|
#> |    U|    -158.54368 |    -1.400 |   -0.7076 |   -0.6757 |    0.7204 |
#> |.....................|  -0.03015 |    0.1159 |...........|...........|
#> |    X|    -158.54368 |    0.2465 |    0.4928 |    0.5088 |     2.055 |
#> |.....................|    0.9703 |    0.1159 |...........|...........|
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> done
#> Calculating residuals/tables
#> done
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod =
#> FALSE, : gradient problems with initial estimate and covariance; see $scaleInfo

m2c_cov <- function() {
  ini({
    tVc    <- log(0.3)
    tVp    <- log(0.3)
    tCl1   <- log(0.75)
    tCl2   <- log(1.5)
    tClD   <- log(1)
    prop.err <- 0.05
  })
  model({
    Vc  = exp(tVc )
    Vp  = exp(tVp)

   ClD = exp(tClD) # crash if you move this line above if / else statement

    if ( SPECIEN == 2){
      Cl  = exp(tCl2)
    }
    if ( SPECIEN == 1){
      Cl  = exp(tCl1)
    }

    d/dt(center) = - Cl / Vc * center - ClD / Vc * center + ClD / Vp * periph
    d/dt(periph) =                      ClD / Vc * center - ClD / Vp * periph
    cp = center / Vc
     cp ~ prop(prop.err)
  })
}

mfit      <- nlmixr( m2c_cov, dat, "focei", list(print=10) )
#> Error in parse(text = paste(c("function() {", rx.txt[1:w], "}"), collapse = "\n")): <text>:11:0: unexpected end of input
#> 9:         Cl = exp(tCl1)
#> 10: }
#>    ^

^{Created on 2021-05-04 by the reprex package (v2.0.0)}

mattfidler commented 3 years ago

Hi @dzianismr,

Since there is a work-around right now, I will work to make this model work with the new parsing routines; In theory, this will also add some more ferjl

nlmixrdevelopment / nlmixr

individual model with a discreate covariate #521

generate test data

fit

comment