`saem` non-informative `ETAs`

mattfidler commented 9 months ago

In saem when running an estimation with non-informative ETAs like a mix of oral and IV dosing, the ETAs for ka bounce around (with no consequence to saem llik), but ETAs should be zero.

In the optimization case (like focei), this doesn't happen for the individual as the gradient keeps the value equal to zero no matter what happens (since there is no information)

This was originally reported by @RichardHooijmaijers (thank you)

library(nlmixr2est)
#> Loading required package: nlmixr2data
# Simulated dataset with 2 IDs IV and 2 IDs PO
dat <- structure(list(ID = c(11, 11, 11, 11, 11, 11, 11, 12, 12, 12, 
                                 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 13, 13, 21, 21, 21, 21, 
                                 21, 21, 21, 22, 22, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 23, 23), 
                         TIME = c(0, 0.05, 0.25, 0.5, 1, 3, 5, 0, 0.05, 0.25, 0.5, 
                                     1, 3, 5, 0, 0.05, 0.25, 0.5, 1, 3, 5, 8, 0, 0.25, 0.5, 1, 3, 
                                     5, 8, 0, 0.25, 0.5, 1, 3, 5, 8, 0, 0.25, 0.5, 1, 3, 5, 8), 
                         DV = c(NA,2017.85, 1323.74, 792.5, 822.72, 36.27, 3.33, NA, 1702, 1290.75, 
                                            1095.95, 907.6, 125.44, 14.44, NA, 1933.04, 1242.43, 661.22, 
                                            193.52, 1.75, NA, NA, NA, 706.58, 1063.14, 2257.62, 941.33, 629.69, 
                                            100, NA, 1462.95, 2217.76, 2739.5, 705.3, 108.47, 8.75, NA, 211.66, 
                                            467.23, 174.24, 153.6, 27.07, 2.81), 
                         AMT = c(1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 
                                    0, 0, 5, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0), 
                         EVID = c(1,0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 
                                       1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0), 
                         CMT = c(2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 1, 2, 
                                    2, 2, 2, 2, 2), 
                         DOSE = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5), 
                         ROUTE = c(1, 1, 1, 1, 1, 1, 1,1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)), 
                     row.names = c(NA, -43L), class = "data.frame")

# First fit with if statement to indicate diff IV/PO
modA <- function() {
ini({
    tka <- 0.45   # Log Ka
    tcl <- -7     # Log Cl
    tv  <- -8     # Log V
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    prop.sd <- 0.7
    })
model({
      if(ROUTE!=1) ka <- exp(tka + eta.ka)
      cl <- exp(tcl + eta.cl)
      v <- exp(tv + eta.v)
      d/dt(depot) = -ka * depot
      d/dt(center) = ka * depot - cl / v * center
      cp = center / v
      cp ~ prop(prop.sd)
      })
}
fitA <- nlmixr(modA, dat, est="saem", control=list(print=0))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: eta.ka
#> as a work-around try putting the mu-referenced expression on a simple line

