From 351248d07f810ccb6c497633a02cd48ee35526e6 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 2 Mar 2022 13:11:49 +0100 Subject: Fix nlmixr fits with parallel metabolites --- R/nlmixr.R | 30 +++++++++++++++++++++++++++++- man/nlmixr.mmkin.Rd | 24 ++++++++++++++++++++++++ test.log | 26 +++++++++++++------------- tests/testthat/test_dmta.R | 17 ++++++++++++++--- 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") }) -- cgit v1.2.1