Let backtransform_odeparms handle nlmixr formation fractions

Also adapt summary.nlmixr.mmkin to correctly handle the way
formation fractions are translated to nlmixr

11 files changed, 243 insertions, 149 deletions

diff --git a/NAMESPACE b/NAMESPACE
index 0f61396d..aa40b570 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -62,6 +62,7 @@ export(f_time_norm_focus)
 export(get_deg_func)
 export(ilr)
 export(invilr)
+export(invtffm0)
 export(loftest)
 export(logistic.solution)
 export(lrtest)
@@ -101,6 +102,7 @@ export(saem)
 export(saemix_data)
 export(saemix_model)
 export(sigma_twocomp)
+export(tffm0)
 export(transform_odeparms)
 import(deSolve)
 import(graphics)
diff --git a/R/dimethenamid_2018.R b/R/dimethenamid_2018.R
index 76b98efe..6e0bda0c 100644
--- a/R/dimethenamid_2018.R
+++ b/R/dimethenamid_2018.R
@@ -31,6 +31,7 @@
 #' dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]])
 #' dmta_ds[["Elliot 1"]] <- NULL
 #' dmta_ds[["Elliot 2"]] <- NULL
+#' \dontrun{
 #' dfop_sfo3_plus <- mkinmod(
 #'   DMTA = mkinsub("DFOP", c("M23", "M27", "M31")),
 #'   M23 = mkinsub("SFO"),
@@ -42,12 +43,15 @@
 #'   list("DFOP-SFO3+" = dfop_sfo3_plus),
 #'   dmta_ds, quiet = TRUE, error_model = "tc")
 #' nlmixr_model(f_dmta_mkin_tc)
-#' f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem",
-#'   control = saemControl(print = 500))
-#' summary(f_dmta_nlmixr_saem)
-#' plot(f_dmta_nlmixr_saem)
 #' f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei",
-#'   control = foceiControl(print = 500))
+#'   control = nlmixr::foceiControl(print = 500))
 #' summary(f_dmta_nlmixr_focei)
 #' plot(f_dmta_nlmixr_focei)
+#' # saem has a problem with this model/data combination, maybe because of the
+#' # overparameterised error model, to be investigated
+#' #f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem",
+#' #  control = saemControl(print = 500))
+#' #summary(f_dmta_nlmixr_saem)
+#' #plot(f_dmta_nlmixr_saem)
+#' }
 "dimethenamid_2018"
diff --git a/R/summary.nlmixr.mmkin.R b/R/summary.nlmixr.mmkin.R
index f2d7c607..a023f319 100644
--- a/R/summary.nlmixr.mmkin.R
+++ b/R/summary.nlmixr.mmkin.R
@@ -85,11 +85,11 @@ summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes
 
   mod_vars <- names(object$mkinmod$diffs)
 
-  pnames <- names(object$mean_dp_start)
-  np <- length(pnames)
-
   conf.int <- confint(object$nm)
-  confint_trans <- as.matrix(conf.int[pnames, c(1, 3, 4)])
+  dpnames <- setdiff(rownames(conf.int), names(object$mean_ep_start))
+  ndp <- length(dpnames)
+
+  confint_trans <- as.matrix(conf.int[dpnames, c(1, 3, 4)])
   colnames(confint_trans) <- c("est.", "lower", "upper")
 
   bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod,
@@ -100,7 +100,7 @@ summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes
   # with the exception of sets of formation fractions (single fractions are OK).
   f_names_skip <- character(0)
   for (box in mod_vars) { # Figure out sets of fractions to skip
-    f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE)
+    f_names <- grep(paste("^f", box, sep = "_"), dpnames, value = TRUE)
     n_paths <- length(f_names)
     if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names)
   }
@@ -109,7 +109,7 @@ summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes
     dimnames = list(bpnames, colnames(confint_trans)))
   confint_back[, "est."] <- bp
 
