`nlmixr2est` focei needs to cache based on theta/omega order.

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:


> source("nlmixr2-194.R", echo=TRUE)

> library(nlmixr2)

> 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)

> 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] 
calculating covariance matrix
[====|====|====|====|====|====|====|====|====|====] 100%; 0:00:00 done
→ Calculating residuals/tables
✔ done
→ compress origData in nlmixr2 object, save 15592
→ compress parHistData in nlmixr2 object, save 4816

> fit2 <- nlmixr(mod2, sim1, est = "focei",
+                control = foceiControl(print=0),
+                table=list(cwres=TRUE))
calculating covariance matrix
[====|====|====|====|====|====|====|====|====|====] 100%; 0:00:00 done
→ Calculating residuals/tables
✔ done
→ compress origData in nlmixr2 object, save 15592
→ compress parHistData in nlmixr2 object, save 4040

> summary(fit1$foceiModel$inner)
rxode2 2.1.1 model named rx_2747c9a2b83e99193adcb9fef6697214 model (✔ ready). 
DLL: /tmp/RtmpWWeyW0/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] -1378.16

> summary(fit2$foceiModel$inner)
rxode2 2.1.1 model named rx_2747c9a2b83e99193adcb9fef6697214 model (✔ ready). 
DLL: /tmp/RtmpWWeyW0/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(fit2)
[1] 8711730
>

However, if you add a rxode2::rxClean() to remove the cache before the second run the two fits they are actually similar:

> 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 15624
→ compress parHistData in nlmixr2 object, save 5160

> rxode2::rxClean()

> 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 15624
→ compress parHistData in nlmixr2 object, save 4856

> summary(fit1$foceiModel$inner)
rxode2 2.1.1 model named rx_2747c9a2b83e99193adcb9fef6697214 model (✖ invalid re-create with rxode2::rxode2). 
DLL: /tmp/RtmpHXWKON/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] -1302.015

> summary(fit2$foceiModel$inner)
rxode2 2.1.1 model named rx_cf0601084d9083bbbb84f8fe6679e0af model (✔ ready). 
DLL: /tmp/RtmpHXWKON/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] -1302.527

> identical(mod1$theta, mod2$theta)
[1] FALSE
>

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 https://github.com/nlmixr2/nlmixr2/discussions/194#discussioncomment-8073709

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
>

nlmixr2 / nlmixr2est

`nlmixr2est` focei needs to cache based on theta/omega order. #415