#> Warning: some etas defaulted to non-mu referenced, possible parsing error: eta.ka
#> as a work-around try putting the mu-referenced expression on a simple line
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of saem model...
#> ✔ done
#> → finding duplicate expressions in saem model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> ✔ done
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> In file included from /usr/share/R/include/R.h:71,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2parse.h:33,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2/include/rxode2.h:9,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2_model_shared.h:3,
#>                  from rx_ab84949ca37377f3e2406ed3a07a5995_.c:117:
#> /usr/share/R/include/R_ext/Complex.h:80:6: warning: ISO C99 doesn’t support unnamed structs/unions [-Wpedantic]
#>    80 |     };
#>       |      ^
#> rxode2 2.1.2 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> Calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of saem model...
#> ✔ done
#> → finding duplicate expressions in saem predOnly model 0...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in saem predOnly model 1...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem predOnly model 1...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in saem predOnly model 2...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem predOnly model 2...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> ✔ done
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> In file included from /usr/share/R/include/R.h:71,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2parse.h:33,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2/include/rxode2.h:9,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2_model_shared.h:3,
#>                  from rx_51f78b882d72698f012fc08565c5b8f9_.c:117:
#> /usr/share/R/include/R_ext/Complex.h:80:6: warning: ISO C99 doesn’t support unnamed structs/unions [-Wpedantic]
#>    80 |     };
#>       |      ^
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 3688
#> → compress phiM in nlmixr2 object, save 84688
#> → compress parHistData in nlmixr2 object, save 15192
#> → compress saem0 in nlmixr2 object, save 25000
unique(fitA$ka[fitA$ROUTE==1])     # KA is 0 as expected
#> [1] 0
unique(fitA$eta.ka[fitA$ROUTE==1]) # eta KA for IV is not 0 (would expect 0 or NA?)
#> [1]  0.02217626  0.02686848 -0.02517917
fitA$parFixed # estimate and back-transformed estimate for Ka is the same (not for other parameters)
#>         Parameter   Est.     SE %RSE       Back-transformed(95%CI)
#> tka        Log Ka -0.285  0.173 60.6       -0.285 (-0.624, 0.0534)
#> tcl        Log Cl  -6.87  0.312 4.54   0.00104 (0.000565, 0.00192)
#> tv          Log V  -7.54 0.0818 1.09 0.000531 (0.000452, 0.000623)
#> prop.sd             0.27                                      0.27
#> eta.ka                                                            
#>         BSV(CV% or SD) Shrink(SD)%
#> tka                               
#> tcl               87.6     -7.88%>
#> tv                1.16      85.3%>
#> prop.sd                           
#> eta.ka           0.284      28.2%=

# Second fit without if statement
modB <- function() {
ini({
    tka <- 0.45   # Log Ka
    tcl <- -7     # Log Cl
    tv  <- -8     # Log V
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    prop.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(center) = ka * depot - cl / v * center
      cp = center / v
      cp ~ prop(prop.sd)
      })
}
fitB <- nlmixr(modB, dat, est="saem", control=list(print=0))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of saem model...
#> ✔ done
#> → finding duplicate expressions in saem model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> ✔ done
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> In file included from /usr/share/R/include/R.h:71,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2parse.h:33,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2/include/rxode2.h:9,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2_model_shared.h:3,
#>                  from rx_8a4f34ae1b3c53747c7e0a1c83eafa3e_.c:117:
#> /usr/share/R/include/R_ext/Complex.h:80:6: warning: ISO C99 doesn’t support unnamed structs/unions [-Wpedantic]
#>    80 |     };
#>       |      ^
#> Calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of saem model...
#> ✔ done
#> → finding duplicate expressions in saem predOnly model 0...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in saem predOnly model 1...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem predOnly model 1...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in saem predOnly model 2...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem predOnly model 2...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> ✔ done
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> In file included from /usr/share/R/include/R.h:71,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2parse.h:33,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2/include/rxode2.h:9,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2_model_shared.h:3,
#>                  from rx_b162c06a99e02078337f456da254e8ed_.c:117:
#> /usr/share/R/include/R_ext/Complex.h:80:6: warning: ISO C99 doesn’t support unnamed structs/unions [-Wpedantic]
#>    80 |     };
#>       |      ^
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 3688
#> → compress phiM in nlmixr2 object, save 42976
#> → compress parHistData in nlmixr2 object, save 14376
#> → compress saem0 in nlmixr2 object, save 24296
unique(fitB$ka[fitB$ID%in%11:13])     # KA is estimated also for IV group
#> [1] 0.6471118 0.6227864 0.6450586
unique(fitB$eta.ka[fitB$ID%in%11:13]) # eta KA for IV is not 0 (would expect 0 or NA?)
#> [1]  0.02667433 -0.01164104  0.02349642
fitB$parFixed # in here the backtransformed parameters are as expected (e.g. ka is backtransformed)
#>         Parameter   Est.    SE %RSE      Back-transformed(95%CI) BSV(CV%)
#> tka        Log Ka -0.462 0.169 36.5          0.63 (0.453, 0.877)     28.2
#> tcl        Log Cl  -6.85  0.31 4.53  0.00106 (0.000575, 0.00194)     87.0
#> tv          Log V  -7.52 0.083  1.1 0.000541 (0.00046, 0.000636)     2.59
#> prop.sd            0.271                                   0.271         
#>         Shrink(SD)%
#> tka          31.0%>
#> tcl         -8.49%>
#> tv           85.7%>
#> prop.sd

