From 21ad91256dc29423bd905de5c298fd23862b1f3b Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Mon, 14 Nov 2022 20:03:42 +0100 Subject: Automatic starting parameters for saem.mmkin For the case of mkin transformations. This gives faster convergence, and appears to avoid problems with numeric ODE solutions --- .../multistart/llhist-for-biphasic-saemix-fit.svg | 41 +- .../multistart/parplot-for-biphasic-saemix-fit.svg | 242 +- ...t-for-saem-object-with-mkin-transformations.svg | 2332 ++++++++++---------- tests/testthat/anova_sfo_saem.txt | 2 +- tests/testthat/print_dfop_saemix_1.txt | 16 +- tests/testthat/summary_hfit_sfo_tc.txt | 6 +- tests/testthat/test_multistart.R | 5 +- tests/testthat/test_saemix_parent.R | 40 +- 8 files changed, 1343 insertions(+), 1341 deletions(-) (limited to 'tests') diff --git a/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg b/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg index fe38865d..6015aed8 100644 --- a/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg +++ b/tests/testthat/_snaps/multistart/llhist-for-biphasic-saemix-fit.svg @@ -21,24 +21,28 @@ Frequency of log likelihoods - - - + + + + --1185 --1180 --1175 --1170 --1165 +-1149.5 +-1149.4 +-1149.3 +-1149.2 +-1149.1 +-1149.0 - - + + + 0 -1 -2 -3 +1 +2 +3 +4 @@ -46,11 +50,12 @@ - - - - - + + + + + + original fit diff --git a/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg b/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg index 75519b07..c0332fd5 100644 --- a/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg +++ b/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg @@ -25,105 +25,109 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -163,26 +167,26 @@ - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + diff --git a/tests/testthat/_snaps/plot/mixed-model-fit-for-saem-object-with-mkin-transformations.svg b/tests/testthat/_snaps/plot/mixed-model-fit-for-saem-object-with-mkin-transformations.svg index 375ab089..69fa6a4d 100644 --- a/tests/testthat/_snaps/plot/mixed-model-fit-for-saem-object-with-mkin-transformations.svg +++ b/tests/testthat/_snaps/plot/mixed-model-fit-for-saem-object-with-mkin-transformations.svg @@ -156,7 +156,7 @@ - + @@ -175,7 +175,7 @@ - + @@ -212,7 +212,7 @@ - + @@ -249,7 +249,7 @@ - + @@ -268,7 +268,7 @@ - + @@ -287,7 +287,7 @@ - + @@ -342,7 +342,7 @@ - + @@ -415,7 +415,7 @@ - + @@ -470,7 +470,7 @@ - + @@ -525,7 +525,7 @@ - + @@ -562,7 +562,7 @@ - + @@ -617,7 +617,7 @@ - + @@ -672,7 +672,7 @@ - + @@ -709,7 +709,7 @@ - + @@ -728,8 +728,8 @@ - - + + @@ -739,34 +739,34 @@ - + - - - - - + + + + + 0 -20 -40 -60 -80 -100 - - - - +20 +40 +60 +80 +100 + + + + - - - --3 --2 --1 + + + +-3 +-2 +-1 0 -1 -2 -3 +1 +2 +3 @@ -780,582 +780,582 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -1419,7 +1419,7 @@ - + @@ -1436,7 +1436,7 @@ - + @@ -1469,7 +1469,7 @@ - + @@ -1502,7 +1502,7 @@ - + @@ -1519,7 +1519,7 @@ - + @@ -1536,7 +1536,7 @@ - + @@ -1585,7 +1585,7 @@ - + @@ -1650,7 +1650,7 @@ - + @@ -1699,7 +1699,7 @@ - + @@ -1748,7 +1748,7 @@ - + @@ -1781,7 +1781,7 @@ - + @@ -1830,7 +1830,7 @@ - + @@ -1879,7 +1879,7 @@ - + @@ -1912,7 +1912,7 @@ - + @@ -1929,8 +1929,8 @@ - - + + @@ -1940,32 +1940,32 @@ - + - - - - + + + + 0 -10 -20 -30 -40 - - - - +10 +20 +30 +40 + + + + - - - --3 --2 --1 + + + +-3 +-2 +-1 0 -1 -2 -3 +1 +2 +3 @@ -1979,518 +1979,518 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/testthat/anova_sfo_saem.txt b/tests/testthat/anova_sfo_saem.txt index 4a21e81f..9e4bf71f 100644 --- a/tests/testthat/anova_sfo_saem.txt +++ b/tests/testthat/anova_sfo_saem.txt @@ -4,4 +4,4 @@ Data: 262 observations of 1 variable(s) grouped in 15 datasets sfo_saem_1_reduced 5 1310 1313 -650 sfo_saem_1_reduced_mkin 5 1310 1313 -650 0 0 sfo_saem_1 6 1312 1316 -650 0 1 1 -sfo_saem_1_mkin 6 1310 1315 -649 1 0 +sfo_saem_1_mkin 6 1312 1316 -650 0 0 diff --git a/tests/testthat/print_dfop_saemix_1.txt b/tests/testthat/print_dfop_saemix_1.txt index f427b3e6..f6fda37c 100644 --- a/tests/testthat/print_dfop_saemix_1.txt +++ b/tests/testthat/print_dfop_saemix_1.txt @@ -13,12 +13,12 @@ Likelihood computed by importance sampling Fitted parameters: estimate lower upper -parent_0 100.04 98.89 101.20 -log_k1 -4.12 -4.24 -4.00 -log_k2 -2.67 -2.90 -2.44 -g_qlogis 0.43 0.22 0.64 -a.1 0.92 0.67 1.16 +parent_0 100.09 98.94 101.25 +log_k1 -2.68 -2.91 -2.45 +log_k2 -4.12 -4.24 -4.00 +g_qlogis -0.41 -0.63 -0.20 +a.1 0.91 0.67 1.15 b.1 0.05 0.04 0.06 -SD.log_k1 0.22 0.14 0.30 -SD.log_k2 0.36 0.21 0.51 -SD.g_qlogis 0.14 -0.11 0.39 +SD.log_k1 0.36 0.21 0.50 +SD.log_k2 0.22 0.13 0.30 +SD.g_qlogis 0.15 -0.09 0.40 diff --git a/tests/testthat/summary_hfit_sfo_tc.txt b/tests/testthat/summary_hfit_sfo_tc.txt index 49765187..41743091 100644 --- a/tests/testthat/summary_hfit_sfo_tc.txt +++ b/tests/testthat/summary_hfit_sfo_tc.txt @@ -32,9 +32,9 @@ Likelihood computed by importance sampling Optimised parameters: est. lower upper -parent_0 101.01 99.57 102.45 +parent_0 101.02 99.58 102.46 log_k_parent -3.32 -3.53 -3.11 -a.1 0.90 0.64 1.17 +a.1 0.91 0.64 1.17 b.1 0.05 0.04 0.06 SD.log_k_parent 0.27 0.11 0.42 @@ -48,7 +48,7 @@ SD.log_k_parent 0.3 0.1 0.4 Variance model: est. lower upper -a.1 0.90 0.64 1.17 +a.1 0.91 0.64 1.17 b.1 0.05 0.04 0.06 Backtransformed parameters: diff --git a/tests/testthat/test_multistart.R b/tests/testthat/test_multistart.R index 56eb140c..c1a10d10 100644 --- a/tests/testthat/test_multistart.R +++ b/tests/testthat/test_multistart.R @@ -21,7 +21,10 @@ test_that("multistart works for saem.mmkin models", { anova_biphasic <- anova(saem_biphasic_m, best(saem_biphasic_m_multi)) - expect_true(anova_biphasic[2, "AIC"] < anova_biphasic[1, "AIC"]) + # With the new starting parameters we do not improve + # with multistart any more + expect_equal(anova_biphasic[2, "AIC"], anova_biphasic[1, "AIC"], + tolerance = 1e-4) skip_on_travis() # Plots are platform dependent llhist_sfo <- function() llhist(saem_sfo_s_multi) diff --git a/tests/testthat/test_saemix_parent.R b/tests/testthat/test_saemix_parent.R index 79c5fa69..20889c6c 100644 --- a/tests/testthat/test_saemix_parent.R +++ b/tests/testthat/test_saemix_parent.R @@ -73,40 +73,36 @@ test_that("Parent fits using saemix are correctly implemented", { # FOMC mmkin_fomc_1 <- mmkin("FOMC", ds_fomc, quiet = TRUE, error_model = "tc", cores = n_cores) fomc_saem_1 <- saem(mmkin_fomc_1, quiet = TRUE, transformations = "saemix", no_random_effect = "parent_0") - fomc_saem_2 <- update(fomc_saem_1, transformations = "mkin") - ci_fomc_s1 <- summary(fomc_saem_1)$confint_back fomc_pop <- as.numeric(fomc_pop) + ci_fomc_s1 <- summary(fomc_saem_1)$confint_back expect_true(all(ci_fomc_s1[, "lower"] < fomc_pop)) expect_true(all(ci_fomc_s1[, "upper"] > fomc_pop)) - expect_equal(endpoints(fomc_saem_1), endpoints(fomc_saem_2), tol = 0.01) - mmkin_fomc_2 <- update(mmkin_fomc_1, state.ini = 100, fixed_initials = "parent") - fomc_saem_2 <- saem(mmkin_fomc_2, quiet = TRUE, transformations = "mkin") + fomc_saem_2 <- update(fomc_saem_1, transformations = "mkin") ci_fomc_s2 <- summary(fomc_saem_2)$confint_back + expect_true(all(ci_fomc_s2[, "lower"] < fomc_pop)) + expect_true(all(ci_fomc_s2[, "upper"] > fomc_pop)) - expect_true(all(ci_fomc_s2[, "lower"] < fomc_pop[2:3])) - expect_true(all(ci_fomc_s2[, "upper"] > fomc_pop[2:3])) + expect_equal(endpoints(fomc_saem_1), endpoints(fomc_saem_2), tol = 0.01) # DFOP dfop_saemix_2 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "saemix", no_random_effect = "parent_0") - s_dfop_s1 <- summary(dfop_saemix_1) - s_dfop_s2 <- summary(dfop_saemix_2) + s_dfop_s1 <- summary(dfop_saemix_1) # mkin transformations + s_dfop_s2 <- summary(dfop_saemix_2) # saemix transformations s_dfop_n <- summary(dfop_nlme_1) dfop_pop <- as.numeric(dfop_pop) - # When using DFOP with mkin transformations, k1 and k2 are sometimes swapped - swap_k1_k2 <- function(p) c(p[1], p[3], p[2], 1 - p[4]) - expect_true(all(s_dfop_s1$confint_back[, "lower"] < swap_k1_k2(dfop_pop))) - expect_true(all(s_dfop_s1$confint_back[, "upper"] > swap_k1_k2(dfop_pop))) + expect_true(all(s_dfop_s1$confint_back[, "lower"] < dfop_pop)) + expect_true(all(s_dfop_s1$confint_back[, "upper"] > dfop_pop)) expect_true(all(s_dfop_s2$confint_back[, "lower"] < dfop_pop)) expect_true(all(s_dfop_s2$confint_back[, "upper"] > dfop_pop)) - # We get < 20% deviations with transformations made in mkin (need to swap k1 and k2) - rel_diff_1 <- (swap_k1_k2(s_dfop_s1$confint_back[, "est."]) - dfop_pop) / dfop_pop + # We get < 20% deviations with transformations made in mkin + rel_diff_1 <- (s_dfop_s1$confint_back[, "est."] - dfop_pop) / dfop_pop expect_true(all(rel_diff_1 < 0.20)) # We get < 20% deviations with transformations made in saemix @@ -123,27 +119,21 @@ test_that("Parent fits using saemix are correctly implemented", { transformations = "saemix") expect_equal( log(endpoints(dfop_saemix_1)$distimes[1:2]), - log(endpoints(sforb_saemix_1)$distimes[1:2]), tolerance = 0.03) + log(endpoints(sforb_saemix_1)$distimes[1:2]), tolerance = 0.01) expect_equal( log(endpoints(sforb_saemix_1)$distimes[1:2]), log(endpoints(sforb_saemix_2)$distimes[1:2]), tolerance = 0.01) 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") + hs_saem_1 <- saem(mmkin_hs_1, quiet = TRUE, no_random_effect = "parent_0") + hs_saem_2 <- saem(mmkin_hs_1, quiet = TRUE, transformations = "mkin", + no_random_effect = "parent_0") 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[, "upper"] > hs_pop)) - - mmkin_hs_2 <- update(mmkin_hs_1, state.ini = 100, fixed_initials = "parent") - hs_saem_3 <- saem(mmkin_hs_2, quiet = TRUE) - ci_hs_s3 <- summary(hs_saem_3)$confint_back - - #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