mattfidler
commented
5 months ago Here is the reprex; the underlying model without etas solves correctly:
library(nlmixr2est)
#> Loading required package: nlmixr2data
#> detected new version of nlmixr2est, cleaning rxode2 cache
## The basic model consists of an ini block that has initial estimates
one.compartment <- function() {
ini({
tka <- log(1.57); label("Ka")
tcl <- log(2.72); label("Cl")
tv <- log(31.5); label("V")
eta.ka ~ 0.6
eta.cl ~ 0.3
eta.v ~ 0.1
add.sd <- 0.7
})
# and a model block with the error specification and model specification
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 ~ add(add.sd)
})
}
dat<-theo_md
dat[dat$EVID==101,]$EVID<-1
dat$MDV<-dat$EVID
dat[6:8,]$MDV=1 # example, set line 6-8 observation MDV as 1
## The fit is performed by the function nlmixr/nlmixr2 specifying the model, data and estimate
fit0 <- nlmixr2(one.compartment,dat, est="saem", saemControl(print=0))
#> → 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’
#> 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
#> ✔ 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 17312
#> → compress phiM in nlmixr2 object, save 63528
#> → compress parHistData in nlmixr2 object, save 14304
#> → compress saem0 in nlmixr2 object, save 29160
fit0[1:10,c("ID", "EVID", "TIME", "tad")] |> as.data.frame()
#> ID EVID TIME tad
#> 1 1 0 0.00 0.00
#> 2 1 0 0.25 0.25
#> 3 1 0 0.57 0.57
#> 4 1 0 1.12 1.12
#> 5 1 2 2.02 0.00
#> 6 1 2 3.82 0.00
#> 7 1 2 5.10 0.00
#> 8 1 0 7.03 1.93
#> 9 1 0 9.05 3.95
#> 10 1 0 12.12 7.02
# Using the underlying rxode2 model with no eta gives the right results:
library(rxode2)
mod <- rxode2({
param(THETA[1], THETA[2], THETA[3], THETA[4], ETA[1], ETA[2],
ETA[3])
cmt(depot)
cmt(center)
rx_pred_ = NA
rx_r_ = NA
tka = THETA[1] + ETA[1]
tcl = THETA[2] + ETA[2]
tv = THETA[3] + ETA[3]
add.sd = THETA[4]
rx_yj_ ~ 2
rx_lambda_ ~ 1
rx_hi_ ~ 1
rx_low_ ~ 0
rx_expr_0 ~ exp(tka)
d/dt(depot) = -rx_expr_0 * depot
d/dt(center) = rx_expr_0 * depot - exp(tcl - tv) * center
rx_pred_ = exp(-tv) * center
rx_r_ = Rx_pow_di(add.sd, 2)
ka = exp(tka)
cl = exp(tcl)
v = exp(tv)
cp = exp(-tv) * center
tad = tad()
dosenum = dosenum()
tka = THETA[1]
tcl = THETA[2]
tv = THETA[3]
add.sd = THETA[4]
eta.ka = ETA[1]
eta.cl = ETA[2]
eta.v = ETA[3]
cmt(cp)
dvid(3)
})
par <- c("THETA[1]" = 0.271646397929738, "THETA[2]" = 1.04837208911781, "THETA[3]" = 3.41786691697643, "THETA[4]" = 1.02739485488012, "ETA[1]" = 0, "ETA[2]" = 0, "ETA[3]" = 0)
dat <- nlmixr2data::theo_md
dat[dat$EVID==101,]$EVID<-1
dat$MDV<-dat$EVID
dat[6:8,]$MDV=1 # example, set line 6-8 observation MDV as 1
df <- rxSolve(mod, par, dat)
df[1:10,c("id", "evid", "time", "tad")]
#> id evid time tad
#> 1 1 0 0.00 0.00
#> 2 1 0 0.25 0.25
#> 3 1 0 0.57 0.57
#> 4 1 0 1.12 1.12
#> 5 1 2 2.02 2.02
#> 6 1 2 3.82 3.82
#> 7 1 2 5.10 5.10
#> 8 1 0 7.03 7.03
#> 9 1 0 9.05 9.05
#> 10 1 0 12.12 12.12Created on 2024-02-05 with reprex v2.1.0
I get a question on the tad (time after the dose) calculation provided by nlmixr2 output. I will be appreciated if you have any idea on it
I find that when I try to set part of observations' MDV to 1 (keep EVID=0) in the PopPK dataset, the tad (time after the last dose) calculated by nlmixr2 become confused. (This wont happen when I use LinCMT(), only happen when it is ODE)
I use theo_md as example, most of the following code are same with the nlmixr2 website's. The only difference is I try to set row 6-8 observation MDV as 1
I select part of columns from nlmxir2 output after running. and it can be seen that the tad is totally changed, and it seems that those observation lines (MDV=1,EVID=0) were considered as a dosing line/event line as well. After line 6-8, all of the tad and dosenum were wrong actually.
The interesting thing is that when I use vpcSim(fit,n = no.sim) to see the simulation results, the tad become normal again.
Have you met this issue before?