-  for (pname in pnames) {
+  for (pname in dpnames) {
     if (!pname %in% f_names_skip) {
       par.lower <- confint_trans[pname, "lower"]
       par.upper <- confint_trans[pname, "upper"]
@@ -131,7 +131,7 @@ summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes
   object$corFixed <- array(
     t(varFix/stdFix)/stdFix,
     dim(varFix),
-    list(pnames, pnames))
+    list(dpnames, dpnames))
 
   object$confint_trans <- confint_trans
   object$confint_back <- confint_back
diff --git a/R/tffm0.R b/R/tffm0.R
index 25787962..bb5f4cf5 100644
--- a/R/tffm0.R
+++ b/R/tffm0.R
@@ -13,7 +13,8 @@
 #'
 #' @param ff Vector of untransformed formation fractions. The sum
 #'   must be smaller or equal to one
-#' @param ff_trans
+#' @param ff_trans Vector of transformed formation fractions that can be
+#'   restricted to the interval from 0 to 1
 #' @return A vector of the transformed formation fractions
 #' @export
 #' @examples
@@ -33,7 +34,8 @@ tffm0 <- function(ff) {
   return(res)
 }
 #' @rdname tffm0
-#' @return
+#' @export
+#' @return A vector of backtransformed formation fractions for natural use in degradation models
 invtffm0 <- function(ff_trans) {
   n <- length(ff_trans)
   res <- numeric(n)
diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R
index 4fe4e5c2..174e7c2d 100644
--- a/R/transform_odeparms.R
+++ b/R/transform_odeparms.R
@@ -229,13 +229,18 @@ backtransform_odeparms <- function(transparms, mkinmod,
     if (length(trans_f) > 0) {
       if(transform_fractions) {
         if (any(grepl("qlogis", names(trans_f)))) {
-          parms[f_names] <- plogis(trans_f)
+          f_tmp  <- plogis(trans_f)
+          if (any(grepl("_tffm0_.*_qlogis$", names(f_tmp)))) {
+            parms[f_names] <- invtffm0(f_tmp)
+          } else {
+            parms[f_names] <- f_tmp
+          }
         } else {
-          f <- invilr(trans_f)
+          f_tmp <- invilr(trans_f)
           if (spec[[box]]$sink) {
-            parms[f_names] <- f[1:length(f)-1]
+            parms[f_names] <- f_tmp[1:length(f_tmp)-1]
           } else {
-            parms[f_names] <- f
+            parms[f_names] <- f_tmp
           }
         }
       } else {
diff --git a/check.log b/check.log
index 2627695d..df344926 100644
--- a/check.log
+++ b/check.log
@@ -57,10 +57,7 @@ Maintainer: ‘Johannes Ranke <jranke@uni-bremen.de>’
 * checking data for ASCII and uncompressed saves ... OK
 * checking installed files from ‘inst/doc’ ... OK
 * checking files in ‘vignettes’ ... OK
-* checking examples ... NOTE
-Examples with CPU (user + system) or elapsed time > 5s
-              user system elapsed
-nlmixr.mmkin 8.129  0.375   5.384
+* checking examples ... OK
 * checking for unstated dependencies in ‘tests’ ... OK
 * checking tests ... SKIPPED
 * checking for unstated dependencies in vignettes ... OK
@@ -71,9 +68,5 @@ nlmixr.mmkin 8.129  0.375   5.384
 * checking for detritus in the temp directory ... OK
 * DONE
 
-Status: 1 NOTE
-See
-  ‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’
-for details.
-
+Status: OK
 
diff --git a/man/dimethenamid_2018.Rd b/man/dimethenamid_2018.Rd
index b6f761e8..93fcad26 100644
--- a/man/dimethenamid_2018.Rd
+++ b/man/dimethenamid_2018.Rd
@@ -31,5 +31,41 @@ specific pieces of information in the comments.
 }
 \examples{
 print(dimethenamid_2018)
+dmta_ds <- lapply(1:8, function(i) {
+  ds_i <- dimethenamid_2018$ds[[i]]$data
+  ds_i[ds_i$name == "DMTAP", "name"] <-  "DMTA"
+  ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i]
+  ds_i
+})
+names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title)
+dmta_ds[["Borstel"]] <- rbind(dmta_ds[["Borstel 1"]], dmta_ds[["Borstel 2"]])
+dmta_ds[["Borstel 1"]] <- NULL
+dmta_ds[["Borstel 2"]] <- NULL
+dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]])
+dmta_ds[["Elliot 1"]] <- NULL
+dmta_ds[["Elliot 2"]] <- NULL
+\dontrun{
+dfop_sfo3_plus <- mkinmod(
+  DMTA = mkinsub("DFOP", c("M23", "M27", "M31")),
+  M23 = mkinsub("SFO"),
+  M27 = mkinsub("SFO"),
+  M31 = mkinsub("SFO", "M27", sink = FALSE),
+  quiet = TRUE
+)
+f_dmta_mkin_tc <- mmkin(
+  list("DFOP-SFO3+" = dfop_sfo3_plus),
+  dmta_ds, quiet = TRUE, error_model = "tc")
+nlmixr_model(f_dmta_mkin_tc)
+f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei",
+  control = nlmixr::foceiControl(print = 500))
+summary(f_dmta_nlmixr_focei)
+plot(f_dmta_nlmixr_focei)
+# saem has a problem with this model/data combination, maybe because of the
+# overparameterised error model, to be investigated
+#f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem",
+#  control = saemControl(print = 500))
+#summary(f_dmta_nlmixr_saem)
+#plot(f_dmta_nlmixr_saem)
+}
 }
 \keyword{datasets}