^{Created on 2024-02-02 with reprex v2.1.0}

This will require a work-around that somehow identifies when ETA values are non-informative for each individual and then sets these ETA values to zero.

Ideally this would be in absence of a gradient :smile:

mattfidler commented 9 months ago

This is the same problem using focei instead; while it doesn't give exactly zero, the gradient from the inner problem makes them approximately zero.

library(nlmixr2est)
#> Loading required package: nlmixr2data
# Simulated dataset with 2 IDs IV and 2 IDs PO
dat <- structure(list(ID = c(11, 11, 11, 11, 11, 11, 11, 12, 12, 12, 
                                 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 13, 13, 21, 21, 21, 21, 
                                 21, 21, 21, 22, 22, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 23, 23), 
                         TIME = c(0, 0.05, 0.25, 0.5, 1, 3, 5, 0, 0.05, 0.25, 0.5, 
                                     1, 3, 5, 0, 0.05, 0.25, 0.5, 1, 3, 5, 8, 0, 0.25, 0.5, 1, 3, 
                                     5, 8, 0, 0.25, 0.5, 1, 3, 5, 8, 0, 0.25, 0.5, 1, 3, 5, 8), 
                         DV = c(NA,2017.85, 1323.74, 792.5, 822.72, 36.27, 3.33, NA, 1702, 1290.75, 
                                            1095.95, 907.6, 125.44, 14.44, NA, 1933.04, 1242.43, 661.22, 
                                            193.52, 1.75, NA, NA, NA, 706.58, 1063.14, 2257.62, 941.33, 629.69, 
                                            100, NA, 1462.95, 2217.76, 2739.5, 705.3, 108.47, 8.75, NA, 211.66, 
                                            467.23, 174.24, 153.6, 27.07, 2.81), 
                         AMT = c(1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 
                                    0, 0, 5, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0), 
                         EVID = c(1,0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 
                                       1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0), 
                         CMT = c(2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 1, 2, 
                                    2, 2, 2, 2, 2), 
                         DOSE = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5), 
                         ROUTE = c(1, 1, 1, 1, 1, 1, 1,1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)), 
                     row.names = c(NA, -43L), class = "data.frame")

# First fit with if statement to indicate diff IV/PO
modA <- function() {
ini({
    tka <- 0.45   # Log Ka
    tcl <- -7     # Log Cl
    tv  <- -8     # Log V
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    prop.sd <- 0.7
    })
model({
      if(ROUTE!=1) ka <- exp(tka + eta.ka)
      cl <- exp(tcl + eta.cl)
      v <- exp(tv + eta.v)
      d/dt(depot) = -ka * depot
      d/dt(center) = ka * depot - cl / v * center
      cp = center / v
      cp ~ prop(prop.sd)
      })
}
fitA <- nlmixr(modA, dat, est="focei", control=list(print=0))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: eta.ka
#> as a work-around try putting the mu-referenced expression on a simple line

