Skip to content

Commit

Permalink
Example 1 in babelmixr2
Browse files Browse the repository at this point in the history
  • Loading branch information
mattfidler committed Sep 12, 2024
1 parent 3128d19 commit d580ab8
Show file tree
Hide file tree
Showing 3 changed files with 204 additions and 3 deletions.
6 changes: 3 additions & 3 deletions inst/poped/ex.1.a.PK.1.comp.oral.md.intro.babelmixr2.R
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ evaluate_design(babel.db)
shrinkage(babel.db)

# Optimization of sample times
output <- poped_optim(babel.db, opt_xt =TRUE, parallel=TRUE)
output <- poped_optim(babel.db, opt_xt =TRUE)

# Evaluate optimization results
summary(output)
Expand All @@ -77,7 +77,7 @@ plot_model_prediction(output$poped.db)


# Optimization of sample times and doses
output_2 <- poped_optim(output$poped.db, opt_xt =TRUE, opt_a = TRUE, parallel = TRUE)
output_2 <- poped_optim(output$poped.db, opt_xt =TRUE, opt_a = TRUE)

summary(output_2)

Expand All @@ -90,7 +90,7 @@ plot_model_prediction(output_2$poped.db)
# faster than continuous optimization in this case
babel.db.discrete <- create.poped.database(babel.db,discrete_xt = list(0:248))

output_discrete <- poped_optim(babel.db.discrete, opt_xt=T, parallel = TRUE)
output_discrete <- poped_optim(babel.db.discrete, opt_xt=T)


summary(output_discrete)
Expand Down
99 changes: 99 additions & 0 deletions inst/poped/ex.1.b.PK.1.comp.oral.md.re-parameterize.babelmixr2.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,99 @@
## using libary models and reparameterizing the problen to KA, KE and V
## optimization of dose and dose interval

library(babelmixr2)

library(PopED)

f <- function() {
ini({
tV <- 72.8
tKa <- 0.25
tKe <- 3.75/72.8
tFavail <- fix(0.9)

eta.v ~ 0.09
eta.ka ~ 0.09
eta.ke ~ 0.25^2

prop.sd <- fix(sqrt(0.04))
add.sd <- fix(sqrt(5e-6))
})
model({
V <- tV*exp(eta.v)
KA <- tKa*exp(eta.ka)
KE <- tKe*exp(eta.ke)
Favail <- tFavail
N <- floor(time/TAU)+1
y <- (DOSE*Favail/V)*(KA/(KA - KE)) *
(exp(-KE * (time - (N - 1) * TAU)) * (1 - exp(-N * KE * TAU))/(1 - exp(-KE * TAU)) -
exp(-KA * (time - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU)))

y ~ prop(prop.sd) + add(add.sd)
})
}

# minxt, maxxt
e <- et(list(c(0, 10),
c(0, 10),
c(0, 10),
c(240, 248),
c(240, 248))) %>%
as.data.frame()

#xt
e$time <- c(1,2,8,240,245)


babel.db <- nlmixr2(f, e, "poped",
popedControl(groupsize=20,
bUseGrouped_xt=TRUE,
a=list(c(DOSE=20,TAU=24),
c(DOSE=40, TAU=24)),
maxa=c(DOSE=200,TAU=24),
mina=c(DOSE=0,TAU=24)))


## create plot of model without variability
plot_model_prediction(babel.db)

## create plot of model with variability
plot_model_prediction(babel.db,IPRED=T,DV=T,separate.groups=T)

## evaluate initial design
evaluate_design(babel.db)

shrinkage(babel.db)

# Optimization of sample times
output <- poped_optim(babel.db, opt_xt =TRUE, parallel=TRUE)

# Evaluate optimization results
summary(output)

get_rse(output$FIM,output$poped.db)

plot_model_prediction(output$poped.db)

# Optimization of sample times, doses and dose intervals
output_2 <- poped_optim(output$poped.db, opt_xt =TRUE, opt_a = TRUE, parallel = TRUE)

summary(output_2)
get_rse(output_2$FIM,output_2$poped.db)
plot_model_prediction(output_2$poped.db)

# Optimization of sample times with only integer time points in design space
# faster than continuous optimization in this case
babel.db.discrete <- create.poped.database(babel.db,discrete_xt = list(0:248))

output_discrete <- poped_optim(babel.db.discrete, opt_xt=T, parallel = TRUE)

summary(output_discrete)

get_rse(output_discrete$FIM,output_discrete$poped.db)

plot_model_prediction(output_discrete$poped.db)


# Efficiency of sampling windows
plot_efficiency_of_windows(output_discrete$poped.db, xt_windows=1)
102 changes: 102 additions & 0 deletions inst/poped/ex.1.c.PK.1.comp.oral.md.ODE.compiled.babelmixr2.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,102 @@
library(babelmixr2)
library(PopED)


## define the ODE
f <- function() {
ini({
tV <- 72.8
tKa <- 0.25
tCl <- 3.75
tF <- fix(0.9)

eta.v ~ 0.09
eta.ka ~ 0.09
eta.cl ~0.25^2

prop.sd <- fix(sqrt(0.04))
add.sd <- fix(sqrt(5e-6))

})
model({
V<-tV*exp(eta.v)
KA<-tKa*exp(eta.ka)
CL<-tCl*exp(eta.cl)
Favail <- tF
d/dt(depot) <- -KA*depot
d/dt(central) <- KA*depot - (CL/V)*central
f(depot) <- Favail*DOSE
y <- central/V
y ~ prop(prop.sd) + add(add.sd)
})
}

# minxt, maxxt
e <- et(list(c(0, 10),
c(0, 10),
c(0, 10),
c(240, 248),
c(240, 248))) %>%
et(amt=1000, ii=24, until=248,cmt="depot") %>%
as.data.frame()

#xt
e$time <- c(0, 1,2,8,240,245)


babel.db <- nlmixr2(f, e, "poped",
popedControl(groupsize=20,
bUseGrouped_xt=TRUE,
a=list(c(DOSE=20,TAU=24),
c(DOSE=40, TAU=24)),
maxa=c(DOSE=200,TAU=24),
mina=c(DOSE=0,TAU=24)))

## create plot of model without variability
plot_model_prediction(babel.db, model_num_points = 300)

## create plot of model with variability
plot_model_prediction(babel.db, IPRED=T, DV=T, separate.groups=T, model_num_points = 300)

## evaluate initial design
evaluate_design(babel.db)

shrinkage(babel.db)

# Optimization of sample times
output <- poped_optim(babel.db, opt_xt =TRUE)

# Evaluate optimization results
summary(output)

get_rse(output$FIM,output$poped.db)

plot_model_prediction(output$poped.db)


# Optimization of sample times and doses
output_2 <- poped_optim(output$poped.db, opt_xt =TRUE, opt_a = TRUE)

summary(output_2)

get_rse(output_2$FIM,output_2$poped.db)

plot_model_prediction(output_2$poped.db)


# Optimization of sample times with only integer time points in design space
# faster than continuous optimization in this case
babel.db.discrete <- create.poped.database(babel.db,discrete_xt = list(0:248))

output_discrete <- poped_optim(babel.db.discrete, opt_xt=T)


summary(output_discrete)

get_rse(output_discrete$FIM,output_discrete$poped.db)

plot_model_prediction(output_discrete$poped.db)


# Efficiency of sampling windows
plot_efficiency_of_windows(output_discrete$poped.db, xt_windows=1)

0 comments on commit d580ab8

Please sign in to comment.