issue with changing model features

[dzianismr]

It is a great functionality to change model features as described in model piping tutorial (https://nlmixrdevelopment.github.io/nlmixr/articles/modelPiping.html). I was experimenting with automating model selection process and run into trouble with this functionality. I was trying to add an additional if/else statement into the model (related to https://github.com/nlmixrdevelopment/nlmixr/issues/521). Below is an example:

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

    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
    cp ~ prop(prop.err)
  })
}

## an attempt to update a model code with an additional if/else statement
m2c_ifelse1upd <- nlmixr(m2c_ifelse1) %>% 
  model( { if ( SPECIEN == 2){ Cl  = exp(tCl2) }}) %>%
  ini(   tCl2   = log(1.5))

## model structure of the base model
nlmixr(m2c_ifelse1)
>> i parameter labels from comments will be replaced by 'label()'
__ RxODE-based ODE model ______________________________________________________________ 
-- Initialization: -------------------------------------------------------------------- 
Fixed Effects ($theta): 
       tVc        tVp       tCl1       tClD 
-1.2039728 -1.2039728 -0.2876821  0.0000000 

 Covariates or Uninitialized Parameters ($all.covs)
[1] "SPECIEN"
-- Model: ----------------------------------------------------------------------------- 
    Vc  = exp(tVc )  
    Vp  = exp(tVp)

    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
    #cp ~  add(add.err)
    cp ~ prop(prop.err) 
_______________________________________________________________________________________ 

## model structure of the updated model
 m1_ifelse1upd
>>__ RxODE-based ODE model ______________________________________________________________ 
-- Initialization: -------------------------------------------------------------------- 
Fixed Effects ($theta): 
       tVc        tVp       tClD       tCl2 
-1.2039728 -1.2039728  0.0000000  0.4054651 

 Covariates or Uninitialized Parameters ($all.covs)
[1] "SPECIEN"
-- Model: ----------------------------------------------------------------------------- 
    Vc  = exp(tVc )  
    Vp  = exp(tVp)

    if ( SPECIEN == 1){
        Cl = exp(tCl2)
    }
    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
    #cp ~  add(add.err)
    cp ~ prop(prop.err) 

---

issue

my expectation was that additional parameter tCl2 will be added to the ini block and that the second if/else statement will be added to the model:

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

However that is not the outcome. Am I using this change model syntax is incorrectly or do I have wrong expectations about this functionality?

[mattfidler]

Hi @dzianismr ,

Currently model piping can change one distinct line in model code, as defined by var <- .... That means this sort of if and else cannot be changed at the moment.

[dzianismr]

Many thanks @mattfidler for the prompt response.

[dzianismr]

A realted question:

Say we have a 2cmt model: m2c <- function() { ini({ tVc <- log(0.3 ) tVp <- log(0.3 ) tCl <- log(1) tClD <- log(1 ) prop.err <- 0.05 }) model({ Vc = exp(tVc )
Vp = exp(tVp ) Cl = exp(tCl ) 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
#cp ~  add(add.err)
cp ~ prop(prop.err)

}) }

The following line then neatly adds bsv for CL parameter: nlmixr(m2c) %>% model({Cl <- exp(tCl + eta.Cl)}) %>% ini(eta.Cl=1)

RxODE-based ODE model ____ -- Initialization: -------------------------------------------------------------------- Fixed Effects ($theta): tVc tVp tCl tClD -1.203973 -1.203973 0.000000 0.000000

Omega ($omega): eta.Cl eta.Cl 3 -- mu-referencing ($muRefTable): ------------------------------------------------------ theta eta 1 tCl eta.Cl

-- Model: ----------------------------------------------------------------------------- Vc = exp(tVc )
Vp = exp(tVp ) Cl <- exp(tCl + eta.Cl) 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
cp ~ prop(prop.err) 

---

However similar as in Example 2 of model piping: one.ka.0.5 <- fit %>% ini(tka=fix(0.5)) %>% nlmixr(est="focei", control=list(print=0), table=list(cwres=TRUE, npde=TRUE)) I request eta to be fixed, however exactly the same output is generated: nlmixr(m2c) %>% model({Cl <- exp(tCl + eta.Cl)}) %>% ini(eta.Cl=fix(2))

RxODE-based ODE model ____ -- Initialization: -------------------------------------------------------------------- Fixed Effects ($theta): tVc tVp tCl tClD -1.203973 -1.203973 0.000000 0.000000

Omega ($omega): eta.Cl eta.Cl 2 -- mu-referencing ($muRefTable): ------------------------------------------------------ theta eta 1 tCl eta.Cl

-- Model: ----------------------------------------------------------------------------- Vc = exp(tVc )
Vp = exp(tVp ) Cl <- exp(tCl + eta.Cl) 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
cp ~ prop(prop.err) 

---

Am I missing something?

[dzianismr]

I did a few more tests and it seems that the issue is only with the way model is printed. m2c_pop2_a<- nlmixr(m2c) %>% model({Cl <- exp(tCl + eta.Cl)}) %>% ini(eta.Cl=fix(2)) m1_pop2_a <- nlmixr(m2c_pop2_a, data = dat, est="focei", control=list(print=0), table=list(cwres=TRUE, npde=TRUE)) gives the same output as m1 <- nlmixr( m2c, dat, "focei", list(print=10) ) m1_pop2_b <- m1 %>% model({Cl <- exp(tCl + eta.Cl)}) %>% ini(eta.Cl=fix(2)) %>% nlmixr(est="focei", control=list(print=0), table=list(cwres=TRUE, npde=TRUE)) and has statement of fixed bsv: Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)% tVc -1.62 0.119 7.36 0.198 (0.157, 0.25)
tVp -0.628 0.0484 7.71 0.534 (0.485, 0.587)
tCl -0.00723 1.15 1.6e+004 0.993 (0.103, 9.54) fix(253.) 50.7%

[mattfidler]

Hi @dzianismr

I agree the problem is with the output

I am reworking the parsing engine, so this will be something to consider in the next iteration.

Thanks for working with this