From 7a36270fe4f5a531f1327d0934c19a7a9d01da76 Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Fri, 4 Oct 2024 03:25:44 +0200 Subject: [PATCH 01/24] Examples added --- inst/poped/ex.14.PK.IOV.babelmixr2.R | 88 +++++++++ .../ex.15.full.covariance.matrix.babelmixr2.R | 120 ++++++++++++ .../poped/ex.6.PK.1.comp.oral.sd.babelmixr2.R | 175 ++++++++++++++++++ 3 files changed, 383 insertions(+) create mode 100644 inst/poped/ex.14.PK.IOV.babelmixr2.R create mode 100644 inst/poped/ex.15.full.covariance.matrix.babelmixr2.R create mode 100644 inst/poped/ex.6.PK.1.comp.oral.sd.babelmixr2.R diff --git a/inst/poped/ex.14.PK.IOV.babelmixr2.R b/inst/poped/ex.14.PK.IOV.babelmixr2.R new file mode 100644 index 00000000..64cd5445 --- /dev/null +++ b/inst/poped/ex.14.PK.IOV.babelmixr2.R @@ -0,0 +1,88 @@ +library(babelmixr2) +library(PopED) +library(ggplot2) + +# +# N = floor(Time/TAU)+1; +# CL = CL_OCC_1; +# if(N>6) CL = CL_OCC_2; + +## 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 + + # IOV + theta.iov.cl <- c(0.09) + + iov.cl1 <- fix(1) + # iov.cl2 <- fix(1) + # iov.cl3 <- fix (1) # as many as occasions as needed + + + prop.sd <- fix(sqrt(0.04)) + add.sd <- fix(sqrt(5e-6)) + + }) + model({ + # No of doses + N = floor(t/TAU)+1 + OCC1 = ifelse(N>6, 1, 0) + iov <- theta.iov.cl * (OCC1 * iov.cl1) + + V<-tV*exp(eta.v) + KA<-tKa*exp(eta.ka) + CL = tCl *exp(eta.cl + iov) + 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) + }) +} + +# careful with the et. the amt is scaled to the F +# minxt, maxxt +e <- et(list(c(0, 10), + c(0, 10), + c(0, 10), + c(240, 248), + c(240, 248))) %>% + et(amt=1/0.9, 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) + coord_cartesian(ylim = c(0,0.5)) + +## Visulaize the IOV +set.seed(12345679) +plot_model_prediction(babel.db, PRED=F,IPRED=F, + separate.groups=T, model_num_points = 500, + groupsize_sim = 1, + IPRED.lines = T, alpha.IPRED.lines=0.6, + sample.times = F +) + geom_vline(xintercept = 24*6,color="red") + coord_cartesian(ylim = c(0, 0.75)) + +## evaluate initial design +evaluate_design(babel.db) + diff --git a/inst/poped/ex.15.full.covariance.matrix.babelmixr2.R b/inst/poped/ex.15.full.covariance.matrix.babelmixr2.R new file mode 100644 index 00000000..ff0470fa --- /dev/null +++ b/inst/poped/ex.15.full.covariance.matrix.babelmixr2.R @@ -0,0 +1,120 @@ +## Warfarin example from software comparison in: +## Nyberg et al., "Methods and software tools for design evaluation +## for population pharmacokinetics-pharmacodynamics studies", +## Br. J. Clin. Pharm., 2014. + +library(PopED) +# library(deSolve) +library(babelmixr2) + +##-- Model: One comp first order absorption +##-- Model: One comp first order absorption, analytic solution +f_without <- function() { + ini({ + tV <- 8 + tKa <- 1 + tCl <- 0.15 + tF <- fix(1) + + eta.v ~ 0.02 + eta.ka ~ 0.6 + eta.cl ~0.07 + + prop.sd <- sqrt(0.01) + + }) + model({ + V<-tV*exp(eta.v) + KA<-tKa*exp(eta.ka) + CL<-tCl*exp(eta.cl) + + KE=CL/V + xt=t # set xt as time + Favail <- tF + + y <- (DOSE*Favail*KA/(V*(KA-KE)))*(exp(-KE*xt)-exp(-KA*xt)) + + y ~ prop(prop.sd) + }) +} + +f_with <- function() { + ini({ + tV <- 8 + tKa <- 1 + tCl <- 0.15 + tF <- fix(1) + + ## For correlated parameters, you specify the names of each + ## correlated parameter separated by a addition operator `+` + ## and the left handed side specifies the lower triangular + ## matrix initial of the covariance matrix. + eta.cl + eta.v + eta.ka ~ c(0.07, + 0.03, 0.02, + 0.1, 0.09, 0.6) + + prop.sd <- sqrt(0.01) + + }) + model({ + V<-tV*exp(eta.v) + KA<-tKa*exp(eta.ka) + CL<-tCl*exp(eta.cl) + + KE=CL/V + xt=t # set xt as time + Favail <- tF + + y <- (DOSE*Favail*KA/(V*(KA-KE)))*(exp(-KE*xt)-exp(-KA*xt)) + + y ~ prop(prop.sd) + }) +} + +# minxt, maxxt +e <- et(list(c(0, 120), + c(0, 120), + c(0, 120), + c(0, 120), + c(0, 120), + c(0, 120), + c(0, 120), + c(0, 120))) %>% + as.data.frame() + +#xt +e$time <- c(0.5,1,2,6,24,36,72,120) + +babel.db.without <- nlmixr2(f_without, e, "poped", + popedControl(groupsize=32, + bUseGrouped_xt=TRUE, + a=c(DOSE=70), + maxa=c(DOSE=200), + mina=c(DOSE=0))) + +babel.db.with <- nlmixr2(f_with, e, "poped", + popedControl(groupsize=32, + bUseGrouped_xt=TRUE, + a=c(DOSE=70), + maxa=c(DOSE=200), + mina=c(DOSE=0))) + + +# What do the covariances mean? +(IIV <- babel.db.with$parameters$param.pt.val$d) +cov2cor(IIV) + +## create plot of model with variability +library(ggplot2) +p1 <- plot_model_prediction(babel.db.without,IPRED=T)+ylim(0,12) +p2 <- plot_model_prediction(babel.db.with,IPRED=T) +ylim(0,12) +gridExtra::grid.arrange(p1, p2, nrow = 1) + +## evaluate initial design +evaluate_design(babel.db.without) +evaluate_design(babel.db.with) + +## Evaluate with full FIM +evaluate_design(babel.db.without, fim.calc.type=0) +evaluate_design(babel.db.with, fim.calc.type=0) + diff --git a/inst/poped/ex.6.PK.1.comp.oral.sd.babelmixr2.R b/inst/poped/ex.6.PK.1.comp.oral.sd.babelmixr2.R new file mode 100644 index 00000000..f1133232 --- /dev/null +++ b/inst/poped/ex.6.PK.1.comp.oral.sd.babelmixr2.R @@ -0,0 +1,175 @@ +library(PopED) +# library(deSolve) +library(babelmixr2) + +##-- Model: One comp first order absorption +##-- Model: One comp first order absorption, analytic solution +f1 <- function() { + ini({ + tV <- 8 + tKa <- 1 + tCl <- 0.15 + tF <- fix(1) + + eta.v ~ 0.02 + eta.ka ~ 0.6 + eta.cl ~0.07 + + prop.sd <- sqrt(0.01) + add.sd <- sqrt(0.25) + + }) + model({ + V<-tV*exp(eta.v) + KA<-tKa*exp(eta.ka) + CL<-tCl*exp(eta.cl) + + KE=CL/V + xt=t # set xt as time + Favail <- tF + + y <- (DOSE*Favail*KA/(V*(KA-KE)))*(exp(-KE*xt)-exp(-KA*xt)) + + y ~ prop(prop.sd) + add(add.sd) + }) +} + +##-- Model: One comp first order absorption, analytic solution +## reparameterization +f2 <- function() { + ini({ + tV <- 8 + tKa <- 1 + tKe <- 0.15/8 + tF <- fix(1) + + eta.v ~ 0.02 + eta.ka ~ 0.6 + eta.ke ~0.07 + + prop.sd <- sqrt(0.01) + add.sd <- sqrt(0.25) + + }) + model({ + V<-tV*exp(eta.v) + KA<-tKa*exp(eta.ka) + KE=tKe*exp(eta.ke) + + Favail <- tF + xt=t # set xt as time + + y <- (DOSE*Favail*KA/(V*(KA-KE)))*(exp(-KE*xt)-exp(-KA*xt)) + + y ~ prop(prop.sd) + add(add.sd) + }) +} + +# ODE solution +f3 <- function() { + ini({ + tV <- 8 + tKa <- 1 + tCl <- 0.15 + tF <- fix(1) + + eta.v ~ 0.02 + eta.ka ~ 0.6 + eta.cl ~0.07 + + prop.sd <- sqrt(0.01) + add.sd <- sqrt(0.25) + + }) + 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, 25), + c(0, 25), + c(0, 25), + c(0, 120), + c(0, 120), + c(0, 120), + c(0, 120), + c(0, 120))) %>% + # users must be careful to not add a dosing record + # et(amt=70, cmt="depot") %>% + as.data.frame() + +#xt +e$time <- c(1,2,3,6,24,36,72,120) + +babel.db.1 <- nlmixr2(f1, e, "poped", + popedControl(groupsize=32, + bUseGrouped_xt=TRUE, + a=c(DOSE=70), + maxa=c(DOSE=200), + mina=c(DOSE=0))) + + +babel.db.2 <- nlmixr2(f2, e, "poped", + popedControl(groupsize=32, + bUseGrouped_xt=TRUE, + a=c(DOSE=70), + maxa=c(DOSE=200), + mina=c(DOSE=0))) + + +# Add a dosing event as placeholder +e2 <- et(list(c(0, 25), + c(0, 25), + c(0, 25), + c(0, 120), + c(0, 120), + c(0, 120), + c(0, 120), + c(0, 120))) %>% + et(amt=1, cmt="depot") %>% + as.data.frame() + +babel.db.3 <- nlmixr2(f3, e2, "poped", + popedControl(groupsize=32, + bUseGrouped_xt=TRUE, + a=c(DOSE=70), + maxa=c(DOSE=200), + mina=c(DOSE=0))) + + +## create plot of models +(plot1 <- plot_model_prediction(babel.db.1,IPRED=T)) +(plot2 <- plot_model_prediction(babel.db.2,IPRED=T)) +(plot3 <- plot_model_prediction(babel.db.3,IPRED=T)) + +# different results for different parameterizations +library(ggplot2) +library(gridExtra) +plot1 <- plot1 + ggtitle("CL Parameterization") +plot2 <- plot2 + ggtitle("KE Parameterization") +grid.arrange(plot1,plot2) + + +## evaluate initial designs +# different results for different parameterizations +evaluate_design(babel.db.1) +evaluate_design(babel.db.2) +evaluate_design(babel.db.3) + + +# Optimization of sample times +# different results for different parameterizations +output.1 <- poped_optim(poped.db.1,opt_xt=T,parallel=T) +output.2 <- poped_optim(poped.db.2,opt_xt=T,parallel=T) + +summary(output.1) +summary(output.2) From fdd6a819740ca49522320c36bd5deb1442811d72 Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Wed, 16 Oct 2024 10:50:16 +0200 Subject: [PATCH 02/24] example file deletions --- R/writeLines.R | 52 ++++++++ inst/poped/ex.14.PK.IOV.babelmixr2.R | 88 ------------- .../ex.15.full.covariance.matrix.babelmixr2.R | 120 ------------------ 3 files changed, 52 insertions(+), 208 deletions(-) create mode 100644 R/writeLines.R delete mode 100644 inst/poped/ex.14.PK.IOV.babelmixr2.R delete mode 100644 inst/poped/ex.15.full.covariance.matrix.babelmixr2.R diff --git a/R/writeLines.R b/R/writeLines.R new file mode 100644 index 00000000..4b7e64dc --- /dev/null +++ b/R/writeLines.R @@ -0,0 +1,52 @@ +# Function to write tree scructure for folders + +# Install fs package if not already installed +if (!requireNamespace("fs", quietly = TRUE)) { + install.packages("fs") +} + +# Load the fs package + + +# Define a function to get and save the folder tree with symbols +save_folder_tree_with_symbols <- function(path = ".", output_file = "folder_tree.txt") { + library(fs) + # Use fs::dir_tree to create the folder tree and get individual lines + folder_tree <- dir_tree(path, recurse = TRUE) + + # Define a function to add symbols based on depth + add_symbols <- function(lines) { + result <- c() + depth <- 0 + + for (i in seq_along(lines)) { + # Calculate the depth based on the number of slashes in the path + depth <- sum(strsplit(lines[i], "/")[[1]] != "") + + # Determine if it's a file or folder and if it's the last in its group + is_last <- (i == length(lines)) || + (sum(strsplit(lines[i + 1], "/")[[1]] != "") <= depth) + + # Add appropriate symbols for the depth level + prefix <- ifelse(depth > 1, paste0(strrep("│ ", depth - 2), ifelse(is_last, "└─ ", "├─ ")), "") + + result <- c(result, paste0(prefix, basename(lines[i]))) + } + return(result) + } + + # Create a visual representation of the folder tree with symbols + visual_tree <- add_symbols(folder_tree) + + # Print the visual folder tree to the console + cat(visual_tree, sep = "\n") + + # Save the visual folder tree to a text file + writeLines(visual_tree, output_file) +} + +# Call the function on the + +save_folder_tree_with_symbols(path = here::here("SCRIPTS"), output_file = "folder_tree_SCRIPTS.txt") +save_folder_tree_with_symbols(path = here::here("ANALYSES"), output_file = "folder_tree_ANALYSES.txt") +save_folder_tree_with_symbols(path = here::here("DATASETS"), output_file = "folder_tree_DATASETS.txt") diff --git a/inst/poped/ex.14.PK.IOV.babelmixr2.R b/inst/poped/ex.14.PK.IOV.babelmixr2.R deleted file mode 100644 index 64cd5445..00000000 --- a/inst/poped/ex.14.PK.IOV.babelmixr2.R +++ /dev/null @@ -1,88 +0,0 @@ -library(babelmixr2) -library(PopED) -library(ggplot2) - -# -# N = floor(Time/TAU)+1; -# CL = CL_OCC_1; -# if(N>6) CL = CL_OCC_2; - -## 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 - - # IOV - theta.iov.cl <- c(0.09) - - iov.cl1 <- fix(1) - # iov.cl2 <- fix(1) - # iov.cl3 <- fix (1) # as many as occasions as needed - - - prop.sd <- fix(sqrt(0.04)) - add.sd <- fix(sqrt(5e-6)) - - }) - model({ - # No of doses - N = floor(t/TAU)+1 - OCC1 = ifelse(N>6, 1, 0) - iov <- theta.iov.cl * (OCC1 * iov.cl1) - - V<-tV*exp(eta.v) - KA<-tKa*exp(eta.ka) - CL = tCl *exp(eta.cl + iov) - 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) - }) -} - -# careful with the et. the amt is scaled to the F -# minxt, maxxt -e <- et(list(c(0, 10), - c(0, 10), - c(0, 10), - c(240, 248), - c(240, 248))) %>% - et(amt=1/0.9, 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) + coord_cartesian(ylim = c(0,0.5)) - -## Visulaize the IOV -set.seed(12345679) -plot_model_prediction(babel.db, PRED=F,IPRED=F, - separate.groups=T, model_num_points = 500, - groupsize_sim = 1, - IPRED.lines = T, alpha.IPRED.lines=0.6, - sample.times = F -) + geom_vline(xintercept = 24*6,color="red") + coord_cartesian(ylim = c(0, 0.75)) - -## evaluate initial design -evaluate_design(babel.db) - diff --git a/inst/poped/ex.15.full.covariance.matrix.babelmixr2.R b/inst/poped/ex.15.full.covariance.matrix.babelmixr2.R deleted file mode 100644 index ff0470fa..00000000 --- a/inst/poped/ex.15.full.covariance.matrix.babelmixr2.R +++ /dev/null @@ -1,120 +0,0 @@ -## Warfarin example from software comparison in: -## Nyberg et al., "Methods and software tools for design evaluation -## for population pharmacokinetics-pharmacodynamics studies", -## Br. J. Clin. Pharm., 2014. - -library(PopED) -# library(deSolve) -library(babelmixr2) - -##-- Model: One comp first order absorption -##-- Model: One comp first order absorption, analytic solution -f_without <- function() { - ini({ - tV <- 8 - tKa <- 1 - tCl <- 0.15 - tF <- fix(1) - - eta.v ~ 0.02 - eta.ka ~ 0.6 - eta.cl ~0.07 - - prop.sd <- sqrt(0.01) - - }) - model({ - V<-tV*exp(eta.v) - KA<-tKa*exp(eta.ka) - CL<-tCl*exp(eta.cl) - - KE=CL/V - xt=t # set xt as time - Favail <- tF - - y <- (DOSE*Favail*KA/(V*(KA-KE)))*(exp(-KE*xt)-exp(-KA*xt)) - - y ~ prop(prop.sd) - }) -} - -f_with <- function() { - ini({ - tV <- 8 - tKa <- 1 - tCl <- 0.15 - tF <- fix(1) - - ## For correlated parameters, you specify the names of each - ## correlated parameter separated by a addition operator `+` - ## and the left handed side specifies the lower triangular - ## matrix initial of the covariance matrix. - eta.cl + eta.v + eta.ka ~ c(0.07, - 0.03, 0.02, - 0.1, 0.09, 0.6) - - prop.sd <- sqrt(0.01) - - }) - model({ - V<-tV*exp(eta.v) - KA<-tKa*exp(eta.ka) - CL<-tCl*exp(eta.cl) - - KE=CL/V - xt=t # set xt as time - Favail <- tF - - y <- (DOSE*Favail*KA/(V*(KA-KE)))*(exp(-KE*xt)-exp(-KA*xt)) - - y ~ prop(prop.sd) - }) -} - -# minxt, maxxt -e <- et(list(c(0, 120), - c(0, 120), - c(0, 120), - c(0, 120), - c(0, 120), - c(0, 120), - c(0, 120), - c(0, 120))) %>% - as.data.frame() - -#xt -e$time <- c(0.5,1,2,6,24,36,72,120) - -babel.db.without <- nlmixr2(f_without, e, "poped", - popedControl(groupsize=32, - bUseGrouped_xt=TRUE, - a=c(DOSE=70), - maxa=c(DOSE=200), - mina=c(DOSE=0))) - -babel.db.with <- nlmixr2(f_with, e, "poped", - popedControl(groupsize=32, - bUseGrouped_xt=TRUE, - a=c(DOSE=70), - maxa=c(DOSE=200), - mina=c(DOSE=0))) - - -# What do the covariances mean? -(IIV <- babel.db.with$parameters$param.pt.val$d) -cov2cor(IIV) - -## create plot of model with variability -library(ggplot2) -p1 <- plot_model_prediction(babel.db.without,IPRED=T)+ylim(0,12) -p2 <- plot_model_prediction(babel.db.with,IPRED=T) +ylim(0,12) -gridExtra::grid.arrange(p1, p2, nrow = 1) - -## evaluate initial design -evaluate_design(babel.db.without) -evaluate_design(babel.db.with) - -## Evaluate with full FIM -evaluate_design(babel.db.without, fim.calc.type=0) -evaluate_design(babel.db.with, fim.calc.type=0) - From bacd64be757f04e39682c767ffdfc095722aa97c Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Wed, 16 Oct 2024 10:53:36 +0200 Subject: [PATCH 03/24] added examples 14 and 15 --- inst/poped/ex.14.PK.IOV.babelmixr2.R | 88 +++++++++++++ .../ex.15.full.covariance.matrix.babelmixr2.R | 120 ++++++++++++++++++ 2 files changed, 208 insertions(+) create mode 100644 inst/poped/ex.14.PK.IOV.babelmixr2.R create mode 100644 inst/poped/ex.15.full.covariance.matrix.babelmixr2.R diff --git a/inst/poped/ex.14.PK.IOV.babelmixr2.R b/inst/poped/ex.14.PK.IOV.babelmixr2.R new file mode 100644 index 00000000..64cd5445 --- /dev/null +++ b/inst/poped/ex.14.PK.IOV.babelmixr2.R @@ -0,0 +1,88 @@ +library(babelmixr2) +library(PopED) +library(ggplot2) + +# +# N = floor(Time/TAU)+1; +# CL = CL_OCC_1; +# if(N>6) CL = CL_OCC_2; + +## 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 + + # IOV + theta.iov.cl <- c(0.09) + + iov.cl1 <- fix(1) + # iov.cl2 <- fix(1) + # iov.cl3 <- fix (1) # as many as occasions as needed + + + prop.sd <- fix(sqrt(0.04)) + add.sd <- fix(sqrt(5e-6)) + + }) + model({ + # No of doses + N = floor(t/TAU)+1 + OCC1 = ifelse(N>6, 1, 0) + iov <- theta.iov.cl * (OCC1 * iov.cl1) + + V<-tV*exp(eta.v) + KA<-tKa*exp(eta.ka) + CL = tCl *exp(eta.cl + iov) + 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) + }) +} + +# careful with the et. the amt is scaled to the F +# minxt, maxxt +e <- et(list(c(0, 10), + c(0, 10), + c(0, 10), + c(240, 248), + c(240, 248))) %>% + et(amt=1/0.9, 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) + coord_cartesian(ylim = c(0,0.5)) + +## Visulaize the IOV +set.seed(12345679) +plot_model_prediction(babel.db, PRED=F,IPRED=F, + separate.groups=T, model_num_points = 500, + groupsize_sim = 1, + IPRED.lines = T, alpha.IPRED.lines=0.6, + sample.times = F +) + geom_vline(xintercept = 24*6,color="red") + coord_cartesian(ylim = c(0, 0.75)) + +## evaluate initial design +evaluate_design(babel.db) + diff --git a/inst/poped/ex.15.full.covariance.matrix.babelmixr2.R b/inst/poped/ex.15.full.covariance.matrix.babelmixr2.R new file mode 100644 index 00000000..ff0470fa --- /dev/null +++ b/inst/poped/ex.15.full.covariance.matrix.babelmixr2.R @@ -0,0 +1,120 @@ +## Warfarin example from software comparison in: +## Nyberg et al., "Methods and software tools for design evaluation +## for population pharmacokinetics-pharmacodynamics studies", +## Br. J. Clin. Pharm., 2014. + +library(PopED) +# library(deSolve) +library(babelmixr2) + +##-- Model: One comp first order absorption +##-- Model: One comp first order absorption, analytic solution +f_without <- function() { + ini({ + tV <- 8 + tKa <- 1 + tCl <- 0.15 + tF <- fix(1) + + eta.v ~ 0.02 + eta.ka ~ 0.6 + eta.cl ~0.07 + + prop.sd <- sqrt(0.01) + + }) + model({ + V<-tV*exp(eta.v) + KA<-tKa*exp(eta.ka) + CL<-tCl*exp(eta.cl) + + KE=CL/V + xt=t # set xt as time + Favail <- tF + + y <- (DOSE*Favail*KA/(V*(KA-KE)))*(exp(-KE*xt)-exp(-KA*xt)) + + y ~ prop(prop.sd) + }) +} + +f_with <- function() { + ini({ + tV <- 8 + tKa <- 1 + tCl <- 0.15 + tF <- fix(1) + + ## For correlated parameters, you specify the names of each + ## correlated parameter separated by a addition operator `+` + ## and the left handed side specifies the lower triangular + ## matrix initial of the covariance matrix. + eta.cl + eta.v + eta.ka ~ c(0.07, + 0.03, 0.02, + 0.1, 0.09, 0.6) + + prop.sd <- sqrt(0.01) + + }) + model({ + V<-tV*exp(eta.v) + KA<-tKa*exp(eta.ka) + CL<-tCl*exp(eta.cl) + + KE=CL/V + xt=t # set xt as time + Favail <- tF + + y <- (DOSE*Favail*KA/(V*(KA-KE)))*(exp(-KE*xt)-exp(-KA*xt)) + + y ~ prop(prop.sd) + }) +} + +# minxt, maxxt +e <- et(list(c(0, 120), + c(0, 120), + c(0, 120), + c(0, 120), + c(0, 120), + c(0, 120), + c(0, 120), + c(0, 120))) %>% + as.data.frame() + +#xt +e$time <- c(0.5,1,2,6,24,36,72,120) + +babel.db.without <- nlmixr2(f_without, e, "poped", + popedControl(groupsize=32, + bUseGrouped_xt=TRUE, + a=c(DOSE=70), + maxa=c(DOSE=200), + mina=c(DOSE=0))) + +babel.db.with <- nlmixr2(f_with, e, "poped", + popedControl(groupsize=32, + bUseGrouped_xt=TRUE, + a=c(DOSE=70), + maxa=c(DOSE=200), + mina=c(DOSE=0))) + + +# What do the covariances mean? +(IIV <- babel.db.with$parameters$param.pt.val$d) +cov2cor(IIV) + +## create plot of model with variability +library(ggplot2) +p1 <- plot_model_prediction(babel.db.without,IPRED=T)+ylim(0,12) +p2 <- plot_model_prediction(babel.db.with,IPRED=T) +ylim(0,12) +gridExtra::grid.arrange(p1, p2, nrow = 1) + +## evaluate initial design +evaluate_design(babel.db.without) +evaluate_design(babel.db.with) + +## Evaluate with full FIM +evaluate_design(babel.db.without, fim.calc.type=0) +evaluate_design(babel.db.with, fim.calc.type=0) + From 30885d0ebb68f4bd4fe02da1880c5e51871c0c42 Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Wed, 16 Oct 2024 20:05:11 +0200 Subject: [PATCH 04/24] add poped example 12 --- ...ex.12.covariate.distributions.babelmixr2.R | 104 ++++++++++++++++++ 1 file changed, 104 insertions(+) create mode 100644 inst/poped/ex.12.covariate.distributions.babelmixr2.R diff --git a/inst/poped/ex.12.covariate.distributions.babelmixr2.R b/inst/poped/ex.12.covariate.distributions.babelmixr2.R new file mode 100644 index 00000000..b2d68728 --- /dev/null +++ b/inst/poped/ex.12.covariate.distributions.babelmixr2.R @@ -0,0 +1,104 @@ +rm(list=ls()) +library(PopED) +library() +## Introduction +# Lets assume that we have a model with a covariate included in the model description. +# Here we define a one-compartment PK model that has weight on both clearance and volume of distribution. + +f <- function() { + ini({ + tCl <- 3.8 + tV <- 20 + WT_CL <- 0.75 + WT_V <- 1 + + eta.cl ~ 0.05 + eta.v ~ 0.05 + + # Residual unexplained variability (RUV) + add.sd <- fix(sqrt(0.0015)) # Additive error + prop.sd <- sqrt(0.015) # Proportional error + + }) + model({ + + CL <- tCl * (WT/70)^(WT_CL) * exp(eta.cl) + V <- tV * (WT/70)^(WT_V) * exp(eta.v) + + DOSE = 1000*(WT/70) + + y = DOSE/V*exp(-CL/V*time) + y ~ prop(prop.sd) + add(add.sd) + + }) +} + +e <- et(c( 1,2,4,6,8,24)) %>% + as.data.frame() + +## -- Define initial design and design space +babel.db <- nlmixr2(f, e, "poped", + control=popedControl( + groupsize=50, + minxt=0, + maxxt=24, + bUseGrouped_xt = T, + a=c(WT=70))) + +# We can create a plot of the model typical predictions: +plot_model_prediction(poped_db) + +# And evaluate the initial design +evaluate_design(poped_db) + +# We see that the covariate parameters can not be estimated +# according to this design calculation (RSE of bpop[3]=0 and bpop[4]=0). +# Why is that? Well, the calculation being done is assuming that every +# individual in the group has the same covariate (to speed +# up the calculation). This is clearly a poor prediction in this case! + +# distribution of covariates: We can improve the computation by assuming a +#distribution of the covariate (WT) in the individuals in the study. We set +#`groupsize=1`, the number of groups to be 50 (`m=50`) and assume that WT is +#sampled from a normal distribution with mean=70 and sd=10 +#(`a=as.list(rnorm(50,mean = 70,sd=10)`). +babel.db.2 <- nlmixr2(f, e, "poped", + control=popedControl( + groupsize=1, + m=50, + minxt=0, + maxxt=24, + bUseGrouped_xt = T, + a = lapply(rnorm(50, mean = 70, sd = 10), function(x) c(WT=x)) )) + + +evaluate_design(babel.db.2) + +#Here we see that, given this distribution of weights, the covariate effect +#parameters (bpop[3] and bpop[4]) would be well estimated. + +#However, we are only looking at one sample of 50 individuals. Maybe a better +#approach is to look at the distribution of RSEs over a number of experiments +#given the expected weight distribution. +nsim <- 10 +rse_list <- c() +for(i in 1:nsim){ + poped_db_tmp <- create.poped.database(ff_fun=mod_1, + fg_fun=par_1, + fError_fun=feps.add.prop, + groupsize=1, + m=50, + sigma=c(0.015,0.0015), + notfixed_sigma = c(1,0), + bpop=c(CL=3.8,V=20,WT_CL=0.75,WT_V=1), + d=c(CL=0.05,V=0.05), + xt=c( 1,2,4,6,8,24), + minxt=0, + maxxt=24, + bUseGrouped_xt=1, + a=lapply(rnorm(50, mean = 70, sd = 10), function(x) c(WT=x))) + rse_tmp <- evaluate_design(poped_db_tmp)$rse + rse_list <- rbind(rse_list,rse_tmp) +} +apply(rse_list,2,quantile) + From 068a1a215bccb6bb0b4d91c774c41d37da555ceb Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Wed, 16 Oct 2024 20:42:43 +0200 Subject: [PATCH 05/24] minor correction --- inst/poped/ex.12.covariate.distributions.babelmixr2.R | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/inst/poped/ex.12.covariate.distributions.babelmixr2.R b/inst/poped/ex.12.covariate.distributions.babelmixr2.R index b2d68728..f5043d2f 100644 --- a/inst/poped/ex.12.covariate.distributions.babelmixr2.R +++ b/inst/poped/ex.12.covariate.distributions.babelmixr2.R @@ -1,6 +1,5 @@ -rm(list=ls()) library(PopED) -library() +library(babelmixr2) ## Introduction # Lets assume that we have a model with a covariate included in the model description. # Here we define a one-compartment PK model that has weight on both clearance and volume of distribution. From 035555d74404e53967c3bdd1fba3b751e8b3d68a Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Wed, 16 Oct 2024 20:43:01 +0200 Subject: [PATCH 06/24] adding poped example 13 --- inst/poped/ex.13.shrinkage.babelmixr2.R | 48 +++++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 inst/poped/ex.13.shrinkage.babelmixr2.R diff --git a/inst/poped/ex.13.shrinkage.babelmixr2.R b/inst/poped/ex.13.shrinkage.babelmixr2.R new file mode 100644 index 00000000..f8088e89 --- /dev/null +++ b/inst/poped/ex.13.shrinkage.babelmixr2.R @@ -0,0 +1,48 @@ +library(PopED) +library(babelmixr2) + +## following example 1.4.2 in PFIM +f <- function() { + ini({ + tKA <- 2 + tK <- 0.25 + tV <- 15 + + eta.ka ~ 1 + eta.k ~ 0.1 + eta.v ~ 0.25 + + # Residual unexplained variability (RUV) + add.sd <- fix(0.5) # Additive error + prop.sd <- fix(0.15) # Proportional error + + }) + model({ + + KA <- tKA * exp(eta.ka) + K <- tK * exp(eta.k) + V <- tV * exp(eta.v) + + y = (DOSE/V*KA/(KA-K)*(exp(-K*time)-exp(-KA*time))) + y ~ prop(prop.sd) + add(add.sd) + + }) +} + +e <- et(c(1,3,8)) %>% + as.data.frame() + +## -- Define initial design and design space +babel.db <- nlmixr2(f, e, "poped", + control=popedControl( + groupsize=1, + minxt=0, + maxxt=10, + a=100)) + + +(shr_out <- shrinkage(babel.db)) +evaluate_design(babel.db) +plot_model_prediction(babel.db) + + From 2c314a8c04a401c52cdaf54a8d3d8af3d1b79af5 Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Wed, 16 Oct 2024 21:03:03 +0200 Subject: [PATCH 07/24] add poped example 4 --- inst/poped/ex.4.PKPD.1.comp.emax.babelmixr2.R | 98 +++++++++++++++++++ 1 file changed, 98 insertions(+) create mode 100644 inst/poped/ex.4.PKPD.1.comp.emax.babelmixr2.R diff --git a/inst/poped/ex.4.PKPD.1.comp.emax.babelmixr2.R b/inst/poped/ex.4.PKPD.1.comp.emax.babelmixr2.R new file mode 100644 index 00000000..fbd2a0c4 --- /dev/null +++ b/inst/poped/ex.4.PKPD.1.comp.emax.babelmixr2.R @@ -0,0 +1,98 @@ +library(PopED) +library(babelmixr2) +library(tidyverse) + +# Define the model +f <- function() { + + ini({ + # Initial estimates for population parameters (fixed effects) + tCL <- 0.5 # Clearance + tV <- 0.2 # Volume of distribution + tE0 <- 1 # Baseline effect + tEMAX <- 1 # Maximum effect + tEC50 <- 1 # Concentration at 50% of maximum effect + + # Variability (random effects) + eta.CL ~ 0.09 + eta.V ~ 0.09 + eta.E0 ~ 0.04 + eta.EC50 ~ 0.09 + + # Residual unexplained variability (RUV) + eps.pk.prop <- sqrt(0.15) # PK proportional error for(SD/mean) + eps.pdadd <- sqrt(0.015) # PD additive error + + }) + + model({ + + # PK model + # Using a simple analytical funtion. Notice the use of the variable "DOSE". This will be defined later in the poped controler + V <- tV * exp(eta.V) + CL <- tCL * exp(eta.CL) + cp = DOSE / V * exp(-CL / V * time) + + cp ~ prop(eps.pk.prop) + + # PD model + E0 <- tE0 * exp(eta.E0) + EMAX <- tEMAX + EC50 <- tEC50 * exp(eta.EC50) + + effect = E0 + cp * EMAX / (EC50 + cp) + effect ~ add(eps.pdadd) + + }) +} + +# e <- et(c(1,3,8)) %>% +# as.data.frame() + +event.table.optim.pk <- + # amt is a placeholder here + et(amt = 1) %>% + # specify observation records with appropriate bounds - these are observation windows + et(replicate(4, c(0, 5), simplify = FALSE)) %>% + # make sure that the appropriate times are kept + mutate(time = rep(c(0, 0.33,0.66,0.9,5), 1)) %>% + mutate(dvid = ifelse(time == 0, "dose", "cp")) + +event.table.optim.pd <- + # specify observation records with appropriate bounds - these are observation windows + et(replicate(4, c(0, 5), simplify = FALSE)) %>% + # make sure that the appropriate times are kept + mutate(time = rep(c(0.1, 1, 2, 5), 1), + dvid = 'effect') + +e <- event.table.optim.pk %>% bind_rows(event.table.optim.pd) %>% arrange(id, time) +# event.table.optim2 <- event.table.optim %>% filter(id ==2) +str(e) + + +## -- Define initial design and design space +babel.db <- nlmixr2(f, e, "poped", + control=popedControl( + groupsize=20, + m = 3, + bUseGrouped_xt = T, + minxt=0, + maxxt=5, + ourzero = 0, + a=list(c(DOSE=0),c(DOSE=1),c(DOSE=2)), + maxa = c(DOSE=10), + mina = c(DOSE=0)) ) + +plot_model_prediction(babel.db,facet_scales="free") +plot_model_prediction(babel.db,IPRED=T,DV=T,facet_scales="free",separate.groups=T) + +## evaluate initial design +evaluate_design(babel.db) +shrinkage(babel.db) + +# Optimization of sample times and doses +output <- poped_optim(babel.db, opt_xt = T, opt_a = T, parallel = T,method = c("LS")) + +get_rse(output$FIM,output$babel.db) +plot_model_prediction(output$babel.db,facet_scales="free") + From b9950a6e1411f0bb0594edcebcb3a625c41812e0 Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Wed, 16 Oct 2024 21:29:45 +0200 Subject: [PATCH 08/24] add poped example 11 --- inst/poped/ex.11.PK.prior.babelmixr2.R | 219 +++++++++++++++++++++++++ 1 file changed, 219 insertions(+) create mode 100644 inst/poped/ex.11.PK.prior.babelmixr2.R diff --git a/inst/poped/ex.11.PK.prior.babelmixr2.R b/inst/poped/ex.11.PK.prior.babelmixr2.R new file mode 100644 index 00000000..f302853d --- /dev/null +++ b/inst/poped/ex.11.PK.prior.babelmixr2.R @@ -0,0 +1,219 @@ +library(babelmixr2) + +library(PopED) +# This example shows how to include a prior FIM into the design evaluation. +# We look at PK assessment in pediatrics, where we are mainly interested +# in assessing if there is a substantial difference of more than 20% in +# clearance (CL) between children and adults. + +# First we setup the general model, then the designs for adults and pediatrics, +# and finally we can evaluate the separate designs and the pooled data. + +##-- Model: One comp first order absorption +## -- Analytic solution for both mutiple and single dosing + +# Define the model +f <- function() { + ini({ + tV <- 72.8 + tKa <- 0.25 + tCl <- 3.75 + tF <- fix(0.9) + pedCL <- 0.8 + + 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) * (pedCL**isPediatric) # add covariate for pediatrics + CL<-tCl*exp(eta.cl) + Favail <- tF + + N <- floor(time/TAU)+1 + y <- (DOSE*Favail/V)*(KA/(KA - CL/V)) * + (exp(-CL/V * (time - (N - 1) * TAU)) * + (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - + exp(-KA * (time - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU))) + + y ~ prop(prop.sd) + add(add.sd) + }) +} + +e <- et(c( 1,8,10,240,245)) + +babel.db <- nlmixr2(f, e, "poped", + popedControl(m = 2, + groupsize=20, + bUseGrouped_xt=TRUE, + a=list(c(DOSE=20,TAU=24,isPediatric = 0), + c(DOSE=40, TAU=24,isPediatric = 0)))) + + +# Note, to be able to use the adults FIM to combine with the pediatrics, +# both have to have the parameter "pedCL" defined and set notfixed_bpop to 1. + +## Define pediatric model/design (isPediatric = 1) +## One arm, 4 time points only + +e.ped <- et(c( 1,2,6,240)) + +babel.db.ped <- nlmixr2(f, e.ped, "poped", + popedControl(m = 1, + groupsize=6, + bUseGrouped_xt=TRUE, + a=list(c(DOSE=40,TAU=24,isPediatric = 1)))) + + +## Create plot of model of adult data without variability +plot_model_prediction(babel.db, model_num_points = 300) + +## Create plot of model of pediatric (single dose level) data without variability +plot_model_prediction(babel.db.ped, model_num_points = 300) + +## To store FIM from adult design - need FIM from evaluate_design +(outAdult <- evaluate_design(babel.db)) +# It is obvious that we cannot estimate the pediatric covariate from adult +# data only - therefore the message from the calculation. +# You can also note the zeros in the 4th column and 4th row of the FIM. + +# We can evaluate the adult design without warning, by setting the pedCL +# parameter to be fixed (i.e., not estimated) +evaluate_design(create.poped.database(babel.db, notfixed_bpop=c(1,1,1,0,0))) +# One obtains good estimates for all parameters for adults (<60% RSE for all). + +## evaluate design of pediatrics only - insufficient +# Similarly as before with only pediatrics we cannot estimate the covariate effect, +# so we fix it. +evaluate_design(create.poped.database(babel.db.ped, notfixed_bpop=c(1,1,1,0,0))) +# Due to having less subjects, less samples per subject, and only one dose level +# the variability in pediatrics cannot be estimated well. + +## Add adult prior +# Now we combined the two sutdies, where we assume that we only assess a difference +# between adults and pediatrics in CL. +# We can set the prior FIM to the adult one: +poped.db.all <- create.poped.database( + babel.db.ped, + prior_fim = outAdult$fim +) + +## evaluate design using prior FIM from adults +(out.all <- evaluate_design(babel.db.all)) +# Obviously, the pooled data leads to much higher precision in parameter estimates +# compared to the pediatrics only. + +# One can also obtain the power for estimating the covariate to be different from 1. +evaluate_power(babel.db.all, bpop_idx=5, h0=1,out=out.all) + + + + + + + + + + + + +# babel.db.ped <- create.poped.database(ff_fun="ff", +# fg_fun="sfg", +# fError_fun="feps", +# bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9,pedCL=0.8), +# notfixed_bpop=c(1,1,1,0,1), +# d=c(V=0.09,KA=0.09,CL=0.25^2), +# sigma=c(0.04,5e-6), +# notfixed_sigma=c(0,0), +# # m=1, +# # groupsize=6, +# # xt=c( 1,2,6,240), +# bUseGrouped_xt=1, +# a=list(c(DOSE=40,TAU=24,isPediatric = 1))) + + + + + +## -- Define design and design space for adults (isPediatric = 0) +## Two arms, 5 time points +poped.db <- create.poped.database(ff_fun="ff", + fg_fun="sfg", + fError_fun="feps", + # bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9,pedCL=0.8), + # notfixed_bpop=c(1,1,1,0,1), + # d=c(V=0.09,KA=0.09,CL=0.25^2), + # sigma=c(0.04,5e-6), + # notfixed_sigma=c(0,0), + # m=2, + # groupsize=20, + xt=c( 1,8,10,240,245), + bUseGrouped_xt=1, + a=list(c(DOSE=20,TAU=24,isPediatric = 0), + c(DOSE=40, TAU=24,isPediatric = 0))) + +# Note, to be able to use the adults FIM to combine with the pediatrics, +# both have to have the parameter "pedCL" defined and set notfixed_bpop to 1. + +## Define pediatric model/design (isPediatric = 1) +## One arm, 4 time points only +poped.db.ped <- create.poped.database(ff_fun="ff", + fg_fun="sfg", + fError_fun="feps", + bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9,pedCL=0.8), + notfixed_bpop=c(1,1,1,0,1), + d=c(V=0.09,KA=0.09,CL=0.25^2), + sigma=c(0.04,5e-6), + notfixed_sigma=c(0,0), + m=1, + groupsize=6, + xt=c( 1,2,6,240), + bUseGrouped_xt=1, + a=list(c(DOSE=40,TAU=24,isPediatric = 1))) + +## Create plot of model of adult data without variability +plot_model_prediction(poped.db, model_num_points = 300) + +## Create plot of model of pediatric (single dose level) data without variability +plot_model_prediction(poped.db.ped, model_num_points = 300) + +## To store FIM from adult design - need FIM from evaluate_design +(outAdult <- evaluate_design(poped.db)) +# It is obvious that we cannot estimate the pediatric covariate from adult +# data only - therefore the message from the calculation. +# You can also note the zeros in the 4th column and 4th row of the FIM. + +# We can evaluate the adult design without warning, by setting the pedCL +# parameter to be fixed (i.e., not estimated) +evaluate_design(create.poped.database(poped.db, notfixed_bpop=c(1,1,1,0,0))) +# One obtains good estimates for all parameters for adults (<60% RSE for all). + +## evaluate design of pediatrics only - insufficient +# Similarly as before with only pediatrics we cannot estimate the covariate effect, +# so we fix it. +evaluate_design(create.poped.database(poped.db.ped, notfixed_bpop=c(1,1,1,0,0))) +# Due to having less subjects, less samples per subject, and only one dose level +# the variability in pediatrics cannot be estimated well. + +## Add adult prior +# Now we combined the two sutdies, where we assume that we only assess a difference +# between adults and pediatrics in CL. +# We can set the prior FIM to the adult one: +poped.db.all <- create.poped.database( + poped.db.ped, + prior_fim = outAdult$fim +) + +## evaluate design using prior FIM from adults +(out.all <- evaluate_design(poped.db.all)) +# Obviously, the pooled data leads to much higher precision in parameter estimates +# compared to the pediatrics only. + +# One can also obtain the power for estimating the covariate to be different from 1. +evaluate_power(poped.db.all, bpop_idx=5, h0=1,out=out.all) + From 18ae0c433c10ea39bb242d4af7b52974ede10677 Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Wed, 16 Oct 2024 21:45:13 +0200 Subject: [PATCH 09/24] poped example 11 corrections. work in progress --- inst/poped/ex.11.PK.prior.babelmixr2.R | 130 ++++--------------------- 1 file changed, 19 insertions(+), 111 deletions(-) diff --git a/inst/poped/ex.11.PK.prior.babelmixr2.R b/inst/poped/ex.11.PK.prior.babelmixr2.R index f302853d..8cee506f 100644 --- a/inst/poped/ex.11.PK.prior.babelmixr2.R +++ b/inst/poped/ex.11.PK.prior.babelmixr2.R @@ -35,11 +35,11 @@ f <- function() { CL<-tCl*exp(eta.cl) Favail <- tF - N <- floor(time/TAU)+1 + N <- floor(t/TAU)+1 y <- (DOSE*Favail/V)*(KA/(KA - CL/V)) * - (exp(-CL/V * (time - (N - 1) * TAU)) * + (exp(-CL/V * (t - (N - 1) * TAU)) * (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - - exp(-KA * (time - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU))) + exp(-KA * (t - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU))) y ~ prop(prop.sd) + add(add.sd) }) @@ -84,13 +84,27 @@ plot_model_prediction(babel.db.ped, model_num_points = 300) # We can evaluate the adult design without warning, by setting the pedCL # parameter to be fixed (i.e., not estimated) -evaluate_design(create.poped.database(babel.db, notfixed_bpop=c(1,1,1,0,0))) +f2 <- f %>% ini(pedCL = fix(0.8)) +babel.db.2 <- nlmixr2(f2, e, "poped", + popedControl(m = 2, + groupsize=20, + bUseGrouped_xt=TRUE, + a=list(c(DOSE=20,TAU=24,isPediatric = 0), + c(DOSE=40, TAU=24,isPediatric = 0)))) + +evaluate_design(babel.db.2) # One obtains good estimates for all parameters for adults (<60% RSE for all). ## evaluate design of pediatrics only - insufficient # Similarly as before with only pediatrics we cannot estimate the covariate effect, # so we fix it. -evaluate_design(create.poped.database(babel.db.ped, notfixed_bpop=c(1,1,1,0,0))) +babel.db.ped.2 <- nlmixr2(f2, e.ped, "poped", + popedControl(m = 1, + groupsize=6, + bUseGrouped_xt=TRUE, + a=list(c(DOSE=40,TAU=24,isPediatric = 1)))) + +evaluate_design(babel.db.ped.2) # Due to having less subjects, less samples per subject, and only one dose level # the variability in pediatrics cannot be estimated well. @@ -111,109 +125,3 @@ poped.db.all <- create.poped.database( # One can also obtain the power for estimating the covariate to be different from 1. evaluate_power(babel.db.all, bpop_idx=5, h0=1,out=out.all) - - - - - - - - - - - -# babel.db.ped <- create.poped.database(ff_fun="ff", -# fg_fun="sfg", -# fError_fun="feps", -# bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9,pedCL=0.8), -# notfixed_bpop=c(1,1,1,0,1), -# d=c(V=0.09,KA=0.09,CL=0.25^2), -# sigma=c(0.04,5e-6), -# notfixed_sigma=c(0,0), -# # m=1, -# # groupsize=6, -# # xt=c( 1,2,6,240), -# bUseGrouped_xt=1, -# a=list(c(DOSE=40,TAU=24,isPediatric = 1))) - - - - - -## -- Define design and design space for adults (isPediatric = 0) -## Two arms, 5 time points -poped.db <- create.poped.database(ff_fun="ff", - fg_fun="sfg", - fError_fun="feps", - # bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9,pedCL=0.8), - # notfixed_bpop=c(1,1,1,0,1), - # d=c(V=0.09,KA=0.09,CL=0.25^2), - # sigma=c(0.04,5e-6), - # notfixed_sigma=c(0,0), - # m=2, - # groupsize=20, - xt=c( 1,8,10,240,245), - bUseGrouped_xt=1, - a=list(c(DOSE=20,TAU=24,isPediatric = 0), - c(DOSE=40, TAU=24,isPediatric = 0))) - -# Note, to be able to use the adults FIM to combine with the pediatrics, -# both have to have the parameter "pedCL" defined and set notfixed_bpop to 1. - -## Define pediatric model/design (isPediatric = 1) -## One arm, 4 time points only -poped.db.ped <- create.poped.database(ff_fun="ff", - fg_fun="sfg", - fError_fun="feps", - bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9,pedCL=0.8), - notfixed_bpop=c(1,1,1,0,1), - d=c(V=0.09,KA=0.09,CL=0.25^2), - sigma=c(0.04,5e-6), - notfixed_sigma=c(0,0), - m=1, - groupsize=6, - xt=c( 1,2,6,240), - bUseGrouped_xt=1, - a=list(c(DOSE=40,TAU=24,isPediatric = 1))) - -## Create plot of model of adult data without variability -plot_model_prediction(poped.db, model_num_points = 300) - -## Create plot of model of pediatric (single dose level) data without variability -plot_model_prediction(poped.db.ped, model_num_points = 300) - -## To store FIM from adult design - need FIM from evaluate_design -(outAdult <- evaluate_design(poped.db)) -# It is obvious that we cannot estimate the pediatric covariate from adult -# data only - therefore the message from the calculation. -# You can also note the zeros in the 4th column and 4th row of the FIM. - -# We can evaluate the adult design without warning, by setting the pedCL -# parameter to be fixed (i.e., not estimated) -evaluate_design(create.poped.database(poped.db, notfixed_bpop=c(1,1,1,0,0))) -# One obtains good estimates for all parameters for adults (<60% RSE for all). - -## evaluate design of pediatrics only - insufficient -# Similarly as before with only pediatrics we cannot estimate the covariate effect, -# so we fix it. -evaluate_design(create.poped.database(poped.db.ped, notfixed_bpop=c(1,1,1,0,0))) -# Due to having less subjects, less samples per subject, and only one dose level -# the variability in pediatrics cannot be estimated well. - -## Add adult prior -# Now we combined the two sutdies, where we assume that we only assess a difference -# between adults and pediatrics in CL. -# We can set the prior FIM to the adult one: -poped.db.all <- create.poped.database( - poped.db.ped, - prior_fim = outAdult$fim -) - -## evaluate design using prior FIM from adults -(out.all <- evaluate_design(poped.db.all)) -# Obviously, the pooled data leads to much higher precision in parameter estimates -# compared to the pediatrics only. - -# One can also obtain the power for estimating the covariate to be different from 1. -evaluate_power(poped.db.all, bpop_idx=5, h0=1,out=out.all) - From baa796bbc2d8d918d395e575761181f61033a7c0 Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Thu, 17 Oct 2024 19:57:54 +0200 Subject: [PATCH 10/24] poped example 10 added - not working right now --- inst/poped/ex.10.PKPD.HCV.babelmixr2.R | 152 +++++++++++++++++++++++++ 1 file changed, 152 insertions(+) create mode 100644 inst/poped/ex.10.PKPD.HCV.babelmixr2.R diff --git a/inst/poped/ex.10.PKPD.HCV.babelmixr2.R b/inst/poped/ex.10.PKPD.HCV.babelmixr2.R new file mode 100644 index 00000000..095dc77a --- /dev/null +++ b/inst/poped/ex.10.PKPD.HCV.babelmixr2.R @@ -0,0 +1,152 @@ +## HCV example from software comparison in: +## Nyberg et al., "Methods and software tools for design evaluation +## for population pharmacokinetics-pharmacodynamics studies", +## Br. J. Clin. Pharm., 2014. + +library(babelmixr2) +library(PopED) + +f <- function() { + ini({ + + tp <- fix(100) + td <- fix(0.001) + te <- fix(1e-7) + ts <- fix(20000) + + tKA <- log(0.8) + tKE <- log(0.15) + tVD <- log(100) #VD=log(100000), + tEC50 <- log(0.12) #EC50=log(0.00012), + tn <- log(2) + tdelta <- log(0.2) + tc <- log(7) + + eta.KA ~ 0.25 + eta.KE ~ 0.25 + eta.VD ~ 0.25 + eta.EC50 ~0.25 + eta.n ~ 0.25 + eta.delta ~ 0.25 + eta.c ~ 0.25 + + add.sd.pk <- sqrt(0.04) # nlmixr2 uses sd + add.sd.pd <- sqrt(0.04) + }) + model({ + + p <- tp + d <- td + e <- te + s <- ts + KA <- exp(tKA + eta.KA) + KE <- exp(tKE + eta.KE) + VD <- exp(tVD + eta.VD) + EC50 <- exp(tEC50 + eta.EC50) + n <- exp(tn + eta.n) + delta <- exp(tdelta + eta.delta) + c <- exp(tc + eta.c) + + # # Apply the conditional logic + if ( (time - (floor(time / TAU) * TAU)) < TINF) { + In <- DOSE / TINF # Infusion rate + } else { + In <- 0 # No infusion + } + + d/dt(depot) = -KA*depot + In + d/dt(central) = KA*depot - KE*central + d/dt(TC) = s - TC*(e*VP+d) # target cells (TC) + d/dt(IC) = e*TC*VP-delta*IC # productively infected cells (IC) + d/dt(VP) = p*(1-(pow(central/VD,n)/(pow(central/VD,n)+pow(EC50,n))))*IC-c*VP # viral particles (VP) + + # Initial conditions + TC(0) =c*delta/(p*e) + IC(0) =(s*e*p-d*c*delta)/(p*delta*e) + VP(0) = (s*e*p-d*c*delta)/(c*delta*e) + + # Concentration and effect are calculated + conc = central/VD + eff = log10(VP) + ## And both are assumed to follow additive error + conc ~ add(add.sd.pk) + eff ~ add(add.sd.pd) + }) +} + +rxClean() +# time%%TAU +# time = seq(0,30,0.1) + +f <- f() + +# Note that design point 240 is repeated +e1 <- et(c( 0,0.25,0.5,1,2,3,4,7,10,14,21,28)) %>% + as.data.frame() %>% + dplyr::mutate(dvid=1) + +e2 <- e1 %>% dplyr::mutate(dvid=2) +e <- rbind(e1, e2) + +babel.db <- nlmixr2(f, e, "poped", + popedControl( + groupsize=30, + a=list(c(DOSE=180,TINF=1,TAU=7)) + )) + +## create plot of model without variability +plot_model_prediction(babel.db, facet_scales = "free") +evaluate_design(babel.db) +shrinkage(babel.db) + +#' ##################################### +#' The reduced FIM +#' #################################### + # computation time +tic(); FIM_compiled <- evaluate.fim(babel.db); toc() + +#' design evaluation +crit <- det(FIM_compiled)^(1/length(get_unfixed_params(babel.db)[["all"]])) +crit +rse <- get_rse(FIM_compiled,babel.db) # this is for the log of the fixed effect parameters +rse_norm <- sqrt(diag(inv(FIM_compiled)))*100 # this is approximately the RSE for the normal scale of the fixed effects +rse[1:7] <- rse_norm[1:7] # replace the log scale for the normal scale RSE values +rse + +#' Evaluation comparison to +#' Nyberg et al., "Methods and software tools for design evaluation +#' for population pharmacokinetics-pharmacodynamics studies", +#' Br. J. Clin. Pharm., 2014. +crit_reference_reduced <- 248.8 +rse_reference_reduced <- c(12.1,10.5,10.0,15.8,10.4,9.4,11.0,40.0,30.8,28.8,60.4,28.8,27.2,32.8,8.5,9.0) + +# the relative differences in percent +(rse - rse_reference_reduced)/rse_reference_reduced * 100 +(crit - crit_reference_reduced)/crit_reference_reduced * 100 + +#' ##################################### +#' The full FIM. +#' ##################################### +FIM_compiled_full <- evaluate.fim(poped_db_compiled,fim.calc.type = 0) + +#' design evaluation +crit_full <- det(FIM_compiled_full)^(1/length(get_unfixed_params(poped_db_compiled)[["all"]])) +crit_full +rse_full <- get_rse(FIM_compiled_full,poped_db_compiled) # this is for the log of the fixed effect parameters +rse_norm_full <- sqrt(diag(inv(FIM_compiled_full)))*100 # this is approximately the RSE for the normal scale of the fixed effects +rse_full[1:7] <- rse_norm_full[1:7] # replace the log scale for the normal scale RSE values +rse_full + +#' Evaluation compared to +#' Nyberg et al., "Methods and software tools for design evaluation +#' for population pharmacokinetics-pharmacodynamics studies", +#' Br. J. Clin. Pharm., 2014. +crit_reference_full <- 318.2 +rse_reference_full <- c(8.6,6.9,8.4,13.5,7.5,8.5,8.7,43.2,37.2,33.2,66.4,32.8,31.6,33.6,8.5,9.3) + +# the relative differences in percent +(rse_full - rse_reference_full)/rse_reference_full * 100 +(crit_full - crit_reference_full)/crit_reference_full * 100 + + + From 602902ec4b42dbd8ea4856d01082cea10c29b2fb Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Thu, 17 Oct 2024 21:49:55 +0200 Subject: [PATCH 11/24] IOV example addtion --- inst/poped/ex.14.PK.IOV.babelmixr2.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/inst/poped/ex.14.PK.IOV.babelmixr2.R b/inst/poped/ex.14.PK.IOV.babelmixr2.R index 64cd5445..abf83cf2 100644 --- a/inst/poped/ex.14.PK.IOV.babelmixr2.R +++ b/inst/poped/ex.14.PK.IOV.babelmixr2.R @@ -2,6 +2,8 @@ library(babelmixr2) library(PopED) library(ggplot2) +# Does not work for now. PopED has a special way to handle IOV. This feature will be added in the future + # # N = floor(Time/TAU)+1; # CL = CL_OCC_1; From 44ff4aa4428c92c9a68a3866f9b8244a05615bbb Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Fri, 18 Oct 2024 11:08:47 +0200 Subject: [PATCH 12/24] updates on ex.10 - still under developemtn --- inst/poped/ex.10.PKPD.HCV.babelmixr2.R | 38 ++++++++++++++++---------- 1 file changed, 23 insertions(+), 15 deletions(-) diff --git a/inst/poped/ex.10.PKPD.HCV.babelmixr2.R b/inst/poped/ex.10.PKPD.HCV.babelmixr2.R index 095dc77a..c1bd5ca4 100644 --- a/inst/poped/ex.10.PKPD.HCV.babelmixr2.R +++ b/inst/poped/ex.10.PKPD.HCV.babelmixr2.R @@ -47,20 +47,24 @@ f <- function() { delta <- exp(tdelta + eta.delta) c <- exp(tc + eta.c) - # # Apply the conditional logic - if ( (time - (floor(time / TAU) * TAU)) < TINF) { - In <- DOSE / TINF # Infusion rate - } else { - In <- 0 # No infusion - } - - d/dt(depot) = -KA*depot + In + # # Apply the conditional logic - + # # This does not work right now - could be addressed in the future + # if ( (time - (floor(time / TAU) * TAU)) < TINF) { + # In <- DOSE / TINF # Infusion rate + # } else { + # In <- 0 # No infusion + # } + + d/dt(depot) = -KA*depot # + In d/dt(central) = KA*depot - KE*central d/dt(TC) = s - TC*(e*VP+d) # target cells (TC) d/dt(IC) = e*TC*VP-delta*IC # productively infected cells (IC) d/dt(VP) = p*(1-(pow(central/VD,n)/(pow(central/VD,n)+pow(EC50,n))))*IC-c*VP # viral particles (VP) # Initial conditions + dur(depot) = TINF + # depot(0) <- DOSE + TC(0) =c*delta/(p*e) IC(0) =(s*e*p-d*c*delta)/(p*delta*e) VP(0) = (s*e*p-d*c*delta)/(c*delta*e) @@ -75,24 +79,28 @@ f <- function() { } rxClean() -# time%%TAU -# time = seq(0,30,0.1) - f <- f() -# Note that design point 240 is repeated +TAU = 7 e1 <- et(c( 0,0.25,0.5,1,2,3,4,7,10,14,21,28)) %>% + add.dosing(dose=180, nbr.doses=4, dosing.interval=TAU) %>% as.data.frame() %>% dplyr::mutate(dvid=1) -e2 <- e1 %>% dplyr::mutate(dvid=2) +e2 <- et(c( 0,0.25,0.5,1,2,3,4,7,10,14,21,28)) %>% + add.dosing(dose=0, nbr.doses=4, dosing.interval=TAU) %>% + as.data.frame() %>% + dplyr::mutate(dvid=2) %>% + dplyr::filter(is.na(ii)) + e <- rbind(e1, e2) babel.db <- nlmixr2(f, e, "poped", popedControl( groupsize=30, - a=list(c(DOSE=180,TINF=1,TAU=7)) - )) + a=list(c(TINF=1)))) + # a=list(c(DOSE=180,TINF=1,TAU=7)) + # )) ## create plot of model without variability plot_model_prediction(babel.db, facet_scales = "free") From 433f82dae616a3eeeef6de34acb280192845167d Mon Sep 17 00:00:00 2001 From: Matthew Fidler Date: Fri, 18 Oct 2024 18:39:53 -0500 Subject: [PATCH 13/24] Fix example 11 --- inst/poped/ex.11.PK.prior.babelmixr2.R | 145 ++++++++++++++++--------- 1 file changed, 94 insertions(+), 51 deletions(-) diff --git a/inst/poped/ex.11.PK.prior.babelmixr2.R b/inst/poped/ex.11.PK.prior.babelmixr2.R index 8cee506f..04362508 100644 --- a/inst/poped/ex.11.PK.prior.babelmixr2.R +++ b/inst/poped/ex.11.PK.prior.babelmixr2.R @@ -2,11 +2,11 @@ library(babelmixr2) library(PopED) # This example shows how to include a prior FIM into the design evaluation. -# We look at PK assessment in pediatrics, where we are mainly interested -# in assessing if there is a substantial difference of more than 20% in +# We look at PK assessment in pediatrics, where we are mainly interested +# in assessing if there is a substantial difference of more than 20% in # clearance (CL) between children and adults. -# First we setup the general model, then the designs for adults and pediatrics, +# First we setup the general model, then the designs for adults and pediatrics, # and finally we can evaluate the separate designs and the pooled data. ##-- Model: One comp first order absorption @@ -14,35 +14,35 @@ library(PopED) # Define the model f <- function() { - ini({ - tV <- 72.8 - tKa <- 0.25 - tCl <- 3.75 - tF <- fix(0.9) - pedCL <- 0.8 - - 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) * (pedCL**isPediatric) # add covariate for pediatrics - CL<-tCl*exp(eta.cl) - Favail <- tF - - N <- floor(t/TAU)+1 - y <- (DOSE*Favail/V)*(KA/(KA - CL/V)) * - (exp(-CL/V * (t - (N - 1) * TAU)) * - (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - - exp(-KA * (t - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU))) - - y ~ prop(prop.sd) + add(add.sd) - }) + ini({ + tV <- 72.8 + tKa <- 0.25 + tCl <- 3.75 + tF <- fix(0.9) + pedCL <- 0.8 + + 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)* (pedCL^(isPediatric)) # add covariate for pediatrics + Favail <- tF + + N <- floor(t/TAU)+1 + y <- (DOSE*Favail/V)*(KA/(KA - CL/V)) * + (exp(-CL/V * (t - (N - 1) * TAU)) * + (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - + exp(-KA * (t - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU))) + + y ~ prop(prop.sd) + add(add.sd) + }) } e <- et(c( 1,8,10,240,245)) @@ -55,7 +55,7 @@ babel.db <- nlmixr2(f, e, "poped", c(DOSE=40, TAU=24,isPediatric = 0)))) -# Note, to be able to use the adults FIM to combine with the pediatrics, +# Note, to be able to use the adults FIM to combine with the pediatrics, # both have to have the parameter "pedCL" defined and set notfixed_bpop to 1. ## Define pediatric model/design (isPediatric = 1) @@ -65,10 +65,10 @@ e.ped <- et(c( 1,2,6,240)) babel.db.ped <- nlmixr2(f, e.ped, "poped", popedControl(m = 1, - groupsize=6, - bUseGrouped_xt=TRUE, - a=list(c(DOSE=40,TAU=24,isPediatric = 1)))) - + groupsize=6, + bUseGrouped_xt=TRUE, + a=list(c(DOSE=40,TAU=24,isPediatric = 1)))) + ## Create plot of model of adult data without variability plot_model_prediction(babel.db, model_num_points = 300) @@ -77,14 +77,25 @@ plot_model_prediction(babel.db, model_num_points = 300) plot_model_prediction(babel.db.ped, model_num_points = 300) ## To store FIM from adult design - need FIM from evaluate_design + +## $rse +## V KA CL pedCL d_V d_KA d_CL +## 6.634931 8.587203 4.354792 NA 33.243601 55.689432 27.133255 + +## Warning message: +## The following parameters are not estimable: +## pedCL +## Is the design adequate to estimate all parameters? + (outAdult <- evaluate_design(babel.db)) -# It is obvious that we cannot estimate the pediatric covariate from adult +# It is obvious that we cannot estimate the pediatric covariate from adult # data only - therefore the message from the calculation. # You can also note the zeros in the 4th column and 4th row of the FIM. -# We can evaluate the adult design without warning, by setting the pedCL -# parameter to be fixed (i.e., not estimated) +# We can evaluate the adult design without warning, by setting the pedCLs +# parameter to be fixed (i.e., not estimated)c f2 <- f %>% ini(pedCL = fix(0.8)) + babel.db.2 <- nlmixr2(f2, e, "poped", popedControl(m = 2, groupsize=20, @@ -92,36 +103,68 @@ babel.db.2 <- nlmixr2(f2, e, "poped", a=list(c(DOSE=20,TAU=24,isPediatric = 0), c(DOSE=40, TAU=24,isPediatric = 0)))) +## PopED +## +## $rse +## V KA CL d_V d_KA d_CL +## 6.634931 8.587203 4.354792 33.243601 55.689432 27.133255 evaluate_design(babel.db.2) # One obtains good estimates for all parameters for adults (<60% RSE for all). ## evaluate design of pediatrics only - insufficient -# Similarly as before with only pediatrics we cannot estimate the covariate effect, +# Similarly as before with only pediatrics we cannot estimate the covariate effect, # so we fix it. babel.db.ped.2 <- nlmixr2(f2, e.ped, "poped", popedControl(m = 1, groupsize=6, bUseGrouped_xt=TRUE, - a=list(c(DOSE=40,TAU=24,isPediatric = 1)))) + a=list(c(DOSE=40,TAU=24,isPediatric = 1)), + literalFix=FALSE)) + + +## PopED +## $rse +## V KA CL d_V d_KA d_CL +## 24.72088 30.84953 11.94767 116.23095 181.19778 77.29188 evaluate_design(babel.db.ped.2) -# Due to having less subjects, less samples per subject, and only one dose level +# Due to having less subjects, less samples per subject, and only one dose level # the variability in pediatrics cannot be estimated well. ## Add adult prior -# Now we combined the two sutdies, where we assume that we only assess a difference -# between adults and pediatrics in CL. -# We can set the prior FIM to the adult one: -poped.db.all <- create.poped.database( - babel.db.ped, - prior_fim = outAdult$fim +# Now we combined the two sutdies, where we assume that we only assess +# a difference between adults and pediatrics in CL. We can set the +# prior FIM to the adult one: Note to setup babelmixr2 models +# appropriately we need to use `babel.poped.database` instead of +# `create.poped.database` +babel.db.all <- babel.poped.database( + babel.db.ped, + prior_fim = outAdult$fim ) + ## evaluate design using prior FIM from adults + +## PopED +## $rse +## V KA CL pedCL d_V d_KA d_CL +## 6.381388 8.222819 4.354761 12.591940 31.808871 52.858399 25.601551 + (out.all <- evaluate_design(babel.db.all)) -# Obviously, the pooled data leads to much higher precision in parameter estimates +# Obviously, the pooled data leads to much higher precision in parameter estimates # compared to the pediatrics only. # One can also obtain the power for estimating the covariate to be different from 1. -evaluate_power(babel.db.all, bpop_idx=5, h0=1,out=out.all) +## PopED +## +## $rse +## V KA CL pedCL d_V d_KA d_CL +## 6.381388 8.222819 4.354761 12.591940 31.808871 52.858399 25.601551 + +## $power +## Value RSE power_pred power_want need_rse min_N_tot +## pedCL 0.8 12.59194 51.01851 80 8.923519 14 +evaluate_power(babel.db.all, + bpop_idx=babelBpopIdx(babel.db.all, "pedCL"), + h0=1,out=out.all) From b1e6ae46c7486c759eaae9a9b8f7a71574daf12b Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Tue, 22 Oct 2024 16:03:45 +0200 Subject: [PATCH 14/24] typo correction --- inst/poped/ex.11.PK.prior.babelmixr2.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inst/poped/ex.11.PK.prior.babelmixr2.R b/inst/poped/ex.11.PK.prior.babelmixr2.R index 04362508..f462e598 100644 --- a/inst/poped/ex.11.PK.prior.babelmixr2.R +++ b/inst/poped/ex.11.PK.prior.babelmixr2.R @@ -132,7 +132,7 @@ evaluate_design(babel.db.ped.2) # the variability in pediatrics cannot be estimated well. ## Add adult prior -# Now we combined the two sutdies, where we assume that we only assess +# Now we combined the two studies, where we assume that we only assess # a difference between adults and pediatrics in CL. We can set the # prior FIM to the adult one: Note to setup babelmixr2 models # appropriately we need to use `babel.poped.database` instead of From 47e1187b823afa289b6f942eff4f79be74fbcbc0 Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Tue, 22 Oct 2024 16:25:19 +0200 Subject: [PATCH 15/24] example 10 corrections added --- inst/poped/ex.10.PKPD.HCV.babelmixr2.R | 35 +++++++++++++------------- 1 file changed, 17 insertions(+), 18 deletions(-) diff --git a/inst/poped/ex.10.PKPD.HCV.babelmixr2.R b/inst/poped/ex.10.PKPD.HCV.babelmixr2.R index c1bd5ca4..90c0fd7a 100644 --- a/inst/poped/ex.10.PKPD.HCV.babelmixr2.R +++ b/inst/poped/ex.10.PKPD.HCV.babelmixr2.R @@ -46,7 +46,7 @@ f <- function() { n <- exp(tn + eta.n) delta <- exp(tdelta + eta.delta) c <- exp(tc + eta.c) - + # # Apply the conditional logic - # # This does not work right now - could be addressed in the future # if ( (time - (floor(time / TAU) * TAU)) < TINF) { @@ -54,7 +54,7 @@ f <- function() { # } else { # In <- 0 # No infusion # } - + d/dt(depot) = -KA*depot # + In d/dt(central) = KA*depot - KE*central d/dt(TC) = s - TC*(e*VP+d) # target cells (TC) @@ -62,7 +62,7 @@ f <- function() { d/dt(VP) = p*(1-(pow(central/VD,n)/(pow(central/VD,n)+pow(EC50,n))))*IC-c*VP # viral particles (VP) # Initial conditions - dur(depot) = TINF + dur(depot) = 1 # depot(0) <- DOSE TC(0) =c*delta/(p*e) @@ -83,24 +83,23 @@ f <- f() TAU = 7 e1 <- et(c( 0,0.25,0.5,1,2,3,4,7,10,14,21,28)) %>% - add.dosing(dose=180, nbr.doses=4, dosing.interval=TAU) %>% + add.dosing(dose=180, nbr.doses=4, dosing.interval=TAU) %>% as.data.frame() %>% dplyr::mutate(dvid=1) e2 <- et(c( 0,0.25,0.5,1,2,3,4,7,10,14,21,28)) %>% - add.dosing(dose=0, nbr.doses=4, dosing.interval=TAU) %>% + add.dosing(dose=0, nbr.doses=4, dosing.interval=TAU) %>% as.data.frame() %>% - dplyr::mutate(dvid=2) %>% + dplyr::mutate(dvid=2) %>% dplyr::filter(is.na(ii)) e <- rbind(e1, e2) babel.db <- nlmixr2(f, e, "poped", popedControl( - groupsize=30, - a=list(c(TINF=1)))) - # a=list(c(DOSE=180,TINF=1,TAU=7)) - # )) + groupsize=30 + # ,a=list(c(DOSE=180,TINF=1,TAU=7)) + )) ## create plot of model without variability plot_model_prediction(babel.db, facet_scales = "free") @@ -111,13 +110,13 @@ shrinkage(babel.db) #' The reduced FIM #' #################################### # computation time -tic(); FIM_compiled <- evaluate.fim(babel.db); toc() +tic(); FIM_babel <- evaluate.fim(babel.db); toc() #' design evaluation -crit <- det(FIM_compiled)^(1/length(get_unfixed_params(babel.db)[["all"]])) +crit <- det(FIM_babel)^(1/length(get_unfixed_params(babel.db)[["all"]])) crit -rse <- get_rse(FIM_compiled,babel.db) # this is for the log of the fixed effect parameters -rse_norm <- sqrt(diag(inv(FIM_compiled)))*100 # this is approximately the RSE for the normal scale of the fixed effects +rse <- get_rse(FIM_babel,babel.db) # this is for the log of the fixed effect parameters +rse_norm <- sqrt(diag(inv(FIM_babel)))*100 # this is approximately the RSE for the normal scale of the fixed effects rse[1:7] <- rse_norm[1:7] # replace the log scale for the normal scale RSE values rse @@ -135,13 +134,13 @@ rse_reference_reduced <- c(12.1,10.5,10.0,15.8,10.4,9.4,11.0,40.0,30.8,28.8,60.4 #' ##################################### #' The full FIM. #' ##################################### -FIM_compiled_full <- evaluate.fim(poped_db_compiled,fim.calc.type = 0) +FIM_babel_full <- evaluate.fim(babel.db,fim.calc.type = 0) #' design evaluation -crit_full <- det(FIM_compiled_full)^(1/length(get_unfixed_params(poped_db_compiled)[["all"]])) +crit_full <- det(FIM_babel_full)^(1/length(get_unfixed_params(babel.db)[["all"]])) crit_full -rse_full <- get_rse(FIM_compiled_full,poped_db_compiled) # this is for the log of the fixed effect parameters -rse_norm_full <- sqrt(diag(inv(FIM_compiled_full)))*100 # this is approximately the RSE for the normal scale of the fixed effects +rse_full <- get_rse(FIM_babel_full,babel.db) # this is for the log of the fixed effect parameters +rse_norm_full <- sqrt(diag(inv(FIM_babel_full)))*100 # this is approximately the RSE for the normal scale of the fixed effects rse_full[1:7] <- rse_norm_full[1:7] # replace the log scale for the normal scale RSE values rse_full From d43e61e3413444f44926a5e5e8e4025bfb8815c8 Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Tue, 22 Oct 2024 16:33:21 +0200 Subject: [PATCH 16/24] Example 10 - corrected infusion in et. --- inst/poped/ex.10.PKPD.HCV.babelmixr2.R | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/inst/poped/ex.10.PKPD.HCV.babelmixr2.R b/inst/poped/ex.10.PKPD.HCV.babelmixr2.R index 90c0fd7a..ace7eac7 100644 --- a/inst/poped/ex.10.PKPD.HCV.babelmixr2.R +++ b/inst/poped/ex.10.PKPD.HCV.babelmixr2.R @@ -82,13 +82,14 @@ rxClean() f <- f() TAU = 7 +TINF = 1 e1 <- et(c( 0,0.25,0.5,1,2,3,4,7,10,14,21,28)) %>% - add.dosing(dose=180, nbr.doses=4, dosing.interval=TAU) %>% + add.dosing(dose=180, nbr.doses=4, dosing.interval=TAU, rate = 180/TINF) %>% as.data.frame() %>% dplyr::mutate(dvid=1) e2 <- et(c( 0,0.25,0.5,1,2,3,4,7,10,14,21,28)) %>% - add.dosing(dose=0, nbr.doses=4, dosing.interval=TAU) %>% + add.dosing(dose=0, nbr.doses=4, dosing.interval=TAU, rate = 180/TINF) %>% as.data.frame() %>% dplyr::mutate(dvid=2) %>% dplyr::filter(is.na(ii)) @@ -102,7 +103,7 @@ babel.db <- nlmixr2(f, e, "poped", )) ## create plot of model without variability -plot_model_prediction(babel.db, facet_scales = "free") +plot_model_prediction(babel.db, facet_scales = "free", model_num_points = 1000) evaluate_design(babel.db) shrinkage(babel.db) From 6729283f747322ec1fbe5665da0ae5549838d5c0 Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Wed, 23 Oct 2024 12:23:03 +0200 Subject: [PATCH 17/24] laplace turned on, and parallelization turned off --- inst/poped/ex.2.d.warfarin.ED.babelmixr2.R | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/inst/poped/ex.2.d.warfarin.ED.babelmixr2.R b/inst/poped/ex.2.d.warfarin.ED.babelmixr2.R index ee0a13b8..582759f0 100644 --- a/inst/poped/ex.2.d.warfarin.ED.babelmixr2.R +++ b/inst/poped/ex.2.d.warfarin.ED.babelmixr2.R @@ -43,7 +43,6 @@ babel.db <- nlmixr2(f, e, "poped", maxa=100)) - # Adding 10% Uncertainty to all fixed effects (not Favail) bpop_vals_ed <- babel.db$parameters$bpop for (n in row.names(bpop_vals_ed)) { @@ -79,15 +78,14 @@ tic();evaluate_design(babel.db,d_switch=FALSE,ED_samp_size=20); toc() ## 4.991010 2.977982 14.014207 29.802546 36.711408 26.754059 31.477157 25.297312 ## optimization with line search search -output_ls <- poped_optim(babel.db, opt_xt=T, parallel=T, method = "LS", d_switch=F, ED_samp_size=20) +## parallelization does not seem to work with babelmixr2 or example: +output_ls <- poped_optim(babel.db, opt_xt=T, parallel=F, method = "LS", d_switch=F, ED_samp_size=20) ## laplace does not seem to work with babelmixr2 or example: -## See - ## ED: E(det(FIM)) using Laplace approximation ## deterministic calculation, relatively fast ## can be more stable for optimization -## tic(); evaluate_design(babel.db,d_switch=FALSE,use_laplace=TRUE); toc() +tic(); evaluate_design(babel.db,d_switch=FALSE,use_laplace=TRUE); toc() ## optimization with Laplace -## output_ls <- poped_optim(babel.db, opt_xt=T, parallel=T, method = "LS", d_switch=F, use_laplace=T) +output_ls <- poped_optim(babel.db, opt_xt=T, parallel=F, method = "LS", d_switch=F, use_laplace=T) From df1e2eea024e87c9ac480d4e7a19e7f152ff1566 Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Wed, 23 Oct 2024 12:23:35 +0200 Subject: [PATCH 18/24] minor corrections --- .../ex.1.a.PK.1.comp.oral.md.intro.babelmixr2.R | 14 ++++++++------ ....PK.1.comp.oral.md.re-parameterize.babelmixr2.R | 6 ++++++ 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/inst/poped/ex.1.a.PK.1.comp.oral.md.intro.babelmixr2.R b/inst/poped/ex.1.a.PK.1.comp.oral.md.intro.babelmixr2.R index ea8a5a87..9166729c 100644 --- a/inst/poped/ex.1.a.PK.1.comp.oral.md.intro.babelmixr2.R +++ b/inst/poped/ex.1.a.PK.1.comp.oral.md.intro.babelmixr2.R @@ -35,7 +35,6 @@ f <- function() { }) } -# minxt, maxxt e <- et(list(c(0, 10), c(0, 10), c(0, 10), @@ -69,7 +68,6 @@ plot_model_prediction(babel.db, IPRED=T, DV=T, separate.groups=T, model_num_poin evaluate_design(babel.db) - ## original: > shrinkage(poped.db) ## # A tibble: 9 × 5 ## d_V d_KA d_CL type group @@ -92,34 +90,38 @@ output <- poped_optim(babel.db, opt_xt =TRUE) # Evaluate optimization results summary(output) +## From original +# V KA CL d_V d_KA d_CL +# 6.281944 7.726279 4.295908 32.416232 49.062880 26.363021 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) +# From original +# V KA CL d_V d_KA d_CL +# 6.252332 7.547072 4.240929 32.205996 47.014629 25.684326 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) +# V KA CL d_V d_KA d_CL +# 6.331614 8.009220 4.297905 32.351741 51.795028 26.386514 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) diff --git a/inst/poped/ex.1.b.PK.1.comp.oral.md.re-parameterize.babelmixr2.R b/inst/poped/ex.1.b.PK.1.comp.oral.md.re-parameterize.babelmixr2.R index a272fd64..0c14759e 100644 --- a/inst/poped/ex.1.b.PK.1.comp.oral.md.re-parameterize.babelmixr2.R +++ b/inst/poped/ex.1.b.PK.1.comp.oral.md.re-parameterize.babelmixr2.R @@ -71,6 +71,9 @@ output <- poped_optim(babel.db, opt_xt =TRUE) # Evaluate optimization results summary(output) +# from original +# V KA KE d_V d_KA d_KE +# 6.303480 7.899813 5.750915 29.156362 43.794939 33.388714 get_rse(output$FIM,output$poped.db) plot_model_prediction(output$poped.db) @@ -90,6 +93,9 @@ output_discrete <- poped_optim(babel.db.discrete, opt_xt=T) summary(output_discrete) +# from original +# V KA KE d_V d_KA d_KE +# 6.367801 8.230615 5.835965 29.067957 45.786136 33.418356 get_rse(output_discrete$FIM,output_discrete$poped.db) plot_model_prediction(output_discrete$poped.db) From 7a1cabd05422a07f1d75694af8557f2d0fe833e5 Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Wed, 23 Oct 2024 12:25:33 +0200 Subject: [PATCH 19/24] turned parallelization off --- inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R b/inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R index b11af149..4c3f280c 100644 --- a/inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R +++ b/inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R @@ -3,8 +3,8 @@ ## for population pharmacokinetics-pharmacodynamics studies", ## Br. J. Clin. Pharm., 2014. -## Optimization using an additive + proportional reidual error to -## avoid sample times at very low concentrations (time 0 or very late samoples). +## Optimization using an additive + proportional residual error to +## avoid sample times at very low concentrations (time 0 or very late samples). library(babelmixr2) library(PopED) @@ -59,7 +59,7 @@ plot_model_prediction(babel.db,IPRED=T,DV=T) evaluate_design(babel.db) # RS+SG+LS optimization of sample times -output <- poped_optim(babel.db, opt_xt=T, parallel=T) +output <- poped_optim(babel.db, opt_xt=T, parallel=F) summary(output) From a8ef8abfc5fe6a4fd2dc9c752b90e764dd539d02 Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Wed, 23 Oct 2024 12:31:50 +0200 Subject: [PATCH 20/24] turn parallel off --- inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R | 1 + 1 file changed, 1 insertion(+) diff --git a/inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R b/inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R index 4c3f280c..45e0a27c 100644 --- a/inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R +++ b/inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R @@ -51,6 +51,7 @@ plot_model_prediction(babel.db) ## create plot of model with variability plot_model_prediction(babel.db,IPRED=T,DV=T) +# Original RSEs ## $rse ## CL V KA d_CL d_V d_KA SIGMA[1,1] SIGMA[2,2] ## 5.096246 3.031164 14.260384 29.761226 36.681388 26.748640 32.011719 25.637971 From 5c066cdecbda17233c8484dc6939010cc3be7a53 Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Wed, 23 Oct 2024 14:18:51 +0200 Subject: [PATCH 21/24] minor corrections --- .../ex.3.a.PKPD.1.comp.oral.md.imax.D-opt.babelmixr2.R | 6 +++++- .../ex.3.b.PKPD.1.comp.oral.md.imax.ED-opt.babelmixr2.R | 4 +--- src/RcppExports.cpp | 1 + 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/inst/poped/ex.3.a.PKPD.1.comp.oral.md.imax.D-opt.babelmixr2.R b/inst/poped/ex.3.a.PKPD.1.comp.oral.md.imax.D-opt.babelmixr2.R index 13c81421..9eb906b3 100644 --- a/inst/poped/ex.3.a.PKPD.1.comp.oral.md.imax.D-opt.babelmixr2.R +++ b/inst/poped/ex.3.a.PKPD.1.comp.oral.md.imax.D-opt.babelmixr2.R @@ -114,10 +114,14 @@ evaluate_design(babel.db) shrinkage(babel.db) # Optimization -output <- poped_optim(babel.db, opt_xt = T, parallel = T) +output <- poped_optim(babel.db, opt_xt = T, parallel = F) summary(output) +# Original: +# V KA CL E0 IMAX IC50 d_V d_KA d_CL d_E0 +# 7.378959 9.224199 4.363698 4.923343 5.945328 24.169815 37.404330 58.360336 27.113994 22.005904 + get_rse(output$FIM,output$poped.db) plot_model_prediction(output$poped.db,facet_scales="free") diff --git a/inst/poped/ex.3.b.PKPD.1.comp.oral.md.imax.ED-opt.babelmixr2.R b/inst/poped/ex.3.b.PKPD.1.comp.oral.md.imax.ED-opt.babelmixr2.R index 0b4c25cc..26acf3a9 100644 --- a/inst/poped/ex.3.b.PKPD.1.comp.oral.md.imax.ED-opt.babelmixr2.R +++ b/inst/poped/ex.3.b.PKPD.1.comp.oral.md.imax.ED-opt.babelmixr2.R @@ -1,8 +1,6 @@ library(babelmixr2) -devtools::load_all() - library(PopED) ##-- Model: One comp first order absorption + inhibitory imax @@ -82,7 +80,7 @@ babel.db <- nlmixr2(f, e, "poped", bpop_vals_ed <- babel.db$parameters$bpop -bpop_vals_ed["tIC50",1] <- 1 # normal distrtibution +bpop_vals_ed["tIC50",1] <- 1 # normal distribution bpop_vals_ed["tIC50",3] <- (bpop_vals_ed["tIC50",2]*0.1)^2 bpop_vals_ed diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index abe37bda..b58c4624 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -1,6 +1,7 @@ // Generated by using Rcpp::compileAttributes() -> do not edit by hand // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 +#include #include #include From cc499190a519d543938848a01d000f50140d5a5d Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Wed, 23 Oct 2024 15:23:37 +0200 Subject: [PATCH 22/24] some corrections added --- ...ex.12.covariate.distributions.babelmixr2.R | 27 ++++++++----------- inst/poped/ex.14.PK.IOV.babelmixr2.R | 18 +++++-------- inst/poped/ex.4.PKPD.1.comp.emax.babelmixr2.R | 16 +++++------ inst/poped/ex.5.PD.emax.hill.babelmixr2.R | 1 - .../poped/ex.6.PK.1.comp.oral.sd.babelmixr2.R | 1 - .../ex.7.PK.1.comp.maturation.babelmixr2.R | 1 - 6 files changed, 25 insertions(+), 39 deletions(-) diff --git a/inst/poped/ex.12.covariate.distributions.babelmixr2.R b/inst/poped/ex.12.covariate.distributions.babelmixr2.R index f5043d2f..cf65ee1d 100644 --- a/inst/poped/ex.12.covariate.distributions.babelmixr2.R +++ b/inst/poped/ex.12.covariate.distributions.babelmixr2.R @@ -45,10 +45,10 @@ babel.db <- nlmixr2(f, e, "poped", a=c(WT=70))) # We can create a plot of the model typical predictions: -plot_model_prediction(poped_db) +plot_model_prediction(babel.db) # And evaluate the initial design -evaluate_design(poped_db) +evaluate_design(babel.db) # We see that the covariate parameters can not be estimated # according to this design calculation (RSE of bpop[3]=0 and bpop[4]=0). @@ -82,20 +82,15 @@ evaluate_design(babel.db.2) nsim <- 10 rse_list <- c() for(i in 1:nsim){ - poped_db_tmp <- create.poped.database(ff_fun=mod_1, - fg_fun=par_1, - fError_fun=feps.add.prop, - groupsize=1, - m=50, - sigma=c(0.015,0.0015), - notfixed_sigma = c(1,0), - bpop=c(CL=3.8,V=20,WT_CL=0.75,WT_V=1), - d=c(CL=0.05,V=0.05), - xt=c( 1,2,4,6,8,24), - minxt=0, - maxxt=24, - bUseGrouped_xt=1, - a=lapply(rnorm(50, mean = 70, sd = 10), function(x) c(WT=x))) + poped_db_tmp <- nlmixr2(f, e, "poped", + control=popedControl( + groupsize=1, + m=50, + minxt=0, + maxxt=24, + bUseGrouped_xt = T, + a = lapply(rnorm(50, mean = 70, sd = 10), function(x) c(WT=x)) )) + rse_tmp <- evaluate_design(poped_db_tmp)$rse rse_list <- rbind(rse_list,rse_tmp) } diff --git a/inst/poped/ex.14.PK.IOV.babelmixr2.R b/inst/poped/ex.14.PK.IOV.babelmixr2.R index abf83cf2..5159cfec 100644 --- a/inst/poped/ex.14.PK.IOV.babelmixr2.R +++ b/inst/poped/ex.14.PK.IOV.babelmixr2.R @@ -2,12 +2,7 @@ library(babelmixr2) library(PopED) library(ggplot2) -# Does not work for now. PopED has a special way to handle IOV. This feature will be added in the future - -# -# N = floor(Time/TAU)+1; -# CL = CL_OCC_1; -# if(N>6) CL = CL_OCC_2; +# Does not work properly for now. PopED has a special way to handle IOV. This feature will be added in the future ## define the ODE f <- function() { @@ -22,21 +17,20 @@ f <- function() { eta.cl ~0.25^2 # IOV - theta.iov.cl <- c(0.09) + theta.iov.cl <- sqrt(c(0.09)) iov.cl1 <- fix(1) # iov.cl2 <- fix(1) # iov.cl3 <- fix (1) # as many as occasions as needed - prop.sd <- fix(sqrt(0.04)) add.sd <- fix(sqrt(5e-6)) }) model({ # No of doses - N = floor(t/TAU)+1 - OCC1 = ifelse(N>6, 1, 0) + N = floor(t/TAU) + 1 + OCC1 = ifelse(N > 6, 1, 0) iov <- theta.iov.cl * (OCC1 * iov.cl1) V<-tV*exp(eta.v) @@ -59,7 +53,7 @@ e <- et(list(c(0, 10), c(0, 10), c(240, 248), c(240, 248))) %>% - et(amt=1/0.9, ii=24, until=248,cmt="depot") %>% + et(amt=1, ii=24, until=248,cmt="depot") %>% as.data.frame() #xt e$time <- c(0, 1,2,8,240,245) @@ -83,7 +77,7 @@ plot_model_prediction(babel.db, PRED=F,IPRED=F, groupsize_sim = 1, IPRED.lines = T, alpha.IPRED.lines=0.6, sample.times = F -) + geom_vline(xintercept = 24*6,color="red") + coord_cartesian(ylim = c(0, 0.75)) +) + geom_vline(xintercept = 24*6,color="red") ## evaluate initial design evaluate_design(babel.db) diff --git a/inst/poped/ex.4.PKPD.1.comp.emax.babelmixr2.R b/inst/poped/ex.4.PKPD.1.comp.emax.babelmixr2.R index fbd2a0c4..5aed746f 100644 --- a/inst/poped/ex.4.PKPD.1.comp.emax.babelmixr2.R +++ b/inst/poped/ex.4.PKPD.1.comp.emax.babelmixr2.R @@ -46,9 +46,6 @@ f <- function() { }) } -# e <- et(c(1,3,8)) %>% -# as.data.frame() - event.table.optim.pk <- # amt is a placeholder here et(amt = 1) %>% @@ -66,9 +63,6 @@ event.table.optim.pd <- dvid = 'effect') e <- event.table.optim.pk %>% bind_rows(event.table.optim.pd) %>% arrange(id, time) -# event.table.optim2 <- event.table.optim %>% filter(id ==2) -str(e) - ## -- Define initial design and design space babel.db <- nlmixr2(f, e, "poped", @@ -91,8 +85,14 @@ evaluate_design(babel.db) shrinkage(babel.db) # Optimization of sample times and doses -output <- poped_optim(babel.db, opt_xt = T, opt_a = T, parallel = T,method = c("LS")) +output <- poped_optim(babel.db, opt_xt = T, opt_a = T, parallel = F,method = c("LS")) + +# Original +# CL V E0 EMAX EC50 d_CL d_V d_E0 d_EC50 SIGMA[1,1] +# 4.713022 4.830290 3.023423 2.471232 9.089469 23.842058 25.111412 20.772844 89.433096 12.217638 +# SIGMA[2,2] +# 12.405893 -get_rse(output$FIM,output$babel.db) +get_rse(output$FIM,output$poped.db) plot_model_prediction(output$babel.db,facet_scales="free") diff --git a/inst/poped/ex.5.PD.emax.hill.babelmixr2.R b/inst/poped/ex.5.PD.emax.hill.babelmixr2.R index fafa284f..12332335 100644 --- a/inst/poped/ex.5.PD.emax.hill.babelmixr2.R +++ b/inst/poped/ex.5.PD.emax.hill.babelmixr2.R @@ -1,4 +1,3 @@ - library(PopED) library(babelmixr2) diff --git a/inst/poped/ex.6.PK.1.comp.oral.sd.babelmixr2.R b/inst/poped/ex.6.PK.1.comp.oral.sd.babelmixr2.R index f1133232..6deb9290 100644 --- a/inst/poped/ex.6.PK.1.comp.oral.sd.babelmixr2.R +++ b/inst/poped/ex.6.PK.1.comp.oral.sd.babelmixr2.R @@ -1,5 +1,4 @@ library(PopED) -# library(deSolve) library(babelmixr2) ##-- Model: One comp first order absorption diff --git a/inst/poped/ex.7.PK.1.comp.maturation.babelmixr2.R b/inst/poped/ex.7.PK.1.comp.maturation.babelmixr2.R index 9d691c60..f1875c08 100644 --- a/inst/poped/ex.7.PK.1.comp.maturation.babelmixr2.R +++ b/inst/poped/ex.7.PK.1.comp.maturation.babelmixr2.R @@ -1,4 +1,3 @@ - library(babelmixr2) library(PopED) From 8caa1b0a4ea9d942988e0e9e73127412aab8ecfb Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Tue, 29 Oct 2024 13:56:52 +0100 Subject: [PATCH 23/24] Added study 1 in gibiansky,JPKPD,2012 --- ....tmdd_qss_one_target_compiled.babelmixr2.R | 54 +++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/inst/poped/ex.8.tmdd_qss_one_target_compiled.babelmixr2.R b/inst/poped/ex.8.tmdd_qss_one_target_compiled.babelmixr2.R index 21c90d7d..13dfac15 100644 --- a/inst/poped/ex.8.tmdd_qss_one_target_compiled.babelmixr2.R +++ b/inst/poped/ex.8.tmdd_qss_one_target_compiled.babelmixr2.R @@ -87,6 +87,50 @@ e <- rbind(e1, e2) %>% e <- rbind(e0, e) + + + + +################################################# +# for study 1 in gibiansky,JPKPD,2012 table 2 +################################################# +db.1 <- nlmixr2(f, e, "poped", + control=popedControl( + groupsize=6, + a=list(c(ID=1, DOSE=100, SC_FLAG=0), + c(ID=1, DOSE=300, SC_FLAG=0), + c(ID=1, DOSE=600, SC_FLAG=0), + c(ID=1, DOSE=1000, SC_FLAG=1)), + discrete_a = list(DOSE=seq(100,1000,by=100), + SC_FLAG=c(0,1)), + )) + +plot_model_prediction(db.1,facet_scales="free") + +tic(); eval <- evaluate_design(db.1); toc() + + +# This is a longer running example, so the choice to check times for +# changes each time may not be ideal. + +# For this first example, we will not check time differences, which +# speeds up solving a bit + +db <- babel.poped.database(db, optTime=FALSE) + +tic();e1 <- evaluate_design(db);toc() + +# Original example: +# $rse +# CL V1 Q V2 FAVAIL KA R0 KSSS KDEG KINT d_CL +# 7.125675 6.867018 7.408125 10.435998 11.471400 19.592490 8.408247 11.028650 8.175673 7.343506 32.695140 +# d_V1 d_Q d_V2 d_FAVAIL d_KA d_R0 d_KSSS d_KDEG d_KINT SIGMA[1,1] SIGMA[2,2] +# 33.548134 74.864908 84.903024 98.434131 78.750729 36.359732 49.607416 66.478853 49.674898 8.935624 9.676856 + +################################################# +# for study 1 + 2 in gibiansky,JPKPD,2012 table 2 +################################################# + db <- nlmixr2(f, e, "poped", control=popedControl( groupsize=rbind(6,6,6,6,100,100), @@ -128,6 +172,16 @@ db <- babel.poped.database(db, optTime=TRUE) tic();e2 <- evaluate_design(db);toc() +# now optimization in parallel for unix/mac +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine +output <- poped_optim(db,opt_xt = F, opt_a = T, parallel=T, method = c("LS")) + +# optimization for windows +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine +# output <- poped_optim(poped.db.2,opt_xt = F, opt_a = T, parallel=T, method = c("LS"), dlls = c('tmdd_qss_one_target')) +plot_model_prediction(output$poped.db,facet_scales="free") # You can see the differences between the two with microbenchmark. m1 <- microbenchmark::microbenchmark(compareTime=evaluate_design(db), times=25) From 2754fa468d7696126c4215fa8a88ac216570aa41 Mon Sep 17 00:00:00 2001 From: Theodoros Papathanasiou Date: Tue, 29 Oct 2024 13:57:57 +0100 Subject: [PATCH 24/24] added parallel option with Windows message --- .../PopED_output_summary_D_cont_opt_1.txt | 229 ++++++++++++++++++ ...x.1.a.PK.1.comp.oral.md.intro.babelmixr2.R | 15 +- ....comp.oral.md.re-parameterize.babelmixr2.R | 12 +- ...K.1.comp.oral.md.ODE.compiled.babelmixr2.R | 14 +- .../ex.2.b.warfarin.optimize.babelmixr2.R | 10 +- .../ex.2.c.warfarin.ODE.compiled.babelmixr2.R | 4 +- inst/poped/ex.2.d.warfarin.ED.babelmixr2.R | 9 +- inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R | 4 +- ...KPD.1.comp.oral.md.imax.D-opt.babelmixr2.R | 4 +- ...PD.1.comp.oral.md.imax.ED-opt.babelmixr2.R | 2 + inst/poped/ex.4.PKPD.1.comp.emax.babelmixr2.R | 4 +- inst/poped/ex.5.PD.emax.hill.babelmixr2.R | 2 + .../poped/ex.6.PK.1.comp.oral.sd.babelmixr2.R | 2 + ...K.2.comp.oral.md.ode.compiled.babelmixr2.R | 4 +- 14 files changed, 294 insertions(+), 21 deletions(-) create mode 100644 inst/poped/PopED_output_summary_D_cont_opt_1.txt diff --git a/inst/poped/PopED_output_summary_D_cont_opt_1.txt b/inst/poped/PopED_output_summary_D_cont_opt_1.txt new file mode 100644 index 00000000..c5ff60fa --- /dev/null +++ b/inst/poped/PopED_output_summary_D_cont_opt_1.txt @@ -0,0 +1,229 @@ +PopED Results + + 2024-10-24 00:44:14.791373 + +============================================================================== +Model description : PopED model + +Model Sizes : +Number of individual model parameters g[j] : Ng = 10 +Number of population model fixed parameters bpop[j] : Nbpop = 8 +Number of population model random effects parameters b[j] : Nb = 8 + +Typical Population Parameters: +bpop[1]: 3.908 +bpop[2]: -2.188 +bpop[3]: 0.558 +bpop[4]: -0.1864 +bpop[5]: 2.261 +bpop[6]: 0.2105 +bpop[7]: 3.708 +bpop[8]: -0.7089 + +Between Subject Variability matrix D (variance units) +0.0625 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 +0.0000 0.0625 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 +0.0000 0.0000 0.0625 0.0000 0.0000 0.0000 0.0000 0.0000 +0.0000 0.0000 0.0000 0.0625 0.0000 0.0000 0.0000 0.0000 +0.0000 0.0000 0.0000 0.0000 0.0625 0.0000 0.0000 0.0000 +0.0000 0.0000 0.0000 0.0000 0.0000 0.0625 0.0000 0.0000 +0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0625 0.0000 +0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0625 + +Diagonal Elements of D [sqrt(param)]: +D[1,1]: 0.0625 [ 0.25] +D[2,2]: 0.0625 [ 0.25] +D[3,3]: 0.0625 [ 0.25] +D[4,4]: 0.0625 [ 0.25] +D[5,5]: 0.0625 [ 0.25] +D[6,6]: 0.0625 [ 0.25] +D[7,7]: 0.0625 [ 0.25] +D[8,8]: 0.0625 [ 0.25] + +Residual Unexplained Variability matrix SIGMA (variance units) : +0.00927944 0.00000000 0.00000000 0.00000000 +0.00000000 0.00100000 0.00000000 0.00000000 +0.00000000 0.00000000 0.02246920 0.00000000 +0.00000000 0.00000000 0.00000000 0.00100000 + +Diagonal Elements of SIGMA [sqrt(param)]: +SIGMA[1,1]: 0.009279 [0.09633] +SIGMA[2,2]: 0.001 [0.03162] +SIGMA[3,3]: 0.02247 [0.1499] +SIGMA[4,4]: 0.001 [0.03162] + +============================================================================== +Experiment description (design and design space) + +Number of individuals: 100 +Number of groups (individuals with same design): 2 +Number of individuals per group: + Group 1: 50 + Group 2: 50 +Number of samples per group: + Number of discrete experimental variables: 0 +Number of model covariates: 2 + +Initial Sampling Schedule +Group 1: Model 1: 0.02 0.25 1 3 10 +Group 1: Model 2: 0.02 0.25 1 3 10 +Group 2: Model 1: 1 7 15 28 42 +Group 2: Model 2: 1 7 15 28 42 + +Minimum allowed sampling values +Group 1: Model 1: 0.02 0.25 1 3 10 +Group 1: Model 2: 0.02 0.25 1 3 10 +Group 2: Model 1: 1 7 15 28 42 +Group 2: Model 2: 1 7 15 28 42 + +Maximum allowed sampling values +Group 1: Model 1: 0.02 0.25 1 3 10 +Group 1: Model 2: 0.02 0.25 1 3 10 +Group 2: Model 1: 1 7 15 28 42 +Group 2: Model 2: 1 7 15 28 42 + +Covariates: +Group 1: 1 : 300 +Group 2: 2 : 10000 + +=============================================================================== +Initial design evaluation + +Initial OFV = 138.56 + +Efficiency criterion [usually defined as OFV^(1/npar)] = 2203.41 + +Initial design +expected relative standard error +(%RSE, rounded to nearest integer) + Parameter Values RSE_0 + tvc 3.91 1 + tk10 -2.19 2 + tk12 0.558 8 + tk21 -0.186 24 + tvm 2.26 2 + tkmc 0.21 48 + tk03 3.71 1 + tk30 -0.709 7 + d_eta.vc 0.0625 17 + d_eta.k10 0.0625 32 + d_eta.k12 0.0625 28 + d_eta.k21 0.0625 32 + d_eta.vm 0.0625 25 + d_eta.kmc 0.0625 102 + d_eta.k03 0.0625 21 + d_eta.k30 0.0625 32 + sig_var_eps1 0.00928 11 + sig_var_eps3 0.0225 18 + +============================================================================== +Criterion Specification + +OFV calculation for FIM: 4 + 1=Determinant of FIM, + 4=log determinant of FIM, + 6=determinant of interesting part of FIM (Ds) + +Approximation method: 0 + 0=FO, + 1=FOCE, + 2=FOCEI, + 3=FOI + +Fisher Information Matrix type: 1 + 0=Full FIM, + 1=Reduced FIM, + 2=weighted models, + 3=Loc models, + 4=reduced FIM with derivative of SD of sigma as pfim, + 5=FULL FIM parameterized with A,B,C matrices & derivative of variance, + 6=Calculate one model switch at a time, good for large matrices, + 7=Reduced FIM parameterized with A,B,C matrices & derivative of variance + +Design family: 1 + D-family design (1) or + ED-family design (0) + (with or without parameter uncertainty) + +============================================================================== +Optimization of design parameters + +* Optimize Sampling Schedule + +******************************* +Initial Value + OFV(mf) = 138.56 +******************************* + +RS - It. : 5 OFV : 138.56 +RS - It. : 10 OFV : 138.56 +RS - It. : 15 OFV : 138.56 +RS - It. : 20 OFV : 138.56 +RS - It. : 25 OFV : 138.56 +RS - It. : 30 OFV : 138.56 +RS - It. : 35 OFV : 138.56 +RS - It. : 40 OFV : 138.56 +RS - It. : 45 OFV : 138.56 +RS - It. : 50 OFV : 138.56 +RS - It. : 55 OFV : 138.56 +RS - It. : 60 OFV : 138.56 +RS - It. : 65 OFV : 138.56 +RS - It. : 70 OFV : 138.56 +RS - It. : 75 OFV : 138.56 +RS - It. : 80 OFV : 138.56 +RS - It. : 85 OFV : 138.56 +RS - It. : 90 OFV : 138.56 +RS - It. : 95 OFV : 138.56 +RS - It. : 100 OFV : 138.56 +RS - It. : 105 OFV : 138.56 +RS - It. : 110 OFV : 138.56 +RS - It. : 115 OFV : 138.56 +RS - It. : 120 OFV : 138.56 +RS - It. : 125 OFV : 138.56 +RS - It. : 130 OFV : 138.56 +RS - It. : 135 OFV : 138.56 +RS - It. : 140 OFV : 138.56 +RS - It. : 145 OFV : 138.56 +RS - It. : 150 OFV : 138.56 +RS - It. : 155 OFV : 138.56 +RS - It. : 160 OFV : 138.56 +RS - It. : 165 OFV : 138.56 +RS - It. : 170 OFV : 138.56 +RS - It. : 175 OFV : 138.56 +RS - It. : 180 OFV : 138.56 +RS - It. : 185 OFV : 138.56 +RS - It. : 190 OFV : 138.56 +RS - It. : 195 OFV : 138.56 +RS - It. : 200 OFV : 138.56 +RS - It. : 205 OFV : 138.56 +RS - It. : 210 OFV : 138.56 +RS - It. : 215 OFV : 138.56 +RS - It. : 220 OFV : 138.56 +RS - It. : 225 OFV : 138.56 +RS - It. : 230 OFV : 138.56 +RS - It. : 235 OFV : 138.56 +RS - It. : 240 OFV : 138.56 +RS - It. : 245 OFV : 138.56 +RS - It. : 250 OFV : 138.56 +RS - It. : 255 OFV : 138.56 +RS - It. : 260 OFV : 138.56 +RS - It. : 265 OFV : 138.56 +RS - It. : 270 OFV : 138.56 +RS - It. : 275 OFV : 138.56 +RS - It. : 280 OFV : 138.56 +RS - It. : 285 OFV : 138.56 +RS - It. : 290 OFV : 138.56 +RS - It. : 295 OFV : 138.56 +RS - It. : 300 OFV : 138.56 + +******************************* +RS Results + OFV(mf) = 138.56 + +Optimized Sampling Schedule +Group 1: Model 1: 0.02 0.25 1 3 10 +Group 1: Model 2: 0.02 0.25 1 3 10 +Group 2: Model 1: 1 7 15 28 42 +Group 2: Model 2: 1 7 15 28 42 +********************************* + diff --git a/inst/poped/ex.1.a.PK.1.comp.oral.md.intro.babelmixr2.R b/inst/poped/ex.1.a.PK.1.comp.oral.md.intro.babelmixr2.R index 9166729c..0418c3b2 100644 --- a/inst/poped/ex.1.a.PK.1.comp.oral.md.intro.babelmixr2.R +++ b/inst/poped/ex.1.a.PK.1.comp.oral.md.intro.babelmixr2.R @@ -42,7 +42,7 @@ e <- et(list(c(0, 10), c(240, 248))) %>% as.data.frame() -#xt +# PopED xt equivalent e$time <- c(1,2,8,240,245) @@ -84,8 +84,9 @@ evaluate_design(babel.db) shrinkage(babel.db) # Optimization of sample times - -output <- poped_optim(babel.db, opt_xt =TRUE) +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine +output <- poped_optim(babel.db, opt_xt =TRUE, parallel=TRUE) # Evaluate optimization results summary(output) @@ -98,7 +99,9 @@ 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) +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine +output_2 <- poped_optim(output$poped.db, opt_xt =TRUE, opt_a = TRUE, parallel = TRUE) summary(output_2) @@ -113,7 +116,9 @@ 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) +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine +output_discrete <- poped_optim(babel.db.discrete, opt_xt=T, parallel = TRUE) summary(output_discrete) diff --git a/inst/poped/ex.1.b.PK.1.comp.oral.md.re-parameterize.babelmixr2.R b/inst/poped/ex.1.b.PK.1.comp.oral.md.re-parameterize.babelmixr2.R index 0c14759e..cd77e395 100644 --- a/inst/poped/ex.1.b.PK.1.comp.oral.md.re-parameterize.babelmixr2.R +++ b/inst/poped/ex.1.b.PK.1.comp.oral.md.re-parameterize.babelmixr2.R @@ -66,7 +66,9 @@ evaluate_design(babel.db) shrinkage(babel.db) # Optimization of sample times -output <- poped_optim(babel.db, opt_xt =TRUE) +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine +output <- poped_optim(babel.db, opt_xt =TRUE, parallel=TRUE) # Evaluate optimization results summary(output) @@ -79,7 +81,9 @@ 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) +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine +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) @@ -89,7 +93,9 @@ 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) +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine +output_discrete <- poped_optim(babel.db.discrete, opt_xt=T, parallel=TRUE) summary(output_discrete) diff --git a/inst/poped/ex.1.c.PK.1.comp.oral.md.ODE.compiled.babelmixr2.R b/inst/poped/ex.1.c.PK.1.comp.oral.md.ODE.compiled.babelmixr2.R index 611da739..aeae688e 100644 --- a/inst/poped/ex.1.c.PK.1.comp.oral.md.ODE.compiled.babelmixr2.R +++ b/inst/poped/ex.1.c.PK.1.comp.oral.md.ODE.compiled.babelmixr2.R @@ -37,7 +37,7 @@ e <- et(list(c(0, 10), c(0, 10), c(240, 248), c(240, 248))) %>% - et(amt=1000, ii=24, until=248,cmt="depot") %>% + et(amt=1/0.9, ii=24, until=248,cmt="depot") %>% as.data.frame() #xt @@ -64,7 +64,9 @@ evaluate_design(babel.db) shrinkage(babel.db) # Optimization of sample times -output <- poped_optim(babel.db, opt_xt =TRUE) +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine +output <- poped_optim(babel.db, opt_xt =TRUE, parallel=TRUE, method = c("LS")) # Evaluate optimization results summary(output) @@ -75,7 +77,9 @@ 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) +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine +output_2 <- poped_optim(output$poped.db, opt_xt =TRUE, opt_a = TRUE, parallel=TRUE, method = c("LS")) summary(output_2) @@ -88,7 +92,9 @@ 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) +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine +output_discrete <- poped_optim(babel.db.discrete, opt_xt=T, parallel=TRUE, method = c("LS")) summary(output_discrete) diff --git a/inst/poped/ex.2.b.warfarin.optimize.babelmixr2.R b/inst/poped/ex.2.b.warfarin.optimize.babelmixr2.R index 1fbb2507..0d861cae 100644 --- a/inst/poped/ex.2.b.warfarin.optimize.babelmixr2.R +++ b/inst/poped/ex.2.b.warfarin.optimize.babelmixr2.R @@ -65,6 +65,8 @@ evaluate_design(babel.db) # below are a number of ways to optimize the problem # Optimization of sample times +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine output <- poped_optim(babel.db, opt_xt = T, parallel = T) get_rse(output$FIM,output$poped.db) plot_model_prediction(output$poped.db) @@ -74,6 +76,8 @@ plot_efficiency_of_windows(output$poped.db,xt_windows=0.5) # Optimization of DOSE and sampling times +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine output_D_T <- poped_optim(babel.db, opt_xt = T, opt_a = T) summary(output_D_T) @@ -107,7 +111,9 @@ babel.db.discrete <- nlmixr2(f, e, "poped", discrete_xt=list(0:120), discrete_a=list(seq(10, 100, 10)))) -output_discrete <- poped_optim(babel.db.discrete, opt_xt = T, opt_a = T) +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine +output_discrete <- poped_optim(babel.db.discrete, opt_xt = T, opt_a = T, parallel = T) get_rse(output_discrete$FIM,output_discrete$poped.db) @@ -115,6 +121,8 @@ plot_model_prediction(output_discrete$poped.db) # Optimization using a genetic algorithm +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine output_ga <- poped_optim(babel.db, opt_xt = T, parallel = T, method = "GA") summary(output_ga) diff --git a/inst/poped/ex.2.c.warfarin.ODE.compiled.babelmixr2.R b/inst/poped/ex.2.c.warfarin.ODE.compiled.babelmixr2.R index bb5659bb..8654aa86 100644 --- a/inst/poped/ex.2.c.warfarin.ODE.compiled.babelmixr2.R +++ b/inst/poped/ex.2.c.warfarin.ODE.compiled.babelmixr2.R @@ -56,5 +56,7 @@ plot_model_prediction(babel.db,IPRED=T,DV=T) ## evaluate initial design (much faster than pure R solution) tic(); design_ode_compiled <- evaluate_design(babel.db); toc() -## making optimization times more resonable +## making optimization times more reasonable +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine output <- poped_optim(babel.db, opt_xt =TRUE, parallel=TRUE, method = "LS") diff --git a/inst/poped/ex.2.d.warfarin.ED.babelmixr2.R b/inst/poped/ex.2.d.warfarin.ED.babelmixr2.R index 582759f0..4bfd4800 100644 --- a/inst/poped/ex.2.d.warfarin.ED.babelmixr2.R +++ b/inst/poped/ex.2.d.warfarin.ED.babelmixr2.R @@ -78,8 +78,9 @@ tic();evaluate_design(babel.db,d_switch=FALSE,ED_samp_size=20); toc() ## 4.991010 2.977982 14.014207 29.802546 36.711408 26.754059 31.477157 25.297312 ## optimization with line search search -## parallelization does not seem to work with babelmixr2 or example: -output_ls <- poped_optim(babel.db, opt_xt=T, parallel=F, method = "LS", d_switch=F, ED_samp_size=20) +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine +output_ls <- poped_optim(babel.db, opt_xt=T, parallel=T, method = "LS", d_switch=F, ED_samp_size=20) ## laplace does not seem to work with babelmixr2 or example: ## ED: E(det(FIM)) using Laplace approximation @@ -88,4 +89,6 @@ output_ls <- poped_optim(babel.db, opt_xt=T, parallel=F, method = "LS", d_switch tic(); evaluate_design(babel.db,d_switch=FALSE,use_laplace=TRUE); toc() ## optimization with Laplace -output_ls <- poped_optim(babel.db, opt_xt=T, parallel=F, method = "LS", d_switch=F, use_laplace=T) +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine +output_ls <- poped_optim(babel.db, opt_xt=T, parallel=T, method = "LS", d_switch=F, use_laplace=T) diff --git a/inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R b/inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R index 45e0a27c..044c4d33 100644 --- a/inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R +++ b/inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R @@ -60,7 +60,9 @@ plot_model_prediction(babel.db,IPRED=T,DV=T) evaluate_design(babel.db) # RS+SG+LS optimization of sample times -output <- poped_optim(babel.db, opt_xt=T, parallel=F) +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine +output <- poped_optim(babel.db, opt_xt=T, parallel=T) summary(output) diff --git a/inst/poped/ex.3.a.PKPD.1.comp.oral.md.imax.D-opt.babelmixr2.R b/inst/poped/ex.3.a.PKPD.1.comp.oral.md.imax.D-opt.babelmixr2.R index 9eb906b3..ae252ea8 100644 --- a/inst/poped/ex.3.a.PKPD.1.comp.oral.md.imax.D-opt.babelmixr2.R +++ b/inst/poped/ex.3.a.PKPD.1.comp.oral.md.imax.D-opt.babelmixr2.R @@ -114,7 +114,9 @@ evaluate_design(babel.db) shrinkage(babel.db) # Optimization -output <- poped_optim(babel.db, opt_xt = T, parallel = F) +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine +output <- poped_optim(babel.db, opt_xt = T, parallel = T) summary(output) diff --git a/inst/poped/ex.3.b.PKPD.1.comp.oral.md.imax.ED-opt.babelmixr2.R b/inst/poped/ex.3.b.PKPD.1.comp.oral.md.imax.ED-opt.babelmixr2.R index 26acf3a9..21ece8b9 100644 --- a/inst/poped/ex.3.b.PKPD.1.comp.oral.md.imax.ED-opt.babelmixr2.R +++ b/inst/poped/ex.3.b.PKPD.1.comp.oral.md.imax.ED-opt.babelmixr2.R @@ -96,6 +96,8 @@ output$E_fim ## optimization with line search +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine output <- poped_optim(babel.db, opt_xt = T, parallel = T, d_switch=F,ED_samp_size=20, method = "LS") diff --git a/inst/poped/ex.4.PKPD.1.comp.emax.babelmixr2.R b/inst/poped/ex.4.PKPD.1.comp.emax.babelmixr2.R index 5aed746f..2753ae24 100644 --- a/inst/poped/ex.4.PKPD.1.comp.emax.babelmixr2.R +++ b/inst/poped/ex.4.PKPD.1.comp.emax.babelmixr2.R @@ -85,7 +85,9 @@ evaluate_design(babel.db) shrinkage(babel.db) # Optimization of sample times and doses -output <- poped_optim(babel.db, opt_xt = T, opt_a = T, parallel = F,method = c("LS")) +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine +output <- poped_optim(babel.db, opt_xt = T, opt_a = T, parallel = T,method = c("LS")) # Original # CL V E0 EMAX EC50 d_CL d_V d_E0 d_EC50 SIGMA[1,1] diff --git a/inst/poped/ex.5.PD.emax.hill.babelmixr2.R b/inst/poped/ex.5.PD.emax.hill.babelmixr2.R index 12332335..dcbab421 100644 --- a/inst/poped/ex.5.PD.emax.hill.babelmixr2.R +++ b/inst/poped/ex.5.PD.emax.hill.babelmixr2.R @@ -50,6 +50,8 @@ plot1 + xlab("Dose") evaluate_design(babel.db) # Optimization of doses +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine output <- poped_optim(babel.db, opt_xt = T, parallel = T) summary(output) diff --git a/inst/poped/ex.6.PK.1.comp.oral.sd.babelmixr2.R b/inst/poped/ex.6.PK.1.comp.oral.sd.babelmixr2.R index 6deb9290..1c21c23d 100644 --- a/inst/poped/ex.6.PK.1.comp.oral.sd.babelmixr2.R +++ b/inst/poped/ex.6.PK.1.comp.oral.sd.babelmixr2.R @@ -167,6 +167,8 @@ evaluate_design(babel.db.3) # Optimization of sample times # different results for different parameterizations +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine output.1 <- poped_optim(poped.db.1,opt_xt=T,parallel=T) output.2 <- poped_optim(poped.db.2,opt_xt=T,parallel=T) diff --git a/inst/poped/ex.9.PK.2.comp.oral.md.ode.compiled.babelmixr2.R b/inst/poped/ex.9.PK.2.comp.oral.md.ode.compiled.babelmixr2.R index fe01125e..8dcb8b0f 100644 --- a/inst/poped/ex.9.PK.2.comp.oral.md.ode.compiled.babelmixr2.R +++ b/inst/poped/ex.9.PK.2.comp.oral.md.ode.compiled.babelmixr2.R @@ -71,5 +71,7 @@ tic(); (eval_compiled <- evaluate_design(babel.db)); toc() #' but a large difference in computation time (8 times faster with the compiled code) (eval_compiled$ofv-eval$ofv)/eval$ofv -#' making optimization times more resonable +#' making optimization times more reasonable +# Note: The parallel option does not work well with Windows machines at this moment. +# Please set parallel = FALSE if you are working on a Windows machine #output <- poped_optim(babel.db,opt_xt=T, opt_a=T, parallel=T)