From 2728910b96f0ec7dd7ccd97fc6c1f6677e5e352d Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Sat, 19 Mar 2022 10:42:07 +0100 Subject: Test saemix without special analytical solutions Also increase the performance a bit (from about 210 s to about 140 s in the case of DFOP with four chains and 300, 100 iterations). --- NEWS.md | 2 ++ R/saem.R | 32 ++++++++++++++++++++------------ test.log | 32 ++++++++++++++++---------------- tests/testthat/test_saemix_parent.R | 12 ++++++++++++ 4 files changed, 50 insertions(+), 28 deletions(-) diff --git a/NEWS.md b/NEWS.md index 1524e1cf..4a3b352d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,8 @@ - 'CAKE_export': Check for validity of the map argument, updates +- 'saem()': Slightly improve speed in the case that no analytical solution for saemix is implemented, activate a test of the respective code + # mkin 1.1.0 (2022-03-14) ## Mixed-effects models diff --git a/R/saem.R b/R/saem.R index d3b23861..36997ad7 100644 --- a/R/saem.R +++ b/R/saem.R @@ -239,10 +239,6 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) { saemix_model <- function(object, solution_type = "auto", transformations = c("mkin", "saemix"), degparms_start = numeric(), test_log_parms = FALSE, conf.level = 0.6, verbose = FALSE, ...) { - if (packageVersion("saemix") < "3.0") { - stop("To use the interface to saemix, you need to install a version >= 3.0\n") - } - if (nrow(object) > 1) stop("Only row objects allowed") mkin_model <- object[[1]]$mkinmod @@ -444,9 +440,22 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk solution_type = "analytical saemix" } else { + if (transformations == "saemix") { + stop("Using saemix transformations is only supported if an analytical solution is implemented for saemix") + } + if (solution_type == "auto") solution_type <- object[[1]]$solution_type + # Define some variables to avoid function calls in model function + transparms_optim_names <- names(degparms_optim) + odeini_optim_names <- gsub('_0$', '', odeini_optim_parm_names) + diff_names <- names(mkin_model$diffs) + ode_transparms_optim_names <- setdiff(transparms_optim_names, odeini_optim_parm_names) + transform_rates <- object[[1]]$transform_rates + transform_fractions <- object[[1]]$transform_fractions + + # Define the model function model_function <- function(psi, id, xidep) { uid <- unique(id) @@ -454,22 +463,21 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk res_list <- lapply(uid, function(i) { transparms_optim <- as.numeric(psi[i, ]) # psi[i, ] is a dataframe when called in saemix.predict - names(transparms_optim) <- names(degparms_optim) + names(transparms_optim) <- transparms_optim_names odeini_optim <- transparms_optim[odeini_optim_parm_names] - names(odeini_optim) <- gsub('_0$', '', odeini_optim_parm_names) + names(odeini_optim) <- odeini_optim_names - odeini <- c(odeini_optim, odeini_fixed)[names(mkin_model$diffs)] + odeini <- c(odeini_optim, odeini_fixed)[diff_names] - ode_transparms_optim_names <- setdiff(names(transparms_optim), odeini_optim_parm_names) odeparms_optim <- backtransform_odeparms(transparms_optim[ode_transparms_optim_names], mkin_model, - transform_rates = object[[1]]$transform_rates, - transform_fractions = object[[1]]$transform_fractions) + transform_rates = transform_rates, + transform_fractions = transform_fractions) odeparms <- c(odeparms_optim, odeparms_fixed) - xidep_i <- subset(xidep, id == i) + xidep_i <- xidep[which(id == i), ] - if (solution_type == "analytical") { + if (solution_type[1] == "analytical") { out_values <- mkin_model$deg_func(xidep_i, odeini, odeparms) } else { diff --git a/test.log b/test.log index 4c273c19..2fc6a43c 100644 --- a/test.log +++ b/test.log @@ -3,18 +3,18 @@ Loading required package: parallel ℹ Testing mkin ✔ | F W S OK | Context ✔ | 5 | AIC calculation -✔ | 5 | Analytical solutions for coupled models [3.4s] +✔ | 5 | Analytical solutions for coupled models [3.3s] ✔ | 5 | Calculation of Akaike weights ✔ | 3 | Export dataset for reading into CAKE ✔ | 12 | Confidence intervals and p-values [1.0s] -✔ | 1 12 | Dimethenamid data from 2018 [32.7s] +✔ | 1 12 | Dimethenamid data from 2018 [31.0s] ──────────────────────────────────────────────────────────────────────────────── Skip (test_dmta.R:98: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 [5.0s] +✔ | 14 | Error model fitting [4.6s] ✔ | 5 | Time step normalisation -✔ | 4 | Calculation of FOCUS chi2 error levels [0.6s] +✔ | 4 | Calculation of FOCUS chi2 error levels [0.5s] ✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.8s] ✔ | 4 | Test fitting the decline of metabolites from their maximum [0.3s] ✔ | 1 | Fitting the logistic model [0.2s] @@ -24,28 +24,28 @@ Skip (test_mixed.R:68:3): saemix results are reproducible for biphasic fits Reason: Fitting with saemix takes around 10 minutes when using deSolve ──────────────────────────────────────────────────────────────────────────────── ✔ | 3 | Test dataset classes mkinds and mkindsg -✔ | 10 | Special cases of mkinfit calls [0.5s] -✔ | 3 | mkinfit features [0.7s] -✔ | 8 | mkinmod model generation and printing [0.3s] -✔ | 3 | Model predictions with mkinpredict [0.4s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.6s] -✔ | 9 | Nonlinear mixed-effects models with nlme [8.6s] -✔ | 16 | Plotting [10.4s] +✔ | 10 | Special cases of mkinfit calls [0.4s] +✔ | 3 | mkinfit features [0.6s] +✔ | 8 | mkinmod model generation and printing [0.2s] +✔ | 3 | Model predictions with mkinpredict [0.3s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.4s] +✔ | 9 | Nonlinear mixed-effects models with nlme [8.0s] +✔ | 16 | Plotting [9.9s] ✔ | 4 | Residuals extracted from mkinfit models -✔ | 23 | saemix parent models [28.0s] +✔ | 25 | saemix parent models [170.8s] ✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s] -✔ | 7 | Fitting the SFORB model [3.7s] +✔ | 7 | Fitting the SFORB model [3.6s] ✔ | 1 | Summaries of old mkinfit objects ✔ | 4 | Summary [0.1s] ✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s] -✔ | 9 | Hypothesis tests [8.1s] +✔ | 9 | Hypothesis tests [8.0s] ✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 112.8 s +Duration: 251.5 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 222 ] +[ FAIL 0 | WARN 0 | SKIP 2 | PASS 224 ] diff --git a/tests/testthat/test_saemix_parent.R b/tests/testthat/test_saemix_parent.R index 340db6f1..731228d9 100644 --- a/tests/testthat/test_saemix_parent.R +++ b/tests/testthat/test_saemix_parent.R @@ -92,3 +92,15 @@ test_that("Parent fits using saemix are correctly implemented", { # HS would likely benefit from implemenation of transformations = "saemix" }) +test_that("We can also use mkin solution methods for saem", { + expect_error(saem(mmkin_dfop_1, quiet = TRUE, transformations = "saemix", solution_type = "analytical"), + "saemix transformations is only supported if an analytical solution is implemented" + ) + skip_on_cran() # This still takes almost 2.5 minutes although we do not solve ODEs + dfop_saemix_3 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "mkin", + solution_type = "analytical") + distimes_dfop <- endpoints(dfop_saemix_1)$distimes + distimes_dfop_analytical <- endpoints(dfop_saemix_3)$distimes + rel_diff <- abs(distimes_dfop_analytical - distimes_dfop) / distimes_dfop + expect_true(all(rel_diff < 0.01)) +}) -- cgit v1.2.1