nlmixr2est focei needs to cache based on theta/omega order.
#415
Comments
|
With the fix I get: > setwd('/home/matt/src/issues')
> source("nlmixr2-194.R", echo=TRUE)
> library(nlmixr2)
Loading required package: nlmixr2data
> mod <- function(){
+ > model({
+ > d/dt(depot) = -ka * depot
+ > d/dt(center) = ka * depot - cl / v * center
+ > })}
> ### piping option 1
> theta_t <- c("tka"=log(1), "tcl" = log(0.1), "v" = 1 )
> modSim <- mod %>%
+ > model(ka <- exp(tka + bsv.ka), append=NA) |>
+ > ini(tka=log(1), bsv.ka ~ 0.1) |>
+ > model(cl <- exp(tcl + bsv.cl), append=NA .... [TRUNCATED]
ℹ promote `tka` to population parameter with initial estimate 1
ℹ promote `bsv.ka` to between subject variability with initial estimate 1
ℹ change initial estimate of `tka` to `0`
ℹ change initial estimate of `bsv.ka` to `0.1`
ℹ promote `tcl` to population parameter with initial estimate 1
ℹ promote `bsv.cl` to between subject variability with initial estimate 1
ℹ change initial estimate of `tcl` to `-0.693147180559945`
ℹ change initial estimate of `bsv.cl` to `0.1`
ℹ change initial estimate of `tka` to `0`
ℹ change initial estimate of `tcl` to `-2.30258509299405`
ℹ promote `v` to population parameter with initial estimate 1
ℹ change initial estimate of `v` to `1`
ℹ add residual parameter `prop.err` and set estimate to 1
ℹ change initial estimate of `prop.err` to `0.1`
> theta <- c("ka"=1, "cl" = 1, "v" = 1)
> inits <- c(depot=0, centr=0 )
> et <- eventTable( ) %>%
+ > add.dosing( dose=1, start.time=0 ) %>%
+ > add.sampling( time = seq(0, 24, 1)) |>
+ > et(id=1:10)
> sim1 <- rxSolve(modSim, events = et, addDosing=TRUE) |>
+ > dplyr::rename(DV=sim)
using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
> mod1 <- mod %>%
+ > model(ka <- exp(tka + bsv.ka), append=NA) |>
+ > ini(tka=log(1), bsv.ka ~ 0.1) |>
+ > model(cl <- exp(tcl + bsv.cl), append=NA) .... [TRUNCATED]
ℹ promote `tka` to population parameter with initial estimate 1
ℹ promote `bsv.ka` to between subject variability with initial estimate 1
ℹ change initial estimate of `tka` to `0`
ℹ change initial estimate of `bsv.ka` to `0.1`
ℹ promote `tcl` to population parameter with initial estimate 1
ℹ promote `bsv.cl` to between subject variability with initial estimate 1
ℹ change initial estimate of `tcl` to `-0.693147180559945`
ℹ change initial estimate of `bsv.cl` to `0.1`
ℹ change initial estimate of `tka` to `0`
ℹ change initial estimate of `tcl` to `-2.30258509299405`
ℹ promote `v` to population parameter with initial estimate 1
ℹ change initial estimate of `v` to `1`
ℹ add residual parameter `prop.err` and set estimate to 1
ℹ change initial estimate of `prop.err` to `0.1`
> ### piping option 2
>
> mod2<- mod |>
+ > model(ka <- exp(tka + bsv.ka), append=NA) |>
+ > ini(tka=log(1), bsv.ka ~ 0.1) |>
+ > model(cl <- exp(tcl .... [TRUNCATED]
ℹ promote `tka` to population parameter with initial estimate 1
ℹ promote `bsv.ka` to between subject variability with initial estimate 1
ℹ change initial estimate of `tka` to `0`
ℹ change initial estimate of `bsv.ka` to `0.1`
ℹ promote `tcl` to population parameter with initial estimate 1
ℹ promote `bsv.cl` to between subject variability with initial estimate 1
ℹ change initial estimate of `tcl` to `-0.693147180559945`
ℹ change initial estimate of `bsv.cl` to `0.1`
ℹ add residual parameter `prop.err` and set estimate to 1
ℹ change initial estimate of `prop.err` to `0.1`
ℹ change initial estimate of `tka` to `0`
ℹ change initial estimate of `tcl` to `-2.30258509299405`
ℹ promote `v` to population parameter with initial estimate 1
ℹ change initial estimate of `v` to `1`
> ### fit the data with two models
> fit1 <- nlmixr(mod1, sim1, est = "focei",
+ > control = foceiControl(print=0),
+ > t .... [TRUNCATED]
→ loading into symengine environment...
→ pruning branches (`if`/`else`) of full model...
✔ done
→ calculate jacobian
→ calculate sensitivities
→ calculate ∂(f)/∂(η)
→ calculate ∂(R²)/∂(η)
→ finding duplicate expressions in inner model...
→ optimizing duplicate expressions in inner model...
→ finding duplicate expressions in EBE model...
→ optimizing duplicate expressions in EBE model...
→ compiling inner model...
using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
✔ done
→ finding duplicate expressions in FD model...
→ optimizing duplicate expressions in FD model...
→ 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.1 using 8 threads (see ?getRxThreads)
no cache: create with `rxCreateCache()`
calculating covariance matrix
[====|====|====|====|====|====|====|====|====|====] 100%; 0:00:00 done
→ Calculating residuals/tables
✔ done
→ compress origData in nlmixr2 object, save 15640
→ compress parHistData in nlmixr2 object, save 4784
> fit2 <- nlmixr(mod2, sim1, est = "focei",
+ > control = foceiControl(print=0),
+ > table=list(cwres=TRUE))
→ loading into symengine environment...
→ pruning branches (`if`/`else`) of full model...
✔ done
→ calculate jacobian
→ calculate sensitivities
→ calculate ∂(f)/∂(η)
→ calculate ∂(R²)/∂(η)
→ finding duplicate expressions in inner model...
→ optimizing duplicate expressions in inner model...
→ finding duplicate expressions in EBE model...
→ optimizing duplicate expressions in EBE model...
→ compiling inner model...
using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
✔ done
→ finding duplicate expressions in FD model...
→ optimizing duplicate expressions in FD model...
→ 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
calculating covariance matrix
[====|====|====|====|====|====|====|====|====|====] 100%; 0:00:00 done
→ Calculating residuals/tables
✔ done
→ compress origData in nlmixr2 object, save 15640
→ compress parHistData in nlmixr2 object, save 4792
> summary(fit1$foceiModel$inner)
rxode2 2.1.1 model named rx_2747c9a2b83e99193adcb9fef6697214 model (⚠ unloaded reload with rxode2::rxLoad).
DLL: /tmp/RtmpboEZRo/rxode2/rx_2747c9a2b83e99193adcb9fef6697214__.rxd/rx_2747c9a2b83e99193adcb9fef6697214_.so
NULL
Calculated Variables:
[1] "rx_pred_" "rx__sens_rx_pred__BY_ETA_1___" "rx__sens_rx_pred__BY_ETA_2___"
[4] "rx_r_" "rx__sens_rx_r__BY_ETA_1___" "rx__sens_rx_r__BY_ETA_2___"
── rxode2 Model Syntax ──
rxode2({
param(THETA[1], THETA[2], THETA[3], THETA[4], ETA[1], ETA[2])
cmt(depot)
cmt(center)
rx_expr_0 ~ ETA[1] + THETA[1]
rx_expr_4 ~ exp(rx_expr_0)
d/dt(depot) = -rx_expr_4 * depot
rx_expr_1 ~ ETA[2] + THETA[2]
rx_expr_5 ~ exp(rx_expr_1)
rx_expr_6 ~ rx_expr_4 * depot
rx_expr_7 ~ rx_expr_5 * center
rx_expr_8 ~ rx_expr_7/THETA[3]
d/dt(center) = rx_expr_6 - rx_expr_8
rx_expr_9 ~ rx_expr_4 * rx__sens_depot_BY_ETA_1___
d/dt(rx__sens_depot_BY_ETA_1___) = -rx_expr_4 * depot - rx_expr_9
d/dt(rx__sens_center_BY_ETA_1___) = rx_expr_6 + rx_expr_9 -
rx_expr_5 * rx__sens_center_BY_ETA_1___/THETA[3]
d/dt(rx__sens_depot_BY_ETA_2___) = -rx_expr_4 * rx__sens_depot_BY_ETA_2___
d/dt(rx__sens_center_BY_ETA_2___) = rx_expr_4 * rx__sens_depot_BY_ETA_2___ -
rx_expr_8 - rx_expr_5 * rx__sens_center_BY_ETA_2___/THETA[3]
rx_yj_ ~ 2
rx_lambda_ ~ 1
rx_hi_ ~ 1
rx_low_ ~ 0
rx_pred_ = center
rx__sens_rx_pred__BY_ETA_1___ = rx__sens_center_BY_ETA_1___
rx__sens_rx_pred__BY_ETA_2___ = rx__sens_center_BY_ETA_2___
rx_expr_2 ~ center * THETA[4]
rx_r_ = Rx_pow_di((rx_expr_2), 2)
rx_expr_3 ~ 2 * (rx_expr_2)
rx__sens_rx_r__BY_ETA_1___ = rx_expr_3 * rx__sens_center_BY_ETA_1___ *
THETA[4]
rx__sens_rx_r__BY_ETA_2___ = rx_expr_3 * rx__sens_center_BY_ETA_2___ *
THETA[4]
dvid(2)
})
> AIC(fit1)
[1] -1178.98
> summary(fit2$foceiModel$inner)
rxode2 2.1.1 model named rx_cf0601084d9083bbbb84f8fe6679e0af model (✔ ready).
DLL: /tmp/RtmpboEZRo/rxode2/rx_cf0601084d9083bbbb84f8fe6679e0af__.rxd/rx_cf0601084d9083bbbb84f8fe6679e0af_.so
NULL
Calculated Variables:
[1] "rx_pred_" "rx__sens_rx_pred__BY_ETA_1___" "rx__sens_rx_pred__BY_ETA_2___"
[4] "rx_r_" "rx__sens_rx_r__BY_ETA_1___" "rx__sens_rx_r__BY_ETA_2___"
── rxode2 Model Syntax ──
rxode2({
param(THETA[1], THETA[2], THETA[3], THETA[4], ETA[1], ETA[2])
cmt(depot)
cmt(center)
rx_expr_0 ~ ETA[1] + THETA[1]
rx_expr_4 ~ exp(rx_expr_0)
d/dt(depot) = -rx_expr_4 * depot
rx_expr_1 ~ ETA[2] + THETA[2]
rx_expr_5 ~ exp(rx_expr_1)
rx_expr_6 ~ rx_expr_4 * depot
rx_expr_7 ~ rx_expr_5 * center
rx_expr_8 ~ rx_expr_7/THETA[4]
d/dt(center) = rx_expr_6 - rx_expr_8
rx_expr_9 ~ rx_expr_4 * rx__sens_depot_BY_ETA_1___
d/dt(rx__sens_depot_BY_ETA_1___) = -rx_expr_4 * depot - rx_expr_9
d/dt(rx__sens_center_BY_ETA_1___) = rx_expr_6 + rx_expr_9 -
rx_expr_5 * rx__sens_center_BY_ETA_1___/THETA[4]
d/dt(rx__sens_depot_BY_ETA_2___) = -rx_expr_4 * rx__sens_depot_BY_ETA_2___
d/dt(rx__sens_center_BY_ETA_2___) = rx_expr_4 * rx__sens_depot_BY_ETA_2___ -
rx_expr_8 - rx_expr_5 * rx__sens_center_BY_ETA_2___/THETA[4]
rx_yj_ ~ 2
rx_lambda_ ~ 1
rx_hi_ ~ 1
rx_low_ ~ 0
rx_pred_ = center
rx__sens_rx_pred__BY_ETA_1___ = rx__sens_center_BY_ETA_1___
rx__sens_rx_pred__BY_ETA_2___ = rx__sens_center_BY_ETA_2___
rx_expr_2 ~ center * THETA[3]
rx_r_ = Rx_pow_di((rx_expr_2), 2)
rx_expr_3 ~ 2 * (rx_expr_2)
rx__sens_rx_r__BY_ETA_1___ = rx_expr_3 * rx__sens_center_BY_ETA_1___ *
THETA[3]
rx__sens_rx_r__BY_ETA_2___ = rx_expr_3 * rx__sens_center_BY_ETA_2___ *
THETA[3]
dvid(2)
})
> AIC(fit2)
[1] -1175.341
> identical(mod1$theta, mod2$theta)
[1] FALSE
> |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
I couldn't reproduce your example so I had to modify it; The problem is rxode2 focei code between the 2 piping options are identical (which you can see at the end of the code), which means the caching doesn't consider the theta order (and should for
focei). I will file a bug with this.Here is the bug reprex:
However, if you add a
rxode2::rxClean()to remove the cache before the second run the two fits they are actually similar:Note that they are not the same. Order of parameters does matter in the final fit, though it shouldn't be as bad as you observed
Originally posted by @mattfidler in nlmixr2/nlmixr2#194 (reply in thread)
The text was updated successfully, but these errors were encountered: