From 22d21cf5efcfb52c59c969d393bb0be077e982dd Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Mon, 8 Aug 2022 10:11:04 +0200 Subject: Fix fitting HS with saemix transformations --- NEWS.md | 2 ++ R/saem.R | 27 +++++++++++++++++++-------- log/test.log | 22 ++++++++++------------ tests/testthat/test_saemix_parent.R | 15 ++++++++------- 4 files changed, 39 insertions(+), 27 deletions(-) diff --git a/NEWS.md b/NEWS.md index 55462fd3..64700704 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # mkin 1.1.2 +- 'R/saem.R': Implement and test saemix transformations for HS. + - 'R/convergence.R': New generic to show convergence information with a method for 'mmin' objects. - 'R/illparms.R': New generic to show ill-defined parameters with methods for 'mkinfit' and 'mmkin' objects. diff --git a/R/saem.R b/R/saem.R index 36997ad7..d2240ea5 100644 --- a/R/saem.R +++ b/R/saem.R @@ -252,7 +252,7 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk } degparms_fixed <- object[[1]]$bparms.fixed - # Transformations are done in the degradation function + # Transformations are done in the degradation function by default transform.par = rep(0, length(degparms_optim)) odeini_optim_parm_names <- grep('_0$', names(degparms_optim), value = TRUE) @@ -339,13 +339,24 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk } } if (parent_type == "HS") { - model_function <- function(psi, id, xidep) { - tb <- exp(psi[id, 4]) - t <- xidep[, "time"] - k1 = exp(psi[id, 2]) - psi[id, 1] * ifelse(t <= tb, - exp(- k1 * t), - exp(- k1 * tb) * exp(- exp(psi[id, 3]) * (t - tb))) + if (transformations == "mkin") { + model_function <- function(psi, id, xidep) { + tb <- exp(psi[id, 4]) + t <- xidep[, "time"] + k1 <- exp(psi[id, 2]) + psi[id, 1] * ifelse(t <= tb, + exp(- k1 * t), + exp(- k1 * tb) * exp(- exp(psi[id, 3]) * (t - tb))) + } + } else { + model_function <- function(psi, id, xidep) { + tb <- exp(psi[id, 4]) + t <- xidep[, "time"] + psi[id, 1] * ifelse(t <= tb, + exp(- psi[id, 2] * t), + exp(- psi[id, 2] * tb) * exp(- psi[id, 3] * (t - tb))) + } + transform.par = c(0, 1, 1, 1) } } } diff --git a/log/test.log b/log/test.log index 266a566f..776b6f49 100644 --- a/log/test.log +++ b/log/test.log @@ -1,13 +1,11 @@ -ℹ Loading mkin -Loading required package: parallel ℹ Testing mkin ✔ | F W S OK | Context ✔ | 5 | AIC calculation -✔ | 5 | Analytical solutions for coupled models [3.2s] +✔ | 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 [31.6s] +✔ | 1 12 | Dimethenamid data from 2018 [32.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 @@ -16,7 +14,7 @@ Reason: Fitting this ODE model with saemix takes about 15 minutes on my system ✔ | 5 | Time step normalisation ✔ | 4 | Calculation of FOCUS chi2 error levels [0.6s] ✔ | 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] +✔ | 4 | Test fitting the decline of metabolites from their maximum [0.4s] ✔ | 1 | Fitting the logistic model [0.2s] ✔ | 1 12 | Nonlinear mixed-effects models [0.2s] ──────────────────────────────────────────────────────────────────────────────── @@ -28,13 +26,13 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve ✔ | 3 | mkinfit features [0.7s] ✔ | 8 | mkinmod model generation and printing [0.2s] ✔ | 3 | Model predictions with mkinpredict [0.4s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.7s] -✔ | 9 | Nonlinear mixed-effects models with nlme [8.3s] -✔ | 16 | Plotting [10.6s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.8s] +✔ | 9 | Nonlinear mixed-effects models with nlme [9.1s] +✔ | 16 | Plotting [10.7s] ✔ | 4 | Residuals extracted from mkinfit models -✔ | 25 | saemix parent models [172.9s] +✔ | 27 | saemix parent models [178.6s] ✔ | 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 ✔ | 5 | Summary [0.2s] ✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.1s] @@ -42,10 +40,10 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve ✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 255.9 s +Duration: 263.0 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 225 ] +[ FAIL 0 | WARN 0 | SKIP 2 | PASS 227 ] diff --git a/tests/testthat/test_saemix_parent.R b/tests/testthat/test_saemix_parent.R index 731228d9..7f550999 100644 --- a/tests/testthat/test_saemix_parent.R +++ b/tests/testthat/test_saemix_parent.R @@ -10,6 +10,7 @@ test_that("Parent fits using saemix are correctly implemented", { sfo_saem_2 <- saem(mmkin_sfo_1, quiet = TRUE, transformations = "mkin") sfo_saem_3 <- expect_error(saem(mmkin_sfo_2, quiet = TRUE), "at least two parameters") + expect_equal(endpoints(sfo_saem_1), endpoints(sfo_saem_2), tol = 0.01) s_sfo_s1 <- summary(sfo_saem_1) s_sfo_s2 <- summary(sfo_saem_2) @@ -76,20 +77,20 @@ test_that("Parent fits using saemix are correctly implemented", { mmkin_hs_1 <- mmkin("HS", ds_hs, quiet = TRUE, error_model = "const", cores = n_cores) hs_saem_1 <- saem(mmkin_hs_1, quiet = TRUE) + hs_saem_2 <- saem(mmkin_hs_1, quiet = TRUE, transformations = "mkin") + expect_equal(endpoints(hs_saem_1), endpoints(hs_saem_2), tol = 0.01) ci_hs_s1 <- summary(hs_saem_1)$confint_back hs_pop <- as.numeric(hs_pop) - # expect_true(all(ci_hs_s1[, "lower"] < hs_pop)) # k1 is overestimated + #expect_true(all(ci_hs_s1[, "lower"] < hs_pop)) # k1 is overestimated expect_true(all(ci_hs_s1[, "upper"] > hs_pop)) mmkin_hs_2 <- update(mmkin_hs_1, state.ini = 100, fixed_initials = "parent") - hs_saem_2 <- saem(mmkin_hs_2, quiet = TRUE) - ci_hs_s2 <- summary(hs_saem_2)$confint_back + hs_saem_3 <- saem(mmkin_hs_2, quiet = TRUE) + ci_hs_s3 <- summary(hs_saem_3)$confint_back - #expect_true(all(ci_hs_s2[, "lower"] < hs_pop[2:4])) # k1 again overestimated - expect_true(all(ci_hs_s2[, "upper"] > hs_pop[2:4])) - - # HS would likely benefit from implemenation of transformations = "saemix" + #expect_true(all(ci_hs_s3[, "lower"] < hs_pop[2:4])) # k1 again overestimated + expect_true(all(ci_hs_s3[, "upper"] > hs_pop[2:4])) }) test_that("We can also use mkin solution methods for saem", { -- cgit v1.2.1