#> Warning: some etas defaulted to non-mu referenced, possible parsing error: eta.ka
#> as a work-around try putting the mu-referenced expression on a simple line
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of full model...
#> ✔ done
#> → calculate jacobian
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate sensitivities
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate ∂(f)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate ∂(R²)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling inner model...
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> In file included from /usr/share/R/include/R.h:71,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2parse.h:33,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2/include/rxode2.h:9,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2_model_shared.h:3,
#>                  from rx_8c54c7bbe5eb107d1c9e3a3cd09be90d_.c:117:
#> /usr/share/R/include/R_ext/Complex.h:80:6: warning: ISO C99 doesn’t support unnamed structs/unions [-Wpedantic]
#>    80 |     };
#>       |      ^
#> ✔ done
#> → finding duplicate expressions in FD model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in FD model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling EBE model...
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> In file included from /usr/share/R/include/R.h:71,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2parse.h:33,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2/include/rxode2.h:9,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2_model_shared.h:3,
#>                  from rx_5b6653fa088810056571e149981aea5b_.c:117:
#> /usr/share/R/include/R_ext/Complex.h:80:6: warning: ISO C99 doesn’t support unnamed structs/unions [-Wpedantic]
#>    80 |     };
#>       |      ^
#> ✔ done
#> → compiling events FD model...
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> In file included from /usr/share/R/include/R.h:71,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2parse.h:33,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2/include/rxode2.h:9,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2_model_shared.h:3,
#>                  from rx_8f1f7d99ebcf008b1c97a36acb990f0f_.c:117:
#> /usr/share/R/include/R_ext/Complex.h:80:6: warning: ISO C99 doesn’t support unnamed structs/unions [-Wpedantic]
#>    80 |     };
#>       |      ^
#> ✔ done
#> rxode2 2.1.2 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 3688
#> → compress parHistData in nlmixr2 object, save 5256
unique(fitA$ka[fitA$ROUTE==1])     # KA is 0 as expected
#> [1] 0
unique(fitA$eta.ka[fitA$ROUTE==1]) # eta KA for IV is not 0 (would expect 0 or NA?)
#> [1] -5.683213e-10  5.233306e-07 -3.404403e-10
fitA$parFixed # estimate and back-transformed estimate for Ka is the same (not for other parameters)
#>         Parameter   Est.    SE %RSE      Back-transformed(95%CI) BSV(CV% or SD)
#> tka        Log Ka -0.439 0.193 43.9     -0.439 (-0.817, -0.0609)               
#> tcl        Log Cl  -6.79 0.293 4.32    0.00113 (0.000635, 0.002)           89.9
#> tv          Log V  -7.44 0.153 2.06 0.00059 (0.000437, 0.000797)           30.1
#> prop.sd            0.328                                   0.328               
#> eta.ka                                                                    0.434
#>         Shrink(SD)%
#> tka                
#> tcl         -4.89%>
#> tv           60.9%>
#> prop.sd            
#> eta.ka       60.7%>

# Second fit without if statement
modB <- function() {
ini({
    tka <- 0.45   # Log Ka
    tcl <- -7     # Log Cl
    tv  <- -8     # Log V
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    prop.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(center) = ka * depot - cl / v * center
      cp = center / v
      cp ~ prop(prop.sd)
      })
}
fitB <- nlmixr(modB, dat, est="focei", control=list(print=0))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of full model...
#> ✔ done
#> → calculate jacobian
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate sensitivities
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate ∂(f)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate ∂(R²)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling inner model...
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> In file included from /usr/share/R/include/R.h:71,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2parse.h:33,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2/include/rxode2.h:9,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2_model_shared.h:3,
#>                  from rx_befdf018ef95caabd7aaaee9b03206af_.c:117:
#> /usr/share/R/include/R_ext/Complex.h:80:6: warning: ISO C99 doesn’t support unnamed structs/unions [-Wpedantic]
#>    80 |     };
#>       |      ^
#> ✔ done
#> → finding duplicate expressions in FD model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in FD model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling EBE model...
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> In file included from /usr/share/R/include/R.h:71,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2parse.h:33,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2/include/rxode2.h:9,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2_model_shared.h:3,
#>                  from rx_f49a7e12ecd7da93c75a499a00195dc9_.c:117:
#> /usr/share/R/include/R_ext/Complex.h:80:6: warning: ISO C99 doesn’t support unnamed structs/unions [-Wpedantic]
#>    80 |     };
#>       |      ^
#> ✔ done
#> → compiling events FD model...
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> In file included from /usr/share/R/include/R.h:71,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2parse.h:33,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2/include/rxode2.h:9,
#>                  from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2_model_shared.h:3,
#>                  from rx_9ea8cbcfb816fb27bf1bcb9bc6560977_.c:117:
#> /usr/share/R/include/R_ext/Complex.h:80:6: warning: ISO C99 doesn’t support unnamed structs/unions [-Wpedantic]
#>    80 |     };
#>       |      ^
#> ✔ done
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 3688
#> → compress parHistData in nlmixr2 object, save 5576
unique(fitB$ka[fitB$ID%in%11:13])     # KA is estimated also for IV group
#> [1] 0.6451421 0.6451447 0.6451447
unique(fitB$eta.ka[fitB$ID%in%11:13]) # eta KA for IV is not 0 (would expect 0 or NA?)
#> [1] -3.898841e-06  9.097824e-11  1.501928e-08
fitB$parFixed # in here the backtransformed parameters are as expected (e.g. ka is backtransformed)
#>         Parameter   Est.     SE  %RSE       Back-transformed(95%CI) BSV(CV%)
#> tka        Log Ka -0.438  0.162    37           0.645 (0.47, 0.886)     26.0
#> tcl        Log Cl  -6.83  0.392  5.74     0.00108 (0.0005, 0.00233)     88.2
#> tv          Log V  -7.54 0.0471 0.625 0.000533 (0.000486, 0.000584)     3.16
#> prop.sd            0.273                                      0.273         
#>         Shrink(SD)%
#> tka          27.2%=
#> tcl         -7.96%>
#> tv           88.5%>
#> prop.sd

