More Advanced Feature Request: Only zero-protect values that need it

[billdenney]

I just generated a model using the following code, and it unnecessarily zero-protected V in the $DES block. It could have known that V is positive because it's exponentiated in the $PK block.

one.cmt <- function() {
    ini({
      tka <- 0.45 ; label("Ka")
      tcl <- log(c(0, 2.7, 100)) ; label("Log Cl")
      tv <- 3.45; label("log V")
      eta.ka ~ 0.6
      eta.cl ~ 0.3
      eta.v ~ 0.1
      add.sd <- 0.7
    })
    model({
      ka <- exp(tka + eta.ka)
      cl <- exp(tcl + eta.cl)
      v <- exp(tv + eta.v)
      d/dt(depot) <- -ka * depot
      d/dt(central) <- ka * depot - cl/v * central
      cp <- central / v
      cp ~ add(add.sd)
    })
  }
  
  ui <- nlmixr(one.cmt, data=nlmixr2data::Oral_1CPT, est="nonmem")

The generated NONMEM model has the following where the RXDZ001 protection of V is unnecessary.

$PK
  MU_1=THETA(1)
  MU_2=THETA(2)
  MU_3=THETA(3)
  KA=DEXP(MU_1+ETA(1)) ; ka <- exp(tka)
  CL=DEXP(MU_2+ETA(2)) ; cl <- exp(tcl)
  V=DEXP(MU_3+ETA(3)) ; v <- exp(tv)

$DES
  DADT(1) = - KA*A(1) ; d/dt(depot) = -ka * depot
  RXDZ001=V
  IF (RXDZ001 .GE. 0.0 .AND. RXDZ001 .LE. 0.000001) THEN
    RXDZ001=0.000001
  END IF
  IF (RXDZ001 .GE. -0.000001 .AND. RXDZ001 .LT. 0.) THEN
    RXDZ001= -0.000001
  END IF
  DADT(2) = KA*A(1)-CL/RXDZ001*A(2) ; d/dt(central) = ka * depot - cl/v * central
  CP=A(2)/RXDZ001 ; cp = central/v

[mattfidler]

If you want to tackle this, that is fine with me.

For now I simply want a working model even if it over-protects at first.

Once I get that done, I will tackle the more advanced features.

[billdenney]

I thought it would be a later thing (and I thought I'd take it on).

[mattfidler]

After thinking about it, this case is actually already parsed. It wouldn't be too much to add, though the general solution will take more.

[billdenney]

Yeah, I wasn't thinking it was terribly hard (but I wasn't sure if I was thinking of everything).

[mattfidler]

I might take on the mu-referenced variables and let you take on the general case; depends on who has time first 😄

[billdenney]

It's a race!

[mattfidler]

It now works for me. If you want to check the code or solution see #40

However, it is not totally general.

Here is the control stream:

$PROBLEM  translated from babelmixr2
; comments show mu referenced model in ui$getSplitMuModel

$DATA one.cmt.csv IGNORE=@

$INPUT ID TIME EVID AMT DV CMT RXROW

$SUBROUTINES ADVAN13 TOL=6

$MODEL NCOMPARTMENTS=2
     COMP(DEPOT, DEFDOSE) ; depot
     COMP(CENTRAL) ; central

$PK
  MU_1=THETA(1)
  MU_2=THETA(2)
  MU_3=THETA(3)
  KA=DEXP(MU_1+ETA(1)) ; ka <- exp(tka)
  CL=DEXP(MU_2+ETA(2)) ; cl <- exp(tcl)
  V=DEXP(MU_3+ETA(3)) ; v <- exp(tv)

$DES
  DADT(1) = - KA*A(1) ; d/dt(depot) = -ka * depot
  DADT(2) = KA*A(1)-CL/V*A(2) ; d/dt(central) = ka * depot - cl/v * central
  CP=A(2)/V ; cp = central/v

$ERROR
  ;Redefine LHS in $DES by prefixing with on RXE_ for $ERROR
  RXE_CP=A(2)/V ; cp = central/v
  RX_PF1=RXE_CP ; rx_pf1 ~ cp
  ; Write out expressions for ipred and w
  RX_IP1 = RX_PF1
  RX_P1 = RX_IP1
  W1=DSQRT((THETA(4))**2) ; W1 ~ sqrt((add.sd)^2)
  IF (W1 .EQ. 0.0) W1 = 1
  IPRED = RX_IP1
  W     = W1
  Y     = IPRED + W*EPS(1)

$THETA (0.45                 ) ; 1 - tka   
       (-INF, 0.99325, 4.6052) ; 2 - tcl   
       (3.45                 ) ; 3 - tv    
       (0,    0.7            ) ; 4 - add.sd

$OMEGA 0.6 ; eta.ka
$OMEGA 0.3 ; eta.cl
$OMEGA 0.1 ; eta.v

$SIGMA 1 FIX

$ESTIMATION METHOD=1 INTER MAXEVALS=10000 SIGDIG=3 SIGL=12 PRINT=1 NOABORT

$COVARIANCE

$TABLE ID ETAS(1:LAST) OBJI FIRSTONLY ONEHEADER NOPRINT
    NOAPPEND FILE=one.cmt.eta

