From 8100ac1d5cb6df315f5c04e38cfa03cfa02982dc Mon Sep 17 00:00:00 2001
From: Matthew Fidler <matthew.fidler@gmail.com>
Date: Thu, 12 Sep 2024 15:12:31 -0500
Subject: [PATCH] Add example 2.e

---
 inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R | 62 ++++++++++++++++++++++
 1 file changed, 62 insertions(+)
 create mode 100644 inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R

diff --git a/inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R b/inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R
new file mode 100644
index 00000000..11232029
--- /dev/null
+++ b/inst/poped/ex.2.e.warfarin.Ds.babelmixr2.R
@@ -0,0 +1,62 @@
+## 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.
+
+## Optimization using an additive + proportional reidual error to
+##   avoid sample times at very low concentrations (time 0 or very late samoples).
+
+library(babelmixr2)
+library(PopED)
+
+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) # nlmixr2 uses sd
+    add.sd <- sqrt(0.25)
+  })
+  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) + add(add.sd)
+  })
+}
+# First define standard controler from nlmixr2
+e <-  et(c(0.5, 1,2,6,24,36,72,120)) %>%
+  as.data.frame()
+
+babel.db <- nlmixr2(f, e, "poped",
+                    popedControl(groupsize=32,
+                                 minxt=0,
+                                 maxxt=120,
+                                 a=70,
+                                 mina=0,
+                                 maxa=100,
+                                 # selecting important/unimportant
+                                 # parameters assumes Ds optimal design.
+                                 important=c("tCl", "tV", "tKa")))
+
+##  create plot of model without variability
+plot_model_prediction(babel.db)
+
+##  create plot of model with variability
+plot_model_prediction(poped.db,IPRED=T,DV=T)
+
+## evaluate initial design
+evaluate_design(babel.db)
+
+# RS+SG+LS optimization of sample times
+output <- poped_optim(babel.db, opt_xt=T, parallel=T)
+
+summary(output)
+
+plot_model_prediction(output$poped.db)