Skip to content

Commit f872c30

Browse files
mattfidlermattfidler
committed
fix(poped): keep the if() guard on adaptive dosing calls
The PopED model translation flattens if/else so that no stray lhs statements are left behind, but it called `.rxPrune()` without the capture environment. Adaptive dosing calls (`evid_()`, `bolus()`, `infuse()`, `infuseDur()`, `reset()`, ...) then fell through to the generic call branch and the condition guarding them was silently dropped, so `if (t <= 0) infuseDur(DOSE, TINF, cmt=depot)` pushed a dose at *every* design point instead of once (the internal mtime model runs straight through `maxExtra=1000`). Use the pruner's `..captureN`/`..capturedEvid` protocol and re-emit each captured call as `if (rxCaptureId#) { ... }` after the branches are flattened. Also add two PopED examples for the HCV model showing how to make the regimen itself optimizable: - dose-and-tinf: keeps the dose record, gets the amount and infusion duration out of `f(depot) <- DOSE` (amt=1) and `dur(depot) <- TINF` (rate=-2) - adaptive-dosing: no dose records at all, `infuseDur()` pushes the regimen, which makes the dosing interval a design variable too (needs rxode2 > 5.1.7 for the pushed-dose timing fix, rxode2#1214) Both agree on the reference design (OFV 88.27) and optimize with `poped_optim(..., opt_a=TRUE)`. Re #131
1 parent 031c1d6 commit f872c30

4 files changed

Lines changed: 310 additions & 1 deletion

File tree

‎NEWS.md‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,29 @@
11
# babelmixr2 0.1.11.9000
22

