Merge pull request #121 from nlmixr2/121-script-output-2designs

Script output for different designs
nlmixr2 · Sep 15, 2024 · dbc95d2 · dbc95d2
2 parents 2de560d + ade9e43
commit dbc95d2
Show file tree

Hide file tree

Showing 13 changed files with 1,452 additions and 194 deletions.
diff --git a/R/monolixReadData.R b/R/monolixReadData.R
@@ -408,7 +408,7 @@ rxUiGet.monolixCovariance <- function(x, ...) {
   .cov <- rxUiGet.monolixCovarianceEstimatesSA(x, ...)
   .ui <- x[[1]]
   .split <- .ui$getSplitMuModel
-  
+
   .muRef <- c(.split$pureMuRef, .split$taintMuRef)
   .sa <- TRUE
   if (is.null(.cov)) {
@@ -507,4 +507,3 @@ rxUiGet.monolixPreds <- function(x, ...) {
        individualAbs=.qai, popAbs=.qap,
        message=.msg)
 }
-
diff --git a/R/poped.R b/R/poped.R
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
@@ -0,0 +1,104 @@
+library(babelmixr2)
+library(PopED)
+
+##-- Model: One comp first order absorption
+## -- Analytic solution for both mutiple and single dosing
+
+f <- function() {
+  ini({
+    tV <- 72.8
+    tKa <- 0.25
+    tCl <- 3.75
+    tF <- fix(0.9)
+
+    eta.v ~ 0.09
+    eta.ka ~ 0.09
+    eta.cl ~0.25^2
+
+    prop.sd <- fix(sqrt(0.04))
+    add.sd <- fix(sqrt(5e-6))
+
+  })
+  model({
+    V<-tV*exp(eta.v)
+    KA<-tKa*exp(eta.ka)
+    CL<-tCl*exp(eta.cl)
+    Favail <- tF
+    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)
+  })
+}
+
+# minxt, maxxt
+e <- et(list(c(0, 10),
+             c(0, 10),
+             c(0, 10),
+             c(240, 248),
+             c(240, 248))) %>%
+  as.data.frame()
+
+#xt
+e$time <-  c(1,2,8,240,245)
+
+
+babel.db <- nlmixr2(f, e, "poped",
+                        popedControl(groupsize=20,
+                                     bUseGrouped_xt=TRUE,
+                                     a=list(c(DOSE=20,TAU=24),
+                                            c(DOSE=40, TAU=24)),
+                                     maxa=c(DOSE=200,TAU=24),
+                                     mina=c(DOSE=0,TAU=24)))
+
+##  create plot of model without variability
+plot_model_prediction(babel.db, model_num_points = 300)
+
+##  create plot of model with variability
+plot_model_prediction(babel.db, IPRED=T, DV=T, separate.groups=T, model_num_points = 300)
+
+## evaluate initial design
+evaluate_design(babel.db)
+
+shrinkage(babel.db)
+
+# Optimization of sample times
+output <- poped_optim(babel.db, opt_xt =TRUE)
+
+# Evaluate optimization results
+summary(output)
+
+get_rse(output$FIM,output$poped.db)
+
+plot_model_prediction(output$poped.db)
+
+
+# Optimization of sample times and doses
+output_2 <- poped_optim(output$poped.db, opt_xt =TRUE, opt_a = TRUE)
+
+summary(output_2)
+
+get_rse(output_2$FIM,output_2$poped.db)
+
+plot_model_prediction(output_2$poped.db)
+
+
+# Optimization of sample times with only integer time points in design space
+# faster than continuous optimization in this case
+babel.db.discrete <- create.poped.database(babel.db,discrete_xt = list(0:248))
+
+output_discrete <- poped_optim(babel.db.discrete, opt_xt=T)
+
+
+summary(output_discrete)
+
+get_rse(output_discrete$FIM,output_discrete$poped.db)
+
+plot_model_prediction(output_discrete$poped.db)
+
+
+# Efficiency of sampling windows
+plot_efficiency_of_windows(output_discrete$poped.db, xt_windows=1)
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
@@ -0,0 +1,99 @@
+## using libary models and reparameterizing the problen to KA, KE and V
+## optimization of dose and dose interval
+
+library(babelmixr2)
+
+library(PopED)
+
+f <- function() {
+  ini({
+    tV <- 72.8
+    tKa <- 0.25
+    tKe <- 3.75/72.8
+    tFavail <- fix(0.9)
+
+    eta.v ~ 0.09
+    eta.ka ~ 0.09
+    eta.ke ~ 0.25^2
+
+    prop.sd <- fix(sqrt(0.04))
+    add.sd <- fix(sqrt(5e-6))
+  })
+  model({
+    V <- tV*exp(eta.v)
+    KA <- tKa*exp(eta.ka)
+    KE <- tKe*exp(eta.ke)
+    Favail <- tFavail
+    N <- floor(time/TAU)+1
+    y <- (DOSE*Favail/V)*(KA/(KA - KE)) *
+      (exp(-KE * (time - (N - 1) * TAU)) * (1 - exp(-N * KE * TAU))/(1 - exp(-KE * TAU)) -
+         exp(-KA * (time - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU)))
+
+    y ~ prop(prop.sd) + add(add.sd)
+  })
+}
+
+# minxt, maxxt
+e <- et(list(c(0, 10),
+             c(0, 10),
+             c(0, 10),
+             c(240, 248),
+             c(240, 248))) %>%
+  as.data.frame()
+
+#xt
+e$time <-  c(1,2,8,240,245)
+
+
+babel.db <- nlmixr2(f, e, "poped",
+                    popedControl(groupsize=20,
+                                 bUseGrouped_xt=TRUE,
+                                 a=list(c(DOSE=20,TAU=24),
+                                        c(DOSE=40, TAU=24)),
+                                 maxa=c(DOSE=200,TAU=24),
+                                 mina=c(DOSE=0,TAU=24)))
+
+
+##  create plot of model without variability
+plot_model_prediction(babel.db)
+
+##  create plot of model with variability
+plot_model_prediction(babel.db,IPRED=T,DV=T,separate.groups=T)
+
+## evaluate initial design
+evaluate_design(babel.db)
+
+shrinkage(babel.db)
+
+# Optimization of sample times
+output <- poped_optim(babel.db, opt_xt =TRUE)
+
+# Evaluate optimization results
+summary(output)
+
+get_rse(output$FIM,output$poped.db)
+
+plot_model_prediction(output$poped.db)
+
+# Optimization of sample times, doses and dose intervals
+output_2 <- poped_optim(output$poped.db, opt_xt =TRUE, opt_a = TRUE)
+
+summary(output_2)
+get_rse(output_2$FIM,output_2$poped.db)
+plot_model_prediction(output_2$poped.db)
+
+# Optimization of sample times with only integer time points in design space
+# faster than continuous optimization in this case
+babel.db.discrete <- create.poped.database(babel.db,discrete_xt = list(0:248))
+
+output_discrete <- poped_optim(babel.db.discrete, opt_xt=T)
+
+summary(output_discrete)
+
+get_rse(output_discrete$FIM,output_discrete$poped.db)
+
+plot_model_prediction(output_discrete$poped.db)
+
+
+# Efficiency of sampling windows
+plot_efficiency_of_windows(output_discrete$poped.db, xt_windows=1)
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
@@ -0,0 +1,102 @@
+library(babelmixr2)
+library(PopED)
+
+
+## define the ODE
+f <- function() {
+  ini({
+    tV <- 72.8
+    tKa <- 0.25
+    tCl <- 3.75
+    tF <- fix(0.9)
+
+    eta.v ~ 0.09
+    eta.ka ~ 0.09
+    eta.cl ~0.25^2
+
+    prop.sd <- fix(sqrt(0.04))
+    add.sd <- fix(sqrt(5e-6))
+
+  })
+  model({
+    V<-tV*exp(eta.v)
+    KA<-tKa*exp(eta.ka)
+    CL<-tCl*exp(eta.cl)
+    Favail <- tF
+    d/dt(depot) <- -KA*depot
+    d/dt(central) <- KA*depot - (CL/V)*central
+    f(depot) <- Favail*DOSE
+    y <- central/V
+    y ~ prop(prop.sd) + add(add.sd)
+  })
+}
+
+# minxt, maxxt
+e <- et(list(c(0, 10),
+             c(0, 10),
+             c(0, 10),
+             c(240, 248),
+             c(240, 248))) %>%
+  et(amt=1000, ii=24, until=248,cmt="depot") %>%
+  as.data.frame()
+
+#xt
+e$time <-  c(0, 1,2,8,240,245)
+
+
+babel.db <- nlmixr2(f, e, "poped",
+                    popedControl(groupsize=20,
+                                 bUseGrouped_xt=TRUE,
+                                 a=list(c(DOSE=20,TAU=24),
+                                        c(DOSE=40, TAU=24)),
+                                 maxa=c(DOSE=200,TAU=24),
+                                 mina=c(DOSE=0,TAU=24)))
+
+##  create plot of model without variability
+plot_model_prediction(babel.db, model_num_points = 300)
+
+##  create plot of model with variability
+plot_model_prediction(babel.db, IPRED=T, DV=T, separate.groups=T, model_num_points = 300)
+
+## evaluate initial design
+evaluate_design(babel.db)
+
+shrinkage(babel.db)
+
+# Optimization of sample times
+output <- poped_optim(babel.db, opt_xt =TRUE)
+
+# Evaluate optimization results
+summary(output)
+
+get_rse(output$FIM,output$poped.db)
+
+plot_model_prediction(output$poped.db)
+
+
+# Optimization of sample times and doses
+output_2 <- poped_optim(output$poped.db, opt_xt =TRUE, opt_a = TRUE)
+
+summary(output_2)
+
+get_rse(output_2$FIM,output_2$poped.db)
+
+plot_model_prediction(output_2$poped.db)
+
+
+# Optimization of sample times with only integer time points in design space
+# faster than continuous optimization in this case
+babel.db.discrete <- create.poped.database(babel.db,discrete_xt = list(0:248))
+
+output_discrete <- poped_optim(babel.db.discrete, opt_xt=T)
+
+
+summary(output_discrete)
+
+get_rse(output_discrete$FIM,output_discrete$poped.db)
+
+plot_model_prediction(output_discrete$poped.db)
+
+
+# Efficiency of sampling windows
+plot_efficiency_of_windows(output_discrete$poped.db, xt_windows=1)
diff --git a/inst/poped/ex.2.a.warfarin.evaluate.babelmixr2.R b/inst/poped/ex.2.a.warfarin.evaluate.babelmixr2.R
@@ -0,0 +1,68 @@
+## 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(babelmixr2)
+library(PopED)
+
+##-- Model: One comp first order absorption
+
+f <- function() {
+  ini({
+    tCl <- 0.15
+    tV <- 8
+    tKA <- 1.0
+    tFavail <- fix(1)
+    eta.cl ~ 0.07
+    eta.v ~ 0.02
+    eta.ka ~ 0.6
+
+    prop.sd <- sqrt(0.01)
+  })
+  model({
+    CL <- tCl*exp(eta.cl)
+    V <- tV*exp(eta.v)
+    KA <- tKA*exp(eta.ka)
+    Favail <- tFavail
+    y <- (DOSE*Favail*KA/(V*(KA-CL/V)))*(exp(-CL/V*time)-exp(-KA*time))
+    y ~ prop(prop.sd)
+
+  })
+}
+
+e <-  et(c(0.5, 1,2,6,24,36,72,120)) %>%
+  as.data.frame()
+
+## -- Define initial design  and design space
+babel.db <- nlmixr2(f, e, "poped",
+                    control=popedControl(
+                      groupsize=32,
+                      minxt=0,
+                      maxxt=120,
+                      a=70))
+
+##  create plot of model without variability
+plot_model_prediction(babel.db)
+
+##  create plot of model with variability
+plot_model_prediction(babel.db,IPRED=T,DV=T)
+
+#########################################
+## NOTE All PopED output for residuals
+## (add or prop) are VARIANCES instead of
+## standard deviations!
+#########################################
+
+## get predictions from model
+model_prediction(babel.db)
+
+## evaluate initial design
+evaluate_design(babel.db)
+shrinkage(babel.db)
+
+## Evaluate with full FIM
+evaluate_design(babel.db, fim.calc.type=0)
+
+# Examine efficiency of sampling windows
+plot_efficiency_of_windows(babel.db,xt_windows=0.5)