kylebaron
commented
2 months ago Hi @kevin199011 -
Thanks for the question and for putting together the reproducible example. The (big) discrepancy in simulation 2 is due to the next observation carry backward behavior which is similar to how NONMEM works. You can control this behavior in mrgsolve with the nocb argument (see example simulation below). I'd say any time you have time-varying covariates, the simulation result will depend on how you construct the data set.
I'd have to study the inputs more to explain the (smaller) discrepancy with simulation 1.
Also, it's possible that we can make changes to the code so that you don't have this dependence on the data set (I'm not sure ... would have to play with it a bit).
But I'm thinking that the outputs are as expected.
Does that clarify things for you?
Kyle
library(tidyverse)
library(mrgsolve)
#>
#> Attaching package: 'mrgsolve'
#> The following object is masked from 'package:stats':
#>
#> filter
code1 <- '
[PARAM]
ALPHA = 0.1
BETA = 0.1
DAR0 = 8
TDOSE = 0
[CMT]
ADC PAYLOAD
[MAIN]
ADC_0 = 10;
[ODE]
double TAD = SOLVERTIME - TDOSE;
double DAR = DAR0*exp(-BETA*TAD);
dxdt_ADC = -ALPHA * ADC;
dxdt_PAYLOAD = DAR*ADC;
[CAPTURE]
TDOSE TAD DAR
'
mod1 <- mcode("time_variant1", code1)
#> Building time_variant1 ...
#> done.
ev <- ev(ID = 1, amt = 0, time = c(0,12)) %>% as.data.frame
obs1 <- expand.grid(ID = 1,
time = c(0:24),
amt = NA, cmt = 1, evid = 0)
obs2 <- expand.grid(ID = 1,
time = c(0,6,12,18,24),
amt = NA, cmt = 1, evid = 0)
df1 <- bind_rows(ev, obs1) %>%
arrange(ID, time,desc(evid)) %>%
group_by(ID) %>%
mutate(TDOSE = ifelse(evid > 0, time, NA)) %>%
fill(TDOSE)
df2 <- bind_rows(ev, obs2) %>%
arrange(ID, time,desc(evid)) %>%
group_by(ID) %>%
mutate(TDOSE = ifelse(evid > 0, time, NA)) %>%
fill(TDOSE)
df3 <- bind_rows(ev, obs1) %>%
arrange(ID, time, evid) %>%
group_by(ID) %>%
mutate(TDOSE = ifelse(evid > 0, time, NA)) %>%
fill(TDOSE, .direction = "downup")
df4 <- bind_rows(ev, obs2) %>%
arrange(ID, time, evid) %>%
group_by(ID) %>%
mutate(TDOSE = ifelse(evid > 0, time, NA)) %>%
fill(TDOSE, .direction = "downup")
out1 <- mod1 %>%
data_set(df1) %>%
mrgsim(carry.out = c("evid")) %>%
as.data.frame() %>%
select(ID, time, evid, PAYLOAD) %>%
rename(PAYLOAD1 = PAYLOAD)
out2 <- mod1 %>%
data_set(df2) %>%
mrgsim(carry.out = c("evid"), nocb = FALSE) %>%
as.data.frame()%>%
select(ID, time, evid, PAYLOAD) %>%
rename(PAYLOAD2 = PAYLOAD)
out3 <- mod1 %>%
data_set(df3) %>%
mrgsim(carry.out = c("evid")) %>%
as.data.frame()%>%
select(ID, time, evid, PAYLOAD) %>%
rename(PAYLOAD3 = PAYLOAD)
out4 <- mod1 %>%
data_set(df4) %>%
mrgsim(carry.out = c("evid")) %>%
as.data.frame()%>%
select(ID, time, evid, PAYLOAD) %>%
rename(PAYLOAD4 = PAYLOAD)
out <- left_join(out1, out2) %>%
left_join(out3) %>%
left_join(out4)
#> Joining with `by = join_by(ID, time, evid)`
#> Joining with `by = join_by(ID, time, evid)`
#> Joining with `by = join_by(ID, time, evid)`
knitr::kable(out)| ID | time | evid | PAYLOAD1 | PAYLOAD2 | PAYLOAD3 | PAYLOAD4 |
|---|---|---|---|---|---|---|
| 1 | 0 | 1 | 0.0000 | 0.0000 | 0.0000 | 0.0000 |
| 1 | 0 | 0 | 0.0000 | 0.0000 | 0.0000 | 0.0000 |
| 1 | 1 | 0 | 72.5077 | NA | 72.5077 | NA |
| 1 | 2 | 0 | 131.8720 | NA | 131.8720 | NA |
| 1 | 3 | 0 | 180.4753 | NA | 180.4753 | NA |
| 1 | 4 | 0 | 220.2684 | NA | 220.2684 | NA |
| 1 | 5 | 0 | 252.8482 | NA | 252.8482 | NA |
| 1 | 6 | 0 | 279.5223 | 279.5223 | 279.5223 | 279.5223 |
| 1 | 7 | 0 | 301.3612 | NA | 301.3612 | NA |
| 1 | 8 | 0 | 319.2414 | NA | 319.2414 | NA |
| 1 | 9 | 0 | 333.8804 | NA | 333.8804 | NA |
| 1 | 10 | 0 | 345.8659 | NA | 345.8659 | NA |
| 1 | 11 | 0 | 355.6787 | NA | 355.6787 | NA |
| 1 | 12 | 1 | 382.3528 | 363.7128 | 363.7128 | 363.7128 |
| 1 | 12 | 0 | 382.3528 | 363.7128 | 363.7128 | 363.7128 |
| 1 | 13 | 0 | 404.1917 | NA | 385.5517 | NA |
| 1 | 14 | 0 | 422.0719 | NA | 403.4319 | NA |
| 1 | 15 | 0 | 436.7109 | NA | 418.0709 | NA |
| 1 | 16 | 0 | 448.6964 | NA | 430.0564 | NA |
| 1 | 17 | 0 | 458.5092 | NA | 439.8692 | NA |
| 1 | 18 | 0 | 466.5433 | 447.9033 | 447.9033 | 447.9033 |
| 1 | 19 | 0 | 473.1211 | NA | 454.4811 | NA |
| 1 | 20 | 0 | 478.5065 | NA | 459.8665 | NA |
| 1 | 21 | 0 | 482.9157 | NA | 464.2757 | NA |
| 1 | 22 | 0 | 486.5256 | NA | 467.8856 | NA |
| 1 | 23 | 0 | 489.4812 | NA | 470.8412 | NA |
| 1 | 24 | 0 | 491.9010 | 473.2610 | 473.2610 | 473.2610 |
Created on 2024-06-21 with reprex v2.0.2
kevin199011
Hi, I'm working on an ADC model, involving modeling the time-varying DAR. Below reprex shows a simplified mrgsolve code, where DAR is expressed by first-order decay, and PAYLOAD is expressed by DAR*ADC. When I was doing the simulation, I noticed how the simulation dataset was constructed has an impact on the simulation result. Please see below example. When there was an observation record with time equals dosing time, the order of this record seems to influence the result. In below example, there was a second dose at 12.
As can be seen from the output, df1 and df2 simulated very different results from time=12, while df3 and df4 produced same output, while slightly than df1.
I'm wondering if you can help explain why it has this behavior, and if it's the correct way to code with SOLVERTIME this way? Thank you!
Created on 2024-06-21 by the reprex package (v2.0.0)