diff --git a/man/nlmixr.mmkin.Rd b/man/nlmixr.mmkin.Rd
index 4ab30272..0f4f41a2 100644
--- a/man/nlmixr.mmkin.Rd
+++ b/man/nlmixr.mmkin.Rd
@@ -16,6 +16,8 @@
   error_model = object[[1]]$err_mod,
   test_log_parms = TRUE,
   conf.level = 0.6,
+  degparms_start = "auto",
+  eta_start = "auto",
   ...,
   save = NULL,
   envir = parent.frame()
@@ -27,7 +29,8 @@ nlmixr_model(
   object,
   est = c("saem", "focei"),
   degparms_start = "auto",
-  test_log_parms = FALSE,
+  eta_start = "auto",
+  test_log_parms = TRUE,
   conf.level = 0.6,
   error_model = object[[1]]$err_mod,
   add_attributes = FALSE
@@ -58,6 +61,13 @@ when calculating mean degradation parameters using \link{mean_degparms}.}
 \item{conf.level}{Possibility to adjust the required confidence level
 for parameter that are tested if requested by 'test_log_parms'.}
 
+\item{degparms_start}{Parameter values given as a named numeric vector will
+be used to override the starting values obtained from the 'mmkin' object.}
+
+\item{eta_start}{Standard deviations on the transformed scale given as a
+named numeric vector will be used to override the starting values obtained
+from the 'mmkin' object.}
+
 \item{\dots}{Passed to \link{nlmixr_model}}
 
 \item{save}{Passed to \link[nlmixr:nlmixr]{nlmixr::nlmixr}}
@@ -68,9 +78,6 @@ for parameter that are tested if requested by 'test_log_parms'.}
 
 \item{digits}{Number of digits to use for printing}
 
-\item{degparms_start}{Parameter values given as a named numeric vector will
-be used to override the starting values obtained from the 'mmkin' object.}
-
 \item{add_attributes}{Should the starting values used for degradation model
 parameters and their distribution and for the error model parameters
 be returned as attributes?}
