Fix nlmixr fits with parallel metabolites

4 files changed, 80 insertions, 17 deletions

diff --git a/R/nlmixr.R b/R/nlmixr.R
index fd12f555..e41777ac 100644
--- a/R/nlmixr.R
+++ b/R/nlmixr.R
@@ -157,6 +157,30 @@ nlmixr::nlmixr
 #' plot(f_nlmixr_fomc_sfo_focei_obs_tc)
 #' summary(f_nlmixr_fomc_sfo_focei_obs_tc)
 #' }
+#'
+#' # Two parallel metabolites
+#' dmta_ds <- lapply(1:7, 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[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]])
+#' dmta_ds[["Elliot 1"]] <- NULL
+#' dmta_ds[["Elliot 2"]] <- NULL
+#' sfo_sfo2 <- mkinmod(
+#'   DMTA = mkinsub("SFO", c("M23", "M27")),
+#'   M23 = mkinsub("SFO"),
+#'   M27 = mkinsub("SFO"),
+#'   quiet = TRUE
+#' )
+#' f_dmta_sfo_sfo2 <- mmkin(
+#'   list("SFO-SFO2" = sfo_sfo2),
+#'   dmta_ds, quiet = TRUE, error_model = "obs")
+#' nlmixr_model(f_dmta_sfo_sfo2)
+#' nlmixr_focei_dmta_sfo_sfo2 <- nlmixr(f_dmta_sfo_sfo2, est = "focei")
+#'
 #' @export
 nlmixr.mmkin <- function(object, data = NULL,
   est = NULL, control = list(),
@@ -452,11 +476,15 @@ nlmixr_model <- function(object,
 
   # Population initial values for logit transformed parameters
   for (parm_name in grep("_qlogis$", names(degparms_start), value = TRUE)) {
-    parm_name_new <- names(
+    if (grepl("_tffm0_", parm_name)) {
+      parm_name_new <- gsub("_qlogis$", "", parm_name)
+    } else {
+      parm_name_new <- names(
       backtransform_odeparms(degparms_start[parm_name],
         object[[1]]$mkinmod,
         object[[1]]$transform_rates,
         object[[1]]$transform_fractions))
+    }
     model_block <- paste0(
       model_block,
       parm_name_new, " = ",
diff --git a/man/nlmixr.mmkin.Rd b/man/nlmixr.mmkin.Rd
index f9349727..e8a9d170 100644
--- a/man/nlmixr.mmkin.Rd
+++ b/man/nlmixr.mmkin.Rd
@@ -216,6 +216,30 @@ AIC(
 plot(f_nlmixr_fomc_sfo_focei_obs_tc)
 summary(f_nlmixr_fomc_sfo_focei_obs_tc)
 }
+
+# Two parallel metabolites
+dmta_ds <- lapply(1:7, 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[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]])
+dmta_ds[["Elliot 1"]] <- NULL
+dmta_ds[["Elliot 2"]] <- NULL
+sfo_sfo2 <- mkinmod(
+  DMTA = mkinsub("SFO", c("M23", "M27")),
+  M23 = mkinsub("SFO"),
+  M27 = mkinsub("SFO"),
+  quiet = TRUE
+)
+f_dmta_sfo_sfo2 <- mmkin(
+  list("SFO-SFO2" = sfo_sfo2),
+  dmta_ds, quiet = TRUE, error_model = "obs")
+nlmixr_model(f_dmta_sfo_sfo2)
+nlmixr_focei_dmta_sfo_sfo2 <- nlmixr(f_dmta_sfo_sfo2, est = "focei")
+
 }
 \seealso{
 \link{summary.nlmixr.mmkin} \link{plot.mixed.mmkin}
diff --git a/test.log b/test.log
index cb7edba3..6dea88fb 100644
--- a/test.log
+++ b/test.log
@@ -3,20 +3,20 @@ Loading required package: parallel
 ℹ Testing mkin
 ✔ | F W S  OK | Context
 ✔ |         5 | AIC calculation
-✔ |         5 | Analytical solutions for coupled models [4.4s]
+✔ |         5 | Analytical solutions for coupled models [4.3s]
 ✔ |         5 | Calculation of Akaike weights
 ✔ |         2 | Export dataset for reading into CAKE
-✔ |        12 | Confidence intervals and p-values [1.0s]
+✔ |        12 | Confidence intervals and p-values [1.1s]
 ⠋ |         1 | Dimethenamid data from 2018                                     
-✔ |     1  25 | Dimethenamid data from 2018 [52.3s]
+✔ |     1  27 | Dimethenamid data from 2018 [62.8s]
 ────────────────────────────────────────────────────────────────────────────────
-Skip (test_dmta.R:147:3): Different backends get consistent results for SFO-SFO3+, dimethenamid data
+Skip (test_dmta.R:160:3): Different backends get consistent results for SFO-SFO3+, dimethenamid data
 Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
 ────────────────────────────────────────────────────────────────────────────────
-✔ |        14 | Error model fitting [7.1s]
+✔ |        14 | Error model fitting [7.0s]
 ✔ |         5 | Time step normalisation
 ✔ |         4 | Calculation of FOCUS chi2 error levels [0.8s]
-✔ |        14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [1.2s]
+✔ |        14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [1.1s]
 ✔ |         4 | Test fitting the decline of metabolites from their maximum [0.6s]
 ✔ |         1 | Fitting the logistic model [0.4s]
 ⠴ |         6 | Nonlinear mixed-effects models                                  
@@ -30,25 +30,25 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve
 ✔ |         3 | mkinfit features [1.1s]
 ✔ |         8 | mkinmod model generation and printing [0.2s]
 ✔ |         3 | Model predictions with mkinpredict [0.3s]
-✔ |        16 | Evaluations according to 2015 NAFTA guidance [2.2s]
+✔ |        16 | Evaluations according to 2015 NAFTA guidance [2.1s]
 ✔ |         9 | Nonlinear mixed-effects models with nlme [9.2s]
 ✔ |        16 | Plotting [1.5s]
 ✔ |         4 | Residuals extracted from mkinfit models
-✔ |        23 | saemix_parent [29.1s]
-✔ |         2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.6s]
-✔ |         7 | Fitting the SFORB model [4.4s]
+✔ |        23 | saemix_parent [29.4s]
+✔ |         2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.7s]
+✔ |         7 | Fitting the SFORB model [4.5s]
 ✔ |         1 | Summaries of old mkinfit objects
 ✔ |         4 | Summary [0.1s]
 ✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.6s]
-✔ |         9 | Hypothesis tests [9.4s]
+✔ |         9 | Hypothesis tests [9.5s]
 ✔ |         2 | tffm0
 ✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 135.1 s
+Duration: 145.9 s
 
 ── Skipped tests  ──────────────────────────────────────────────────────────────
 • Fitting this ODE model with saemix takes about 15 minutes on my system (1)
 • Fitting with saemix takes around 10 minutes when using deSolve (1)
 
-[ FAIL 0 | WARN 0 | SKIP 2 | PASS 239 ]
+[ FAIL 0 | WARN 0 | SKIP 2 | PASS 241 ]
diff --git a/tests/testthat/test_dmta.R b/tests/testthat/test_dmta.R
index 7b130999..71e68f51 100644
--- a/tests/testthat/test_dmta.R
+++ b/tests/testthat/test_dmta.R
@@ -144,6 +144,19 @@ test_that("Different backends get consistent results for SFO-SFO3+, dimethenamid
     "Iteration 5, LME step.*not converge")
   ints_nlme_mets <- intervals(nlme_sfo_sfo3p_tc, which = "fixed")
 
+  # The saem fit with nlmixr takes only about 15 seconds
+  tmp <- capture.output(
+    saem_nlmixr_sfo_sfo3p_tc <- nlmixr(dmta_sfo_sfo3p_tc, est = "saem",
+      control = nlmixr::saemControl(print = 0)))
+  ints_nlmixr_saem_mets <- intervals(saem_nlmixr_sfo_sfo3p_tc)
+
+  # We need to exclude the ilr transformed formation fractions in these
+  # tests, as they do not have a one to one relation in the transformations
+  expect_true(all(ints_nlmixr_saem_mets$fixed[, "est."][-c(6, 7, 8)] >
+      backtransform_odeparms(ints_nlme_mets$fixed[, "lower"][-c(6, 7, 8)], sfo_sfo3p)))
+  expect_true(all(ints_nlmixr_saem_mets$fixed[, "est."][-c(6, 7, 8)] <
+      backtransform_odeparms(ints_nlme_mets$fixed[, "upper"], sfo_sfo3p)[-c(6, 7, 8)]))
+
   skip("Fitting this ODE model with saemix takes about 15 minutes on my system")
   # As DFOP is overparameterised and leads to instabilities and errors, we
   # need to use SFO.
@@ -156,15 +169,13 @@ test_that("Different backends get consistent results for SFO-SFO3+, dimethenamid
     fixed.estim = c(0, rep(1, 7)))
   ints_saemix_mets <- intervals(saem_saemix_sfo_sfo3p_tc_DMTA_0_fixed)
 
-  # We need to exclude the ilr transformed formation fractions in these
+  # Again, we need to exclude the ilr transformed formation fractions in these
   # tests, as they do not have a one to one relation in the transformations
   expect_true(all(ints_saemix_mets$fixed[, "est."][-c(6, 7, 8)] >
       backtransform_odeparms(ints_nlme_mets$fixed[, "lower"][-c(6, 7, 8)], sfo_sfo3p)))
   expect_true(all(ints_saemix_mets$fixed[, "est."][-c(6, 7, 8)] <
       backtransform_odeparms(ints_nlme_mets$fixed[, "upper"], sfo_sfo3p)[-c(6, 7, 8)]))
 
-  # The model is not supported by nlmixr.mmkin
-  #saem_nlmixr_sfo_sfo3p_tc <- nlmixr(dmta_sfo_sfo3p_tc, est = "saem")
 
 })