^{Created on 2024-02-02 with reprex v2.1.0}

mattfidler commented 9 months ago

This means:

If the eta is non-informative for saem, you can identify it manually with a solution like @RichardHooijmaijers with ROUTE==1.
The second solution is use focei

RichardHooijmaijers commented 8 months ago

@mattfidler thanks for looking into this! Will continue to use the initial solution. On a side note; it seems when implementing this solution, the ka is not back-transformed in the parameter table. Do you think there is a way to achieve this?

mattfidler commented 8 months ago

Yes, you can use backTransform()

In this case you would use:

    tka <- 0.45; backTransform("exp") # Log Ka

The reprex shows the back-transformed estimate:

library(nlmixr2est)
#> Loading required package: nlmixr2data
# Simulated dataset with 2 IDs IV and 2 IDs PO
dat <- structure(list(ID = c(11, 11, 11, 11, 11, 11, 11, 12, 12, 12, 
                                 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 13, 13, 21, 21, 21, 21, 
                                 21, 21, 21, 22, 22, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 23, 23), 
                         TIME = c(0, 0.05, 0.25, 0.5, 1, 3, 5, 0, 0.05, 0.25, 0.5, 
                                     1, 3, 5, 0, 0.05, 0.25, 0.5, 1, 3, 5, 8, 0, 0.25, 0.5, 1, 3, 
                                     5, 8, 0, 0.25, 0.5, 1, 3, 5, 8, 0, 0.25, 0.5, 1, 3, 5, 8), 
                         DV = c(NA,2017.85, 1323.74, 792.5, 822.72, 36.27, 3.33, NA, 1702, 1290.75, 
                                            1095.95, 907.6, 125.44, 14.44, NA, 1933.04, 1242.43, 661.22, 
                                            193.52, 1.75, NA, NA, NA, 706.58, 1063.14, 2257.62, 941.33, 629.69, 
                                            100, NA, 1462.95, 2217.76, 2739.5, 705.3, 108.47, 8.75, NA, 211.66, 
                                            467.23, 174.24, 153.6, 27.07, 2.81), 
                         AMT = c(1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 
                                    0, 0, 5, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0), 
                         EVID = c(1,0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 
                                       1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0), 
                         CMT = c(2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 1, 2, 
                                    2, 2, 2, 2, 2), 
                         DOSE = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5), 
                         ROUTE = c(1, 1, 1, 1, 1, 1, 1,1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)), 
                     row.names = c(NA, -43L), class = "data.frame")

# First fit with if statement to indicate diff IV/PO
modA <- function() {
ini({
    tka <- 0.45; backTransform("exp") # Log Ka
    tcl <- -7     # Log Cl
    tv  <- -8     # Log V
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    prop.sd <- 0.7
    })
model({
      if(ROUTE!=1) ka <- exp(tka + eta.ka)
      cl <- exp(tcl + eta.cl)
      v <- exp(tv + eta.v)
      d/dt(depot) = -ka * depot
      d/dt(center) = ka * depot - cl / v * center
      cp = center / v
      cp ~ prop(prop.sd)
      })
}
fitA <- nlmixr(modA, dat, est="focei", control=list(print=0))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: eta.ka
#> as a work-around try putting the mu-referenced expression on a simple line