$TABLE ID TIME IPRED PRED RXROW ONEHEADER NOPRINT
    NOAPPEND FILE=one.cmt.pred

[billdenney]

You won that race, by far!

[mattfidler]

Well it isn't quite complete, I just did the "low hanging fruit" so to speak.

Consider the model:

one.cmt <- function() {
    ini({
      tka <- 0.45 ; label("Ka")
      tcl <- log(c(0, 2.7, 100)) ; label("Log Cl")
      tv <- 3.45; label("log V")
      eta.ka ~ 0.6
      eta.cl ~ 0.3
      eta.v ~ 0.1
      add.sd <- 0.7
    })
    model({
      ka <- exp(tka + eta.ka)
      cl <- exp(tcl + eta.cl)
      v <- exp(tv + eta.v)
      b1 <- exp(ka+cl+v)
      d/dt(depot) <- -ka * depot
      d/dt(central) <- ka * depot - cl/v * central
      cp <- central / v
      nonsense <- central / b1
      cp ~ add(add.sd)
    })
  }

While it is nonsense, it will protect b1, though you can assume from the code that it in fact positive:

$PROBLEM  translated from babelmixr2
; comments show mu referenced model in ui$getSplitMuModel

$DATA one.cmt.csv IGNORE=@

$INPUT ID TIME EVID AMT DV CMT RXROW

$SUBROUTINES ADVAN13 TOL=6

$MODEL NCOMPARTMENTS=2
     COMP(DEPOT, DEFDOSE) ; depot
     COMP(CENTRAL) ; central

$PK
  MU_1=THETA(1)
  MU_2=THETA(2)
  MU_3=THETA(3)
  KA=DEXP(MU_1+ETA(1)) ; ka <- exp(tka)
  CL=DEXP(MU_2+ETA(2)) ; cl <- exp(tcl)
  V=DEXP(MU_3+ETA(3)) ; v <- exp(tv)

$DES
  RXR1=DEXP(KA+CL+V) ; b1 = exp(ka + cl + v)
  DADT(1) = - KA*A(1) ; d/dt(depot) = -ka * depot
  DADT(2) = KA*A(1)-CL/V*A(2) ; d/dt(central) = ka * depot - cl/v * central
  CP=A(2)/V ; cp = central/v
  RXDZ001=RXR1
  IF (RXDZ001 .GE. 0.0 .AND. RXDZ001 .LE. 0.000001) THEN
    RXDZ001=0.000001
  END IF
  IF (RXDZ001 .GE. -0.000001 .AND. RXDZ001 .LT. 0.) THEN
    RXDZ001= -0.000001
  END IF
  NONSENSE=A(2)/RXDZ001 ; nonsense = central/b1

$ERROR
  ;Redefine LHS in $DES by prefixing with on RXE_ for $ERROR
  RXE_RXR1=DEXP(KA+CL+V) ; b1 = exp(ka + cl + v)
  RXE_CP=A(2)/V ; cp = central/v
  RXDZE001=RXE_RXR1
  IF (RXDZE001 .GE. 0.0 .AND. RXDZE001 .LE. 0.000001) THEN
    RXDZE001=0.000001
  END IF
  IF (RXDZE001 .GE. -0.000001 .AND. RXDZE001 .LT. 0.) THEN
    RXDZE001= -0.000001
  END IF
  RXE_NONSENSE=A(2)/RXDZE001 ; nonsense = central/b1
  RX_PF1=RXE_CP ; rx_pf1 ~ cp
  ; Write out expressions for ipred and w
  RX_IP1 = RX_PF1
  RX_P1 = RX_IP1
  W1=DSQRT((THETA(4))**2) ; W1 ~ sqrt((add.sd)^2)
  IF (W1 .EQ. 0.0) W1 = 1
  IPRED = RX_IP1
  W     = W1
  Y     = IPRED + W*EPS(1)

$THETA (0.45                 ) ; 1 - tka   
       (-INF, 0.99325, 4.6052) ; 2 - tcl   
       (3.45                 ) ; 3 - tv    
       (0,    0.7            ) ; 4 - add.sd

$OMEGA 0.6 ; eta.ka
$OMEGA 0.3 ; eta.cl
$OMEGA 0.1 ; eta.v

$SIGMA 1 FIX

$ESTIMATION METHOD=1 INTER MAXEVALS=10000 SIGDIG=3 SIGL=12 PRINT=1 NOABORT

$COVARIANCE

$TABLE ID ETAS(1:LAST) OBJI FIRSTONLY ONEHEADER NOPRINT
    NOAPPEND FILE=one.cmt.eta

$TABLE ID TIME IPRED PRED RXROW ONEHEADER NOPRINT
    NOAPPEND FILE=one.cmt.pred
```

[billdenney]

Another part to this (and maybe another low-hanging fruit): do not zero-protect dividing by a number. In a model that I'm testing now, WT/70 had the 70 zero-protected.

[mattfidler]

I think I saw this. I justified this because rxode2 may do this since its parsing isn't quite as smart here. (Worth mentioning since it may change the correlation between the nonmem output and the nlmixr2 output and may be a future feature request. It is a bit harder to implement though since it is entirely in C)