diff --git a/man/tffm0.Rd b/man/tffm0.Rd
new file mode 100644
index 00000000..46978d5e
--- /dev/null
+++ b/man/tffm0.Rd
@@ -0,0 +1,42 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/tffm0.R
+\name{tffm0}
+\alias{tffm0}
+\alias{invtffm0}
+\title{Transform formation fractions as in the first published mkin version}
+\usage{
+tffm0(ff)
+
+invtffm0(ff_trans)
+}
+\arguments{
+\item{ff}{Vector of untransformed formation fractions. The sum
+must be smaller or equal to one}
+
+\item{ff_trans}{Vector of transformed formation fractions that can be
+restricted to the interval from 0 to 1}
+}
+\value{
+A vector of the transformed formation fractions
+
+A vector of backtransformed formation fractions for natural use in degradation models
+}
+\description{
+The transformed fractions can be restricted between 0 and 1 in model
+optimisations. Therefore this transformation was used originally in mkin. It
+was later replaced by the \link{ilr} transformation because the ilr transformed
+fractions can assumed to follow normal distribution. As the ilr
+transformation is not available in \link{RxODE} and can therefore not be used in
+the nlmixr modelling language, this transformation is currently used for
+translating mkin models with formation fractions to more than one target
+compartment for fitting with nlmixr in \link{nlmixr_model}. However,
+this implementation cannot be used there, as it is not accessible
+from RxODE.
+}
+\examples{
+ff_example <- c(
+  0.10983681, 0.09035905, 0.08399383
+)
+ff_example_trans <- tffm0(ff_example)
+invtffm0(ff_example_trans)
+}
diff --git a/test.log b/test.log
index f2a60729..6ef8191f 100644
--- a/test.log
+++ b/test.log
@@ -3,17 +3,17 @@ Loading required package: parallel
 ℹ Testing mkin
 ✔ |  OK F W S | Context
 ✔ |   5       | AIC calculation
-✔ |   5       | Analytical solutions for coupled models [3.3 s]
+✔ |   5       | Analytical solutions for coupled models [3.5 s]
 ✔ |   5       | Calculation of Akaike weights
 ✔ |   2       | Export dataset for reading into CAKE
-✔ |  12       | Confidence intervals and p-values [1.3 s]
+✔ |  12       | Confidence intervals and p-values [1.0 s]
 ✔ |  14       | Error model fitting [4.7 s]
 ✔ |   5       | Time step normalisation
-✔ |   4       | Calculation of FOCUS chi2 error levels [0.5 s]
+✔ |   4       | Calculation of FOCUS chi2 error levels [0.6 s]
 ✔ |  14       | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.8 s]
 ✔ |   4       | Test fitting the decline of metabolites from their maximum [0.3 s]
 ✔ |   1       | Fitting the logistic model [0.2 s]