#> Warning: some etas defaulted to non-mu referenced, possible parsing error: eta.ka
#> as a work-around try putting the mu-referenced expression on a simple line
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of full model...
#> ✔ done
#> → calculate jacobian
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate sensitivities
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate ∂(f)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate ∂(R²)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling inner model...
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> ✔ done
#> → finding duplicate expressions in FD model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in FD model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling EBE model...
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> ✔ done
#> → compiling events FD model...
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> ✔ done
#> rxode2 2.1.2 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 3688
#> → compress parHistData in nlmixr2 object, save 5256

print(fitA)
#> ── nlmixr² FOCEi (outer: nlminb) ──
#> 
#>           OBJF      AIC      BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 385.9319 464.2576 475.1451      -225.1288         7.41983        4.163486
#> 
#> ── Time (sec $time): ──
#> 
#>            setup optimize covariance table compress    other
#> elapsed 0.011995 0.224088   0.224089 0.027    0.007 3.350828
#> 
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#> 
#>         Parameter   Est.    SE %RSE      Back-transformed(95%CI) BSV(CV% or SD)
#> tka               -0.439 0.193 43.9 0.644721 (0.441769, 0.94091)               
#> tcl        Log Cl  -6.79 0.293 4.32    0.00113 (0.000635, 0.002)           89.9
#> tv          Log V  -7.44 0.153 2.06 0.00059 (0.000437, 0.000797)           30.1
#> prop.sd            0.328                                   0.328               
#> eta.ka                                                                    0.434
#>         Shrink(SD)%
#> tka                
#> tcl         -4.89% 
#> tv           60.9% 
#> prop.sd            
#> eta.ka       60.7% 
#>  
#>   Covariance Type ($covMethod): r,s
#>   No correlations in between subject variability (BSV) matrix
#>   Full BSV covariance ($omega) or correlation ($omegaR; diagonals=SDs) 
#>   Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink 
#>   Information about run found ($runInfo):
#>    • gradient problems with initial estimate and covariance; see $scaleInfo 
#>    • last objective function was not at minimum, possible problems in optimization 
#>    • ETAs were reset to zero during optimization; (Can control by foceiControl(resetEtaP=.)) 
#>    • initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=)) 
#>   Censoring ($censInformation): No censoring
#>   Minimization message ($message):  
#>     false convergence (8) 
#>   In an ODE system, false convergence may mean "useless" evaluations were performed.
#>   See https://tinyurl.com/yyrrwkce
#>   It could also mean the convergence is poor, check results before accepting fit
#>   You may also try a good derivative free optimization:
#>     nlmixr2(...,control=list(outerOpt="bobyqa"))
#> 
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 37 × 24
#>   ID     EVID  TIME    DV  PRED   RES  WRES IPRED   IRES  IWRES CPRED  CRES
#>   <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl>
#> 1 11        0  0.05 2018. 1540.  478. 0.723 1780.  238.   0.407 1532.  486.
#> 2 11        0  0.25 1324. 1051.  273. 0.503 1376.  -52.1 -0.115 1044.  279.
#> 3 11        0  0.5   792.  651.  141. 0.269  997. -204.  -0.625  630.  162.
#> # ℹ 34 more rows
#> # ℹ 12 more variables: CWRES <dbl>, eta.ka <dbl>, eta.cl <dbl>, eta.v <dbl>,
#> #   depot <dbl>, center <dbl>, ka <dbl>, cl <dbl>, v <dbl>, tad <dbl>,
#> #   dosenum <dbl>, ROUTE <dbl>

^{Created on 2024-02-08 with reprex v2.1.0}

RichardHooijmaijers commented 8 months ago