3+
* The PopED model translation no longer drops the `if ()` condition that
4+
guards an adaptive dosing call (`evid_()`, `bolus()`, `infuse()`,
5+
`infuseDur()`, `reset()`, ...). The branch pruner used to flatten the
6+
model unconditionally, so a model like `if (t <= 0) infuseDur(DOSE,
7+
TINF, cmt=1)` pushed a dose at *every* design point instead of once
8+
(#131). The pruner's capture protocol is now used and the guarded call
9+
is restored after the branches are flattened.
10+
11+
* Added two PopED examples showing how to make the dosing regimen itself
12+
optimizable (#131):
13+
14+
- `inst/poped/ex.10.PKPD.HCV.dose-and-tinf.babelmixr2.R` keeps the dose
15+
record and makes the amount and infusion duration design (`a`)
16+
variables via `f(depot) <- DOSE` (with `amt=1`) and
17+
`dur(depot) <- TINF` (with `rate=-2`).
18+
19+
- `inst/poped/ex.10.PKPD.HCV.adaptive-dosing.babelmixr2.R` drops the
20+
dose records entirely and pushes the regimen from inside the model
21+
with `infuseDur()`, which makes the dosing *interval* a design
22+
variable as well. This one needs rxode2 > 5.1.7 (rxode2#1214).
23+
24+
Both are optimized with `poped_optim(..., opt_a=TRUE)` and agree on the
25+
reference design (OFV 88.27).
26+
327
* The NONMEM/Monolix fit cache is now written with `saveRDS()` as
428
`<model>.rds` / `nlmixr.rds` instead of `qs2`, so `qs2` moved from
529
`Imports` to `Suggests`. Existing run directories keep working: a

‎R/poped.R‎

Lines changed: 24 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -875,7 +875,19 @@ rxUiGet.popedRxmodelBase <- function(x, ...) {
875875
.mod <- .mod[.w]
876876
.errDf <- .iniDf[!is.na(.iniDf$err), ,drop=FALSE]
877877
# remove if/else so extra lhs statements are not hanging around
878-
.mod <- str2lang(paste0("{", rxode2::.rxPrune(as.call(c(quote(`{`), .mod))), "}"))
878+
#
879+
# Adaptive dosing calls (`evid_()`, `bolus()`, `infuse()`, `infuseDur()`,
880+
# ...) cannot be flattened this way; the pruner captures them (and the
881+
# conditions they were nested under) when `..captureN`/`..capturedEvid`
882+
# are present in the pruning environment. Without this the enclosing
883+
# `if ()` is dropped and the dose fires at every event time.
884+
.pruneEnv <- new.env(parent=emptyenv())
885+
.pruneEnv$.if <- NULL
886+
.pruneEnv$.def1 <- NULL
887+
.pruneEnv$..captureN <- 0L
888+
.pruneEnv$..capturedEvid <- list()
889+
.mod <- str2lang(paste0("{", rxode2::.rxPrune(as.call(c(quote(`{`), .mod)),
890+
envir=.pruneEnv), "}"))
879891
.mod <- lapply(seq_along(.mod)[-1],
880892
function(i) { .mod[[i]]})
881893
.mod <- lapply(seq_along(.mod),
@@ -896,6 +908,17 @@ rxUiGet.popedRxmodelBase <- function(x, ...) {
896908
}
897909
.replaceErrWithConst(.cur, .errDf)
898910
})
911+
# restore the adaptive dosing calls captured by the pruner; the
912+
# `rxCaptureId#` variables hold the (flattened) condition they were
913+
# nested under.
914+
if (length(.pruneEnv$..capturedEvid) > 0L) {
915+
.mod <- c(.mod,
916+
lapply(.pruneEnv$..capturedEvid,
917+
function(cap) {
918+
str2lang(paste0("if (", cap$capVar, ") { ",
919+
cap$original, " }"))
920+
}))
921+
}
899922
.mod
900923
}
901924
attr(rxUiGet.popedRxmodelBase, "desc") <- "This gets the base rxode2 model for PopED"

‎inst/poped/ex.10.PKPD.HCV.adaptive-dosing.babelmixr2.R‎

Lines changed: 130 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,130 @@
1+
## HCV example (Nyberg et al., Br. J. Clin. Pharm., 2014) where the whole
2+
## regimen -- dose amount, infusion duration *and* dosing interval -- is
3+
## registered from inside the model with rxode2's adaptive dosing
4+
## functions, and all three are PopED design ("a") variables.
5+
## See https://github.com/nlmixr2/babelmixr2/issues/131
6+
##
7+
## Compare with ex.10.PKPD.HCV.dose-and-tinf.babelmixr2.R, which gets the
8+
## amount and the duration out of a dose record with `f()`/`dur()`. That
9+
## approach is simpler and needs nothing from the solver, but the dosing
10+
## *interval* stays baked into the event table. Here the design dataset
11+
## carries no dose records at all: `infuseDur()` pushes them, so `ii` is
12+
## just another covariate and can be optimized too.
13+
##
14+
## Requires rxode2 > 5.1.7 (the pushed-dose timing fix, rxode2#1214).
15+
16+
library(babelmixr2)
17+
library(PopED)
18+
19+
f <- function() {
20+
ini({
21+
tp <- fix(100)
22+
td <- fix(0.001)
23+
te <- fix(1e-7)
24+
ts <- fix(20000)
25+
26+
tKA <- log(0.8)
27+
tKE <- log(0.15)
28+
tVD <- log(100)
29+
tEC50 <- log(0.12)
30+
tn <- log(2)
31+
tdelta <- log(0.2)
32+
tc <- log(7)
33+
34+
eta.KA ~ 0.25
35+
eta.KE ~ 0.25
36+
eta.VD ~ 0.25
37+
eta.EC50 ~ 0.25
38+
eta.n ~ 0.25
39+
eta.delta ~ 0.25
40+
eta.c ~ 0.25
41+
42+
add.sd.pk <- sqrt(0.04) # nlmixr2 uses sd
43+
add.sd.pd <- sqrt(0.04)
44+
})
45+
model({
46+
p <- tp
47+
d <- td
48+
e <- te
49+
s <- ts
50+
KA <- exp(tKA + eta.KA)
51+
KE <- exp(tKE + eta.KE)
52+
VD <- exp(tVD + eta.VD)
53+
EC50 <- exp(tEC50 + eta.EC50)
54+
n <- exp(tn + eta.n)
55+
delta <- exp(tdelta + eta.delta)
56+
c <- exp(tc + eta.c)
57+
58+
## The design variables have to be copied into model variables before
59+
## they are handed to infuseDur(); a covariate that appears *only* as
60+
## an adaptive dosing argument is not registered as a parameter and
61+
## the model will not compile (rxode2#1231).
62+
amtI <- DOSE
63+
durI <- TINF
64+
iiI <- TAU
65+
66+
d/dt(depot) <- -KA*depot
67+
d/dt(central) <- KA*depot - KE*central
68+
d/dt(TC) <- s - TC*(e*VP + d) # target cells (TC)
69+
d/dt(IC) <- e*TC*VP - delta*IC # productively infected cells
70+
d/dt(VP) <- p*(1 - (pow(central/VD, n)/(pow(central/VD, n) +
71+
pow(EC50, n))))*IC - c*VP
72+
73+
TC(0) <- c*delta/(p*e)
74+
IC(0) <- (s*e*p - d*c*delta)/(p*delta*e)
75+
VP(0) <- (s*e*p - d*c*delta)/(c*delta*e)
76+
77+
## four infusions of DOSE over TINF, every TAU
78+
if (t <= 0) {
79+
infuseDur(amtI, durI, cmt=depot, ii=iiI, addl=3)
80+
}
81+
82+
conc <- central/VD
83+
eff <- log10(VP)
84+
conc ~ add(add.sd.pk)
85+
eff ~ add(add.sd.pd)
86+
})
87+
}
88+
89+
## observation-only design: there are no dose records to write down
90+
tms <- c(0, 0.25, 0.5, 1, 2, 3, 4, 7, 10, 14, 21, 28)
91+
e1 <- as.data.frame(et(tms))
92+
e1$dvid <- 1
93+
e2 <- e1
94+
e2$dvid <- 2
95+
e <- rbind(e1, e2)
96+
97+
babel.db <- nlmixr2(f, e, "poped",
98+
popedControl(groupsize=30,
99+
a=list(c(DOSE=180, TINF=1, TAU=7)),
100+
mina=c(DOSE=20, TINF=0.1, TAU=1),
101+
maxa=c(DOSE=600, TINF=12, TAU=14)))
102+
103+
plot_model_prediction(babel.db, facet_scales="free")
104+
105+
evaluate_design(babel.db)$ofv
106+
#> [1] 88.27007
107+
## ...which is the same design as the `f()`/`dur()` version of this
108+
## example (88.26974); the two parameterizations agree.
109+
110+
## the interval is now a design variable like any other
111+
for (ta in c(1, 2, 3, 5, 7, 10, 14)) {
112+
.db <- babel.db
113+
.db$design$a[1, "TAU"] <- ta
114+
message(sprintf("TAU = %5.1f ofv = %.4f", ta, evaluate_design(.db)$ofv))
115+
}
116+
#> TAU = 1.0 ofv = 91.0976
117+
#> TAU = 2.0 ofv = 90.6507
118+
#> TAU = 3.0 ofv = 90.2115
119+
#> TAU = 5.0 ofv = 89.4562
120+
#> TAU = 7.0 ofv = 88.2701
121+
#> TAU = 10.0 ofv = 89.1156
122+
#> TAU = 14.0 ofv = 88.8393
123+
124+
r <- poped_optim(babel.db, opt_a=TRUE, opt_xt=FALSE, parallel=FALSE)
125+
#> Optimized Covariates:
126+
#> Group 1: 1 : 600 : 7.1429 : 1
127+
#>
128+
#> OFV = 100.785
129+
#
130+
# (a random search, so the exact optimum moves between runs)

‎inst/poped/ex.10.PKPD.HCV.dose-and-tinf.babelmixr2.R‎

Lines changed: 132 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,132 @@
1+
## HCV example (Nyberg et al., Br. J. Clin. Pharm., 2014) with the dose
2+
## amount *and* the infusion duration treated as PopED design ("a")
3+
## variables. See https://github.com/nlmixr2/babelmixr2/issues/131
4+
##
5+
## The trick is to keep the dose *record* in the design dataset, but to
6+
## make both of its interesting quantities model quantities:
7+
##
8+
## - `amt=1` in the data + `f(depot) <- DOSE` in the model
9+
## => the administered amount is the covariate `DOSE`
10+
## - `rate=-2` in the data + `dur(depot) <- TINF` in the model
11+
## => the infusion duration is the covariate `TINF`
12+
##
13+
## Because `DOSE` and `TINF` are ordinary covariates as far as
14+
## rxode2/babelmixr2 are concerned, PopED sees them as elements of `a`
15+
## and can optimize over them with `opt_a=TRUE`.
16+
##
17+
## The dosing *interval* is still fixed by the event table here. See
18+
## ex.10.PKPD.HCV.adaptive-dosing.babelmixr2.R for the variant that
19+
## pushes the whole regimen from inside the model with `infuseDur()`,
20+
## which makes `ii` a design variable too.
21+
22+
library(babelmixr2)
23+
library(PopED)
24+
25+
f <- function() {
26+
ini({
27+
tp <- fix(100)
28+
td <- fix(0.001)
29+
te <- fix(1e-7)
30+
ts <- fix(20000)
31+
32+
tKA <- log(0.8)
33+
tKE <- log(0.15)
34+
tVD <- log(100)
35+
tEC50 <- log(0.12)
36+
tn <- log(2)
37+
tdelta <- log(0.2)
38+
tc <- log(7)
39+
40+
eta.KA ~ 0.25
41+
eta.KE ~ 0.25
42+
eta.VD ~ 0.25
43+
eta.EC50 ~ 0.25
44+
eta.n ~ 0.25
45+
eta.delta ~ 0.25
46+
eta.c ~ 0.25
47+
48+
add.sd.pk <- sqrt(0.04) # nlmixr2 uses sd
49+
add.sd.pd <- sqrt(0.04)
50+
})
51+
model({
52+
p <- tp
53+
d <- td
54+
e <- te
55+
s <- ts
56+
KA <- exp(tKA + eta.KA)
57+
KE <- exp(tKE + eta.KE)
58+
VD <- exp(tVD + eta.VD)
59+
EC50 <- exp(tEC50 + eta.EC50)
60+
n <- exp(tn + eta.n)
61+
delta <- exp(tdelta + eta.delta)
62+
c <- exp(tc + eta.c)
63+
64+
d/dt(depot) <- -KA*depot
65+
d/dt(central) <- KA*depot - KE*central
66+
d/dt(TC) <- s - TC*(e*VP + d) # target cells (TC)
67+
d/dt(IC) <- e*TC*VP - delta*IC # productively infected cells
68+
d/dt(VP) <- p*(1 - (pow(central/VD, n)/(pow(central/VD, n) +
69+
pow(EC50, n))))*IC - c*VP
70+
71+
## the two design variables enter here
72+
f(depot) <- DOSE # the dataset carries amt=1
73+
dur(depot) <- TINF # the dataset carries rate=-2
74+
75+
TC(0) <- c*delta/(p*e)
76+
IC(0) <- (s*e*p - d*c*delta)/(p*delta*e)
77+
VP(0) <- (s*e*p - d*c*delta)/(c*delta*e)
78+
79+
conc <- central/VD
80+
eff <- log10(VP)
81+
conc ~ add(add.sd.pk)
82+
eff ~ add(add.sd.pd)
83+
})
84+
}
85+
86+
TAU <- 7
87+
tms <- c(0, 0.25, 0.5, 1, 2, 3, 4, 7, 10, 14, 21, 28)
88+
89+
## note `amt=1` and `rate=-2`; the actual amount/duration come from the model
90+
e1 <- et(amt=1, rate=-2, ii=TAU, addl=3, cmt="depot") |>
91+
et(tms) |>
92+
as.data.frame()
93+
e1$dvid <- 1
94+
95+
e2 <- e1[e1$evid == 0, ]
96+
e2$dvid <- 2
97+
98+
e <- rbind(e1, e2)
99+
100+
babel.db <- nlmixr2(f, e, "poped",
101+
popedControl(groupsize=30,
102+
a=list(c(DOSE=180, TINF=1)),
103+
mina=c(DOSE=20, TINF=0.1),
104+
maxa=c(DOSE=600, TINF=12)))
105+
106+
plot_model_prediction(babel.db, facet_scales="free")
107+
108+
evaluate_design(babel.db)
109+
#> $ofv
110+
#> [1] 88.26974
111+
112+
## the criterion really does depend on the infusion duration
113+
for (ti in c(0.1, 0.5, 1, 2, 4, 8, 12)) {
114+
.db <- babel.db
115+
.db$design$a[1, "TINF"] <- ti
116+
message(sprintf("TINF = %5.2f ofv = %.4f", ti, evaluate_design(.db)$ofv))
117+
}
118+
119+
## ...so it can be optimized over
120+
r <- poped_optim(babel.db, opt_a=TRUE, opt_xt=FALSE, parallel=FALSE)
121+
#> Optimized Covariates:
122+
#> Group 1: 1 : 446.115 : 0.456544
123+
#>
124+
#> OFV = 95.5746
125+
#
126+
# (the default algorithm is a random search, so the exact DOSE/TINF the
127+
# optimizer lands on varies from run to run; the point is that both are
128+
# now optimized instead of being fixed features of the event table)
129+
130+
## Note: `cmt` in the design dataset must be given by *name*
131+
## (`cmt="depot"`); a bare compartment number is not translated in the
132+
## PopED path and silently produces a design with no dosing.

0 commit comments

Comments
 (0)