-✔ |  35     1 | Nonlinear mixed-effects models [27.1 s]
+✔ |  35     1 | Nonlinear mixed-effects models [26.9 s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip (test_mixed.R:161:3): saem results are reproducible for biphasic fits
 Reason: Fitting with saemix takes around 10 minutes when using deSolve
@@ -21,33 +21,36 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve
 ✔ |   2       | Test dataset classes mkinds and mkindsg
 ✔ |  10       | Special cases of mkinfit calls [0.4 s]
 ✔ |   1       | mkinfit features [0.3 s]
-✔ |   8       | mkinmod model generation and printing [0.3 s]
+✔ |   8       | mkinmod model generation and printing [0.2 s]
 ✔ |   3       | Model predictions with mkinpredict [0.3 s]
-✔ |  16       | Evaluations according to 2015 NAFTA guidance [1.7 s]
-✔ |   9       | Nonlinear mixed-effects models with nlme [8.1 s]
-✖ |  14 2     | Plotting [1.9 s]
+✔ |  14     2 | Evaluations according to 2015 NAFTA guidance [1.3 s]
 ────────────────────────────────────────────────────────────────────────────────
-Failure (test_plot.R:40:5): Plotting mkinfit, mmkin and mixed model objects is reproducible
-Figures don't match: mixed-model-fit-for-saem-object-with-saemix-transformations.svg
-
-
-Failure (test_plot.R:55:5): Plotting mkinfit, mmkin and mixed model objects is reproducible
-Figures don't match: mixed-model-fit-for-saem-object-with-mkin-transformations.svg
+Skip (test_nafta.R:25:5): Test data from Appendix B are correctly evaluated
+Reason: getRversion() >= "4.1.0" is TRUE
 
+Skip (test_nafta.R:53:5): Test data from Appendix D are correctly evaluated
+Reason: getRversion() >= "4.1.0" is TRUE
+────────────────────────────────────────────────────────────────────────────────
+✔ |   9       | Nonlinear mixed-effects models with nlme [8.1 s]
+✔ |   0     1 | Plotting [0.8 s]
+────────────────────────────────────────────────────────────────────────────────
+Skip (test_plot.R:18:3): Plotting mkinfit, mmkin and mixed model objects is reproducible
+Reason: getRversion() >= "4.1.0" is TRUE
 ────────────────────────────────────────────────────────────────────────────────
 ✔ |   4       | Residuals extracted from mkinfit models
 ✔ |   2       | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5 s]
-✔ |   7       | Fitting the SFORB model [3.9 s]
+✔ |   7       | Fitting the SFORB model [3.8 s]
 ✔ |   1       | Summaries of old mkinfit objects
 ✔ |   4       | Summary [0.1 s]