@mattfidler Thanks again, the backTransform() works as you said. Maybe one last thing; when I try to simulate with the model that include the if statement (if(ROUTE!=1) ka <- exp(tka + eta.ka)) I run into problems. When setting the ROUTE to a value of 2 an oral profile is simulated. However when using a value of 1 it results in all NA values (I used the model as stated above):

et <- et(amount.units='mg', time.units='hours') %>%
  et(dose=10) %>% et(seq(0,8,0.2)) %>% mutate(ROUTE=1)
sims <- rxSolve(modA,et)
plot(sims)

mattfidler commented 8 months ago

@mattfidler Thanks again, the backTransform() works as you said. Maybe one last thing; when I try to simulate with the model that include the if statement (if(ROUTE!=1) ka <- exp(tka + eta.ka)) I run into problems. When setting the ROUTE to a value of 2 an oral profile is simulated. However when using a value of 1 it results in all NA values (I used the model as stated above):
et <- et(amount.units='mg', time.units='hours') %>%
  et(dose=10) %>% et(seq(0,8,0.2)) %>% mutate(ROUTE=1)
sims <- rxSolve(modA,et)
plot(sims)

@RichardHooijmaijers

This is actually expected behavior: See:

https://github.com/nlmixr2/rxode2/discussions/655

Also based on that behavior, I would think that you should probably update the model a bit further to:

modA <- function() {
ini({
    tka <- 0.45; backTransform("exp") # Log Ka
    tcl <- -7     # Log Cl
    tv  <- -8     # Log V
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    prop.sd <- 0.7
    })
model({
     ka <- 0
      if(ROUTE!=1) ka <- exp(tka + eta.ka)
      cl <- exp(tcl + eta.cl)
      v <- exp(tv + eta.v)
      d/dt(depot) = -ka * depot
      d/dt(center) = ka * depot - cl / v * center
      cp = center / v
      cp ~ prop(prop.sd)
      })
}

To work around the issue.

RichardHooijmaijers commented 8 months ago

@mattfidler Thanks that makes sense! I will update the model as indicated

mattfidler commented 5 months ago

@RichardHooijmaijers

This is now taken care of in #418, though it is only applied in saem.

You can handle it in the old way or the new way that accounts for non-informative etas:

  library(nlmixr2)
#> Loading required package: nlmixr2data

  data <- data.frame(ID = c(11, 11, 11, 11, 11, 11, 11, 12, 12, 12,
                           12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 13, 13, 21, 21, 21, 21,
                           21, 21, 21, 22, 22, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 23, 23),
                    TIME = c(0, 0.05, 0.25, 0.5, 1, 3, 5, 0, 0.05, 0.25, 0.5,
                             1, 3, 5, 0, 0.05, 0.25, 0.5, 1, 3, 5, 8, 0, 0.25, 0.5, 1, 3,
                             5, 8, 0, 0.25, 0.5, 1, 3, 5, 8, 0, 0.25, 0.5, 1, 3, 5, 8),
                    DV = c(NA,2017.85, 1323.74, 792.5, 822.72, 36.27, 3.33, NA, 1702, 1290.75,
                           1095.95, 907.6, 125.44, 14.44, NA, 1933.04, 1242.43, 661.22,
                           193.52, 1.75, NA, NA, NA, 706.58, 1063.14, 2257.62, 941.33, 629.69,
                           100, NA, 1462.95, 2217.76, 2739.5, 705.3, 108.47, 8.75, NA, 211.66,
                           467.23, 174.24, 153.6, 27.07, 2.81),
                    AMT = c(1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0,
                            0, 0, 5, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0),
                    EVID = c(1,0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
                             1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0),
                    CMT = c(2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 1, 2,
                            2, 2, 2, 2, 2),
                    DOSE = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5),
                    ROUTE = c(1, 1, 1, 1, 1, 1, 1,1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2))

  # First fit with if statement to indicate diff IV/PO
  modA <- function() {
    ini({
      tka <- 0.45   # Log Ka
      tcl <- -7     # Log Cl
      tv  <- -8     # Log V
      eta.ka ~ 0.6
      eta.cl ~ 0.3
      eta.v ~ 0.1
      prop.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(center) <- ka * depot - cl/v * center
      cp <- center/v
      cp ~ prop(prop.sd)
    })
  }

  f <- nlmixr2(modA, data, "saem", control=saemControl(print=0))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of saem model...
#> ✔ done
#> → finding duplicate expressions in saem model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> ✔ done
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> ℹ calculate uninformed etas
#> ℹ done
#> rxode2 2.1.2.9000 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> Calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of saem model...
#> ✔ done
#> → finding duplicate expressions in saem predOnly model 0...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in saem predOnly model 1...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem predOnly model 1...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in saem predOnly model 2...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem predOnly model 2...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> ✔ done
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 3696
#> → compress phiM in nlmixr2 object, save 48976
#> → compress parHistData in nlmixr2 object, save 14768
#> → compress saem0 in nlmixr2 object, save 24544


  print(f$etaObf)
#>   ID        eta.ka      eta.cl        eta.v OBJI
#> 1 11  0.0000000000 -0.41570114  0.006239803   NA
#> 2 12  0.0000000000 -0.68575973  0.029676164   NA
#> 3 13  0.0000000000  0.16399354 -0.018635434   NA
#> 4 21  0.0003296939 -0.61660470  0.083404192   NA
#> 5 22  0.0006148233  0.02105591 -0.079070131   NA
#> 6 23 -0.0009119386  1.54935134  0.010834432   NA

  print(f$parFixed)
#>         Parameter   Est.     SE %RSE      Back-transformed(95%CI) BSV(CV%)
#> tka        Log Ka -0.392 0.0584 14.9         0.676 (0.602, 0.758)    0.444
#> tcl        Log Cl  -6.88  0.318 4.63  0.00103 (0.000553, 0.00193)     89.3
#> tv          Log V  -7.53  0.107 1.42 0.000535 (0.000434, 0.00066)     10.3
#> prop.sd             0.33                                     0.33         
#>         Shrink(SD)%
#> tka          88.4%>
#> tcl         -8.46%>
#> tv           47.6%>
#> prop.sd


  f2 <- nlmixr2(modA, data, "saem", control=saemControl(handleUninformativeEtas=FALSE, print=0))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of saem model...
#> ✔ done
#> → finding duplicate expressions in saem model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> ✔ done
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> Calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of saem model...
#> ✔ done
#> → finding duplicate expressions in saem predOnly model 0...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in saem predOnly model 1...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem predOnly model 1...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in saem predOnly model 2...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem predOnly model 2...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> ✔ done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 3696
#> → compress phiM in nlmixr2 object, save 42976
#> → compress parHistData in nlmixr2 object, save 14376
#> → compress saem0 in nlmixr2 object, save 24328


  print(f2$etaObf)
#>   ID       eta.ka     eta.cl         eta.v OBJI
#> 1 11  0.026674327 -0.4397220 -0.0030638187   NA
#> 2 12 -0.011641037 -0.7248965  0.0039533600   NA
#> 3 13  0.023496425  0.1568911 -0.0010848195   NA
#> 4 21 -0.277940988 -0.4199011  0.0014437556   NA
#> 5 22  0.323803630 -0.1291742 -0.0067886909   NA
#> 6 23  0.004207461  1.5435251 -0.0005657952   NA

  print(f2$parFixed)
#>         Parameter   Est.    SE %RSE      Back-transformed(95%CI) BSV(CV%)
#> tka        Log Ka -0.462 0.169 36.5          0.63 (0.453, 0.877)     28.2
#> tcl        Log Cl  -6.85  0.31 4.53  0.00106 (0.000575, 0.00194)     87.0
#> tv          Log V  -7.52 0.083  1.1 0.000541 (0.00046, 0.000636)     2.59
#> prop.sd            0.271                                   0.271         
#>         Shrink(SD)%
#> tka          31.0%>
#> tcl         -8.49%>
#> tv           85.7%>
#> prop.sd

^{Created on 2024-05-17 with reprex v2.1.0}

mattfidler commented 5 months ago

This has been integrated into the development version of nlmixr2est. Currently the un-informative etas are only detected for saem which is consistent with the observed behavior.

nlmixr2 / nlmixr2est

`saem` non-informative `ETAs` #418