-✔ |   4       | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2 s]
-✔ |   9       | Hypothesis tests [8.2 s]
-✔ |   4       | Calculation of maximum time weighted average concentrations (TWAs) [2.4 s]
+✔ |   4       | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3 s]
+✔ |   9       | Hypothesis tests [8.5 s]
+✔ |   4       | Calculation of maximum time weighted average concentrations (TWAs) [2.2 s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 70.0 s
+Duration: 68.2 s
 
 ── Skipped tests  ──────────────────────────────────────────────────────────────
-● Fitting with saemix takes around 10 minutes when using deSolve (1)
+• Fitting with saemix takes around 10 minutes when using deSolve (1)
+• getRversion() >= "4.1.0" is TRUE (3)
 
-[ FAIL 2 | WARN 0 | SKIP 1 | PASS 204 ]
+[ FAIL 0 | WARN 0 | SKIP 4 | PASS 188 ]
diff --git a/tests/testthat/test_nlmixr.r b/tests/testthat/test_nlmixr.r
index e3bd3d66..dcbb50ac 100644
--- a/tests/testthat/test_nlmixr.r
+++ b/tests/testthat/test_nlmixr.r
@@ -1,99 +1,99 @@
 
 
-dmta_ds <- lapply(1:8, function(i) {
-  ds_i <- dimethenamid_2018$ds[[i]]$data
-  ds_i[ds_i$name == "DMTAP", "name"] <-  "DMTA"
-  ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i]
-  ds_i
-})
-names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title)
-dmta_ds[["Borstel"]] <- rbind(dmta_ds[["Borstel 1"]], dmta_ds[["Borstel 2"]])
-dmta_ds[["Borstel 1"]] <- NULL
-dmta_ds[["Borstel 2"]] <- NULL
-dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]])
-dmta_ds[["Elliot 1"]] <- NULL
-dmta_ds[["Elliot 2"]] <- NULL
-dfop_sfo3_plus <- mkinmod(
-  DMTA = mkinsub("DFOP", c("M23", "M27", "M31")),
-  M23 = mkinsub("SFO"),
-  M27 = mkinsub("SFO"),
-  M31 = mkinsub("SFO", "M27", sink = FALSE),
-  quiet = TRUE
-)
-f_dmta_mkin_tc <- mmkin(
-  list("DFOP-SFO3+" = dfop_sfo3_plus),
-  dmta_ds, quiet = TRUE, error_model = "tc")
-
-d_dmta_nlmixr <- nlmixr_data(f_dmta_mkin_tc)
-m_dmta_nlmixr <- function ()
-{
-    ini({
-        DMTA_0 = 98.7697627680706
-        eta.DMTA_0 ~ 2.35171765917765
-        log_k_M23 = -3.92162409637283
-        eta.log_k_M23 ~ 0.549278519419884
-        log_k_M27 = -4.33774620773911
-        eta.log_k_M27 ~ 0.864474956685295
-        log_k_M31 = -4.24767627688461
-        eta.log_k_M31 ~ 0.750297149164171
-        f_DMTA_tffm0_1_qlogis = -2.092409
-        eta.f_DMTA_tffm0_1_qlogis ~ 0.3
-        f_DMTA_tffm0_2_qlogis = -2.180576
-        eta.f_DMTA_tffm0_2_qlogis ~ 0.3
-        f_DMTA_tffm0_3_qlogis = -2.142672
-        eta.f_DMTA_tffm0_3_qlogis ~ 0.3
-        log_k1 = -2.2341008812259
-        eta.log_k1 ~ 0.902976221565793
-        log_k2 = -3.7762779983269
-        eta.log_k2 ~ 1.57684519529298
-        g_qlogis = 0.450175725479389
-        eta.g_qlogis ~ 3.0851335687675
-        sigma_low_DMTA = 0.697933852349996
-        rsd_high_DMTA = 0.0257724286053519
-        sigma_low_M23 = 0.697933852349996
-        rsd_high_M23 = 0.0257724286053519
-        sigma_low_M27 = 0.697933852349996
-        rsd_high_M27 = 0.0257724286053519
-        sigma_low_M31 = 0.697933852349996
-        rsd_high_M31 = 0.0257724286053519
-    })
-    model({
-        DMTA_0_model = DMTA_0 + eta.DMTA_0
-        DMTA(0) = DMTA_0_model
-        k_M23 = exp(log_k_M23 + eta.log_k_M23)
-        k_M27 = exp(log_k_M27 + eta.log_k_M27)
-        k_M31 = exp(log_k_M31 + eta.log_k_M31)
-        k1 = exp(log_k1 + eta.log_k1)
-        k2 = exp(log_k2 + eta.log_k2)
-        g = expit(g_qlogis + eta.g_qlogis)
-        f_DMTA_tffm0_1 = expit(f_DMTA_tffm0_1_qlogis + eta.f_DMTA_tffm0_1_qlogis)
-        f_DMTA_tffm0_2 = expit(f_DMTA_tffm0_2_qlogis + eta.f_DMTA_tffm0_2_qlogis)
-        f_DMTA_tffm0_3 = expit(f_DMTA_tffm0_3_qlogis + eta.f_DMTA_tffm0_3_qlogis)
-        f_DMTA_to_M23 = f_DMTA_tffm0_1
-        f_DMTA_to_M27 = (1 - f_DMTA_tffm0_1) * f_DMTA_tffm0_2
-        f_DMTA_to_M31 = (1 - f_DMTA_tffm0_1) * (1 - f_DMTA_tffm0_2) * f_DMTA_tffm0_3
-        d/dt(DMTA) = -((k1 * g * exp(-k1 * time) + k2 * (1 -
-            g) * exp(-k2 * time))/(g * exp(-k1 * time) + (1 -
-            g) * exp(-k2 * time))) * DMTA
-        d/dt(M23) = +f_DMTA_to_M23 * ((k1 * g * exp(-k1 * time) +
-            k2 * (1 - g) * exp(-k2 * time))/(g * exp(-k1 * time) +
-            (1 - g) * exp(-k2 * time))) * DMTA - k_M23 * M23
-        d/dt(M27) = +f_DMTA_to_M27 * ((k1 * g * exp(-k1 * time) +
-            k2 * (1 - g) * exp(-k2 * time))/(g * exp(-k1 * time) +
-            (1 - g) * exp(-k2 * time))) * DMTA - k_M27 * M27 +
-            k_M31 * M31
-        d/dt(M31) = +f_DMTA_to_M31 * ((k1 * g * exp(-k1 * time) +
-            k2 * (1 - g) * exp(-k2 * time))/(g * exp(-k1 * time) +
-            (1 - g) * exp(-k2 * time))) * DMTA - k_M31 * M31
-        DMTA ~ add(sigma_low_DMTA) + prop(rsd_high_DMTA)
-        M23 ~ add(sigma_low_M23) + prop(rsd_high_M23)
-        M27 ~ add(sigma_low_M27) + prop(rsd_high_M27)
-        M31 ~ add(sigma_low_M31) + prop(rsd_high_M31)
-    })
-}
-m_dmta_nlmixr_mkin <- nlmixr_model(f_dmta_mkin_tc, test_log_parms = TRUE)
-f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem", control = saemControl(print = 250))
-f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei", control = foceiControl(print = 250))
-plot(f_dmta_nlmixr_saem)
-plot(f_dmta_nlmixr_focei)
-
+# dmta_ds <- lapply(1:8, function(i) {
+#   ds_i <- dimethenamid_2018$ds[[i]]$data
+#   ds_i[ds_i$name == "DMTAP", "name"] <-  "DMTA"
+#   ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i]
+#   ds_i
+# })
+# names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title)
+# dmta_ds[["Borstel"]] <- rbind(dmta_ds[["Borstel 1"]], dmta_ds[["Borstel 2"]])
+# dmta_ds[["Borstel 1"]] <- NULL
+# dmta_ds[["Borstel 2"]] <- NULL
+# dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]])
+# dmta_ds[["Elliot 1"]] <- NULL
+# dmta_ds[["Elliot 2"]] <- NULL
+# dfop_sfo3_plus <- mkinmod(
+#   DMTA = mkinsub("DFOP", c("M23", "M27", "M31")),
+#   M23 = mkinsub("SFO"),
+#   M27 = mkinsub("SFO"),
+#   M31 = mkinsub("SFO", "M27", sink = FALSE),
+#   quiet = TRUE
+# )
+# f_dmta_mkin_tc <- mmkin(
+#   list("DFOP-SFO3+" = dfop_sfo3_plus),
+#   dmta_ds, quiet = TRUE, error_model = "tc")
+# 
+# d_dmta_nlmixr <- nlmixr_data(f_dmta_mkin_tc)
+# m_dmta_nlmixr <- function ()
+# {
+#     ini({
+#         DMTA_0 = 98.7697627680706
+#         eta.DMTA_0 ~ 2.35171765917765
+#         log_k_M23 = -3.92162409637283
+#         eta.log_k_M23 ~ 0.549278519419884
+#         log_k_M27 = -4.33774620773911
+#         eta.log_k_M27 ~ 0.864474956685295
+#         log_k_M31 = -4.24767627688461
+#         eta.log_k_M31 ~ 0.750297149164171
+#         f_DMTA_tffm0_1_qlogis = -2.092409
+#         eta.f_DMTA_tffm0_1_qlogis ~ 0.3
+#         f_DMTA_tffm0_2_qlogis = -2.180576
+#         eta.f_DMTA_tffm0_2_qlogis ~ 0.3
+#         f_DMTA_tffm0_3_qlogis = -2.142672
+#         eta.f_DMTA_tffm0_3_qlogis ~ 0.3
+#         log_k1 = -2.2341008812259
+#         eta.log_k1 ~ 0.902976221565793
+#         log_k2 = -3.7762779983269
+#         eta.log_k2 ~ 1.57684519529298
+#         g_qlogis = 0.450175725479389
+#         eta.g_qlogis ~ 3.0851335687675
+#         sigma_low_DMTA = 0.697933852349996
+#         rsd_high_DMTA = 0.0257724286053519
+#         sigma_low_M23 = 0.697933852349996
+#         rsd_high_M23 = 0.0257724286053519
+#         sigma_low_M27 = 0.697933852349996
+#         rsd_high_M27 = 0.0257724286053519
+#         sigma_low_M31 = 0.697933852349996
+#         rsd_high_M31 = 0.0257724286053519
+#     })
+#     model({
+#         DMTA_0_model = DMTA_0 + eta.DMTA_0
+#         DMTA(0) = DMTA_0_model
+#         k_M23 = exp(log_k_M23 + eta.log_k_M23)
+#         k_M27 = exp(log_k_M27 + eta.log_k_M27)
+#         k_M31 = exp(log_k_M31 + eta.log_k_M31)
+#         k1 = exp(log_k1 + eta.log_k1)
+#         k2 = exp(log_k2 + eta.log_k2)
+#         g = expit(g_qlogis + eta.g_qlogis)
+#         f_DMTA_tffm0_1 = expit(f_DMTA_tffm0_1_qlogis + eta.f_DMTA_tffm0_1_qlogis)
+#         f_DMTA_tffm0_2 = expit(f_DMTA_tffm0_2_qlogis + eta.f_DMTA_tffm0_2_qlogis)
+#         f_DMTA_tffm0_3 = expit(f_DMTA_tffm0_3_qlogis + eta.f_DMTA_tffm0_3_qlogis)
+#         f_DMTA_to_M23 = f_DMTA_tffm0_1
+#         f_DMTA_to_M27 = (1 - f_DMTA_tffm0_1) * f_DMTA_tffm0_2
+#         f_DMTA_to_M31 = (1 - f_DMTA_tffm0_1) * (1 - f_DMTA_tffm0_2) * f_DMTA_tffm0_3
+#         d/dt(DMTA) = -((k1 * g * exp(-k1 * time) + k2 * (1 -
+#             g) * exp(-k2 * time))/(g * exp(-k1 * time) + (1 -
+#             g) * exp(-k2 * time))) * DMTA
+#         d/dt(M23) = +f_DMTA_to_M23 * ((k1 * g * exp(-k1 * time) +
+#             k2 * (1 - g) * exp(-k2 * time))/(g * exp(-k1 * time) +
+#             (1 - g) * exp(-k2 * time))) * DMTA - k_M23 * M23
+#         d/dt(M27) = +f_DMTA_to_M27 * ((k1 * g * exp(-k1 * time) +
+#             k2 * (1 - g) * exp(-k2 * time))/(g * exp(-k1 * time) +
+#             (1 - g) * exp(-k2 * time))) * DMTA - k_M27 * M27 +
+#             k_M31 * M31
+#         d/dt(M31) = +f_DMTA_to_M31 * ((k1 * g * exp(-k1 * time) +
+#             k2 * (1 - g) * exp(-k2 * time))/(g * exp(-k1 * time) +
+#             (1 - g) * exp(-k2 * time))) * DMTA - k_M31 * M31
+#         DMTA ~ add(sigma_low_DMTA) + prop(rsd_high_DMTA)
+#         M23 ~ add(sigma_low_M23) + prop(rsd_high_M23)
+#         M27 ~ add(sigma_low_M27) + prop(rsd_high_M27)
+#         M31 ~ add(sigma_low_M31) + prop(rsd_high_M31)
+#     })
+# }
+# m_dmta_nlmixr_mkin <- nlmixr_model(f_dmta_mkin_tc, test_log_parms = TRUE)
+# f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem", control = saemControl(print = 250))
+# f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei", control = foceiControl(print = 250))
+# plot(f_dmta_nlmixr_saem)
+# plot(f_dmta_nlmixr_focei)
+#