From a8ff8bed72dc537fe70cf2995ea769d3f519f877 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 25 Oct 2022 21:45:30 +0200 Subject: Change DFOP mixed model data in tests, updates --- .../plot/mixed-model-fit-for-nlme-object.svg | 2404 ++++++++++---------- tests/testthat/anova_sfo_saem.txt | 0 tests/testthat/print_dfop_nlme_1.txt | 12 +- tests/testthat/setup_script.R | 33 +- tests/testthat/summary_dfop_nlme_1.txt | 36 +- tests/testthat/test_confidence.R | 14 + tests/testthat/test_saemix_parent.R | 56 +- tests/testthat/test_tests.R | 1 + 8 files changed, 1286 insertions(+), 1270 deletions(-) create mode 100644 tests/testthat/anova_sfo_saem.txt (limited to 'tests') diff --git a/tests/testthat/_snaps/plot/mixed-model-fit-for-nlme-object.svg b/tests/testthat/_snaps/plot/mixed-model-fit-for-nlme-object.svg index 051b46b1..76fed0dc 100644 --- a/tests/testthat/_snaps/plot/mixed-model-fit-for-nlme-object.svg +++ b/tests/testthat/_snaps/plot/mixed-model-fit-for-nlme-object.svg @@ -113,19 +113,19 @@ 80 100 120 - + - - - - - + + + + + 0 -20 -40 -60 -80 -100 +20 +40 +60 +80 +100 @@ -138,598 +138,598 @@ Residues parent - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -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 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/testthat/anova_sfo_saem.txt b/tests/testthat/anova_sfo_saem.txt new file mode 100644 index 00000000..e69de29b diff --git a/tests/testthat/print_dfop_nlme_1.txt b/tests/testthat/print_dfop_nlme_1.txt index 1092821d..435ff409 100644 --- a/tests/testthat/print_dfop_nlme_1.txt +++ b/tests/testthat/print_dfop_nlme_1.txt @@ -8,17 +8,23 @@ d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * Data: 270 observations of 1 variable(s) grouped in 15 datasets -Log-likelihood: -612 +Log-likelihood: -695 Fixed effects: list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1) parent_0 log_k1 log_k2 g_qlogis - 100.1 -2.7 -4.1 -0.5 + 100.1 -2.7 -4.1 -0.4 Random effects: Formula: list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1) Level: ds Structure: Diagonal parent_0 log_k1 log_k2 g_qlogis Residual -StdDev: 0.4 0.3 0.2 0.2 2 +StdDev: 2 0.3 0.2 0.2 1 +Variance function: + Structure: Constant plus proportion of variance covariate + Formula: ~fitted(.) + Parameter estimates: + const prop + 0.92275475 -0.04804649 diff --git a/tests/testthat/setup_script.R b/tests/testthat/setup_script.R index af059668..4de2674f 100644 --- a/tests/testthat/setup_script.R +++ b/tests/testthat/setup_script.R @@ -14,25 +14,12 @@ if (Sys.getenv("TRAVIS") != "") n_cores = 2 # On Windows we would need to make a cluster first if (Sys.info()["sysname"] == "Windows") n_cores = 1 -# We set up some models and fits with nls for comparisons -SFO_trans <- function(t, parent_0, log_k_parent_sink) { - parent_0 * exp(- exp(log_k_parent_sink) * t) -} -SFO_notrans <- function(t, parent_0, k_parent_sink) { - parent_0 * exp(- k_parent_sink * t) -} -f_1_nls_trans <- nls(value ~ SFO_trans(time, parent_0, log_k_parent_sink), - data = FOCUS_2006_A, - start = list(parent_0 = 100, log_k_parent_sink = log(0.1))) -f_1_nls_notrans <- nls(value ~ SFO_notrans(time, parent_0, k_parent_sink), - data = FOCUS_2006_A, - start = list(parent_0 = 100, k_parent_sink = 0.1)) - +# Very simple example fits f_1_mkin_trans <- mkinfit("SFO", FOCUS_2006_A, quiet = TRUE) f_1_mkin_notrans <- mkinfit("SFO", FOCUS_2006_A, quiet = TRUE, transform_rates = FALSE) -# mmkin object of parent fits for tests +# mmkin object of parent fits models <- c("SFO", "FOMC", "DFOP", "HS") fits <- suppressWarnings( # FOCUS A FOMC was, it seems, in testthat output mmkin(models, @@ -72,9 +59,8 @@ DFOP_par_c <- synthetic_data_for_UBA_2014[[12]]$data f_2_mkin <- mkinfit("DFOP", DFOP_par_c, quiet = TRUE) f_2_nls <- nls(value ~ SSbiexp(time, A1, lrc1, A2, lrc2), data = subset(DFOP_par_c, name == "parent")) -f_2_anova <- lm(value ~ as.factor(time), data = subset(DFOP_par_c, name == "parent")) -# Two metabolites +# mkinfit with two metabolites m_synth_SFO_lin <- mkinmod( parent = mkinsub("SFO", "M1"), M1 = mkinsub("SFO", "M2"), @@ -137,7 +123,7 @@ set.seed(123456) ds_dfop <- lapply(1:n, function(i) { ds_mean <- mkinpredict(DFOP, dfop_parms[i, ], c(parent = dfop_pop$parent_0), sampling_times) - add_err(ds_mean, const, n = 1)[[1]] + add_err(ds_mean, tc, n = 1)[[1]] }) set.seed(123456) @@ -184,20 +170,21 @@ ds_biphasic <- lapply(ds_biphasic_mean, function(ds) { # Mixed model fits mmkin_sfo_1 <- mmkin("SFO", ds_sfo, quiet = TRUE, error_model = "tc", cores = n_cores) -mmkin_dfop_1 <- mmkin("DFOP", ds_dfop, quiet = TRUE, cores = n_cores) +mmkin_dfop_1 <- mmkin("DFOP", ds_dfop, quiet = TRUE, cores = n_cores, + error_model = "tc") + mmkin_biphasic <- mmkin(list("DFOP-SFO" = DFOP_SFO), ds_biphasic, quiet = TRUE, cores = n_cores, control = list(eval.max = 500, iter.max = 400), error_model = "tc") # nlme -dfop_nlme_1 <- nlme(mmkin_dfop_1) +dfop_nlme_1 <- suppressWarnings(nlme(mmkin_dfop_1)) nlme_biphasic <- suppressWarnings(nlme(mmkin_biphasic)) # saemix sfo_saem_1 <- saem(mmkin_sfo_1, quiet = TRUE, transformations = "saemix") - -dfop_saemix_1 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "mkin") -dfop_saemix_2 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "saemix") +dfop_saemix_1 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "mkin", + no_random_effect = "parent_0") saem_biphasic_m <- saem(mmkin_biphasic, transformations = "mkin", quiet = TRUE) saem_biphasic_s <- saem(mmkin_biphasic, transformations = "saemix", quiet = TRUE) diff --git a/tests/testthat/summary_dfop_nlme_1.txt b/tests/testthat/summary_dfop_nlme_1.txt index 9da3fb1b..773be488 100644 --- a/tests/testthat/summary_dfop_nlme_1.txt +++ b/tests/testthat/summary_dfop_nlme_1.txt @@ -14,13 +14,13 @@ Data: Model predictions using solution type analytical -Fitted in test time 0 s using 3 iterations +Fitted in test time 0 s using 5 iterations -Variance model: Constant variance +Variance model: Two-component variance function Mean of starting values for individual parameters: parent_0 log_k1 log_k2 g_qlogis - 100.06 -2.68 -4.16 0.01 + 100.2 -2.6 -4.2 0.1 Fixed degradation parameter values: None @@ -28,36 +28,42 @@ None Results: AIC BIC logLik - 1242 1274 -612 + 1410 1446 -695 Optimised, transformed parameters with symmetric confidence intervals: lower est. upper -parent_0 99.6 100.1 100.6 -log_k1 -2.9 -2.7 -2.4 +parent_0 98.7 100.1 101.5 +log_k1 -2.9 -2.7 -2.5 log_k2 -4.2 -4.1 -4.0 -g_qlogis -0.7 -0.5 -0.2 +g_qlogis -0.7 -0.4 -0.2 Correlation: pr_0 lg_1 lg_2 -log_k1 0.2 +log_k1 0.3 log_k2 0.1 0.2 -g_qlogis -0.2 -0.5 -0.4 +g_qlogis -0.1 -0.5 -0.4 Random effects: Formula: list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1) Level: ds Structure: Diagonal parent_0 log_k1 log_k2 g_qlogis Residual -StdDev: 0.4 0.3 0.2 0.2 2 +StdDev: 2 0.3 0.2 0.2 1 +Variance function: + Structure: Constant plus proportion of variance covariate + Formula: ~fitted(.) + Parameter estimates: + const prop + 0.92275475 -0.04804649 Backtransformed parameters with asymmetric confidence intervals: lower est. upper -parent_0 1e+02 1e+02 1e+02 -k1 6e-02 7e-02 9e-02 -k2 1e-02 2e-02 2e-02 -g 3e-01 4e-01 5e-01 +parent_0 98.69 1e+02 1e+02 +k1 0.06 7e-02 9e-02 +k2 0.01 2e-02 2e-02 +g 0.34 4e-01 5e-01 Estimated disappearance times: DT50 DT90 DT50back DT50_k1 DT50_k2 -parent 23 111 33 10 42 +parent 23 111 33 10 43 diff --git a/tests/testthat/test_confidence.R b/tests/testthat/test_confidence.R index 36e9738d..6c645ca4 100644 --- a/tests/testthat/test_confidence.R +++ b/tests/testthat/test_confidence.R @@ -1,5 +1,19 @@ context("Confidence intervals and p-values") +# We set up some models and fits with nls for comparisons +SFO_trans <- function(t, parent_0, log_k_parent_sink) { + parent_0 * exp(- exp(log_k_parent_sink) * t) +} +SFO_notrans <- function(t, parent_0, k_parent_sink) { + parent_0 * exp(- k_parent_sink * t) +} +f_1_nls_trans <- nls(value ~ SFO_trans(time, parent_0, log_k_parent_sink), + data = FOCUS_2006_A, + start = list(parent_0 = 100, log_k_parent_sink = log(0.1))) +f_1_nls_notrans <- nls(value ~ SFO_notrans(time, parent_0, k_parent_sink), + data = FOCUS_2006_A, + start = list(parent_0 = 100, k_parent_sink = 0.1)) + test_that("The confint method 'quadratic' is consistent with the summary", { expect_equivalent( confint(fit_nw_1, parm = "parent_0", method = "quadratic"), diff --git a/tests/testthat/test_saemix_parent.R b/tests/testthat/test_saemix_parent.R index 39f69f51..4504e573 100644 --- a/tests/testthat/test_saemix_parent.R +++ b/tests/testthat/test_saemix_parent.R @@ -5,6 +5,7 @@ test_that("Parent fits using saemix are correctly implemented", { skip_on_cran() expect_error(saem(fits), "Only row objects") + # SFO # mmkin_sfo_1 was generated in the setup script # We did not introduce variance of parent_0 in the data generation # This is correctly detected @@ -53,9 +54,16 @@ test_that("Parent fits using saemix are correctly implemented", { expect_equal(round(s_sfo_nlme_1$confint_back["k_parent", "est."], 3), round(s_sfo_saem_1$confint_back["k_parent", "est."], 3)) + # Compare fits + expect_known_output(anova(sfo_saem_1, sfo_saem_1_reduced, + sfo_saem_1_mkin, sfo_saem_1_reduced_mkin, test = TRUE), + file = "anova_sfo_saem.txt" + ) + + # 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") - fomc_saem_2 <- saem(mmkin_fomc_1, quiet = TRUE, transformations = "mkin") + 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) @@ -70,45 +78,39 @@ test_that("Parent fits using saemix are correctly implemented", { expect_true(all(ci_fomc_s2[, "lower"] < fomc_pop[2:3])) expect_true(all(ci_fomc_s2[, "upper"] > fomc_pop[2:3])) + # 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_n <- summary(dfop_nlme_1) dfop_pop <- as.numeric(dfop_pop) - expect_true(all(s_dfop_s1$confint_back[, "lower"] < dfop_pop)) - expect_true(all(s_dfop_s1$confint_back[, "upper"] > 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_s2$confint_back[, "lower"] < dfop_pop)) expect_true(all(s_dfop_s2$confint_back[, "upper"] > dfop_pop)) - dfop_mmkin_means_trans_tested <- mean_degparms(mmkin_dfop_1, test_log_parms = TRUE) - dfop_mmkin_means_trans <- apply(parms(mmkin_dfop_1, transformed = TRUE), 1, mean) - - dfop_mmkin_means_tested <- backtransform_odeparms(dfop_mmkin_means_trans_tested, mmkin_dfop_1$mkinmod) - dfop_mmkin_means <- backtransform_odeparms(dfop_mmkin_means_trans, mmkin_dfop_1$mkinmod) - - # We get < 20% deviations for parent_0 and k1 by averaging the transformed parameters - # If we average only parameters passing the t-test, the deviation for k2 is also < 20% - rel_diff_mmkin <- (dfop_mmkin_means - dfop_pop) / dfop_pop - rel_diff_mmkin_tested <- (dfop_mmkin_means_tested - dfop_pop) / dfop_pop - expect_true(all(rel_diff_mmkin[c("parent_0", "k1")] < 0.20)) - expect_true(all(rel_diff_mmkin_tested[c("parent_0", "k1", "k2")] < 0.20)) - - # We get < 20% deviations with transformations made in mkin - rel_diff_1 <- (s_dfop_s1$confint_back[, "est."] - dfop_pop) / 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 expect_true(all(rel_diff_1 < 0.20)) # We get < 20% deviations with transformations made in saemix rel_diff_2 <- (s_dfop_s2$confint_back[, "est."] - dfop_pop) / dfop_pop expect_true(all(rel_diff_2 < 0.2)) - # We use constant error for SFORB because tc is overparameterised (b.1 is ill-defined in saem) - mmkin_sforb_2 <- mmkin("SFORB", ds_dfop, quiet = TRUE, error_model = "const", cores = n_cores) - sforb_saemix_1 <- saem(mmkin_sforb_2, quiet = TRUE, - no_random_effect = c("parent_free_0", "k_parent_free_bound"), - transformations = "saemix") - sforb_saemix_2 <- saem(mmkin_sforb_2, quiet = TRUE, - no_random_effect = c("parent_free_0", "log_k_parent_free_bound"), + # SFORB + mmkin_sforb_1 <- mmkin("SFORB", ds_dfop, quiet = TRUE, cores = n_cores) + sforb_saemix_1 <- saem(mmkin_sforb_1, quiet = TRUE, + no_random_effect = c("parent_free_0"), transformations = "mkin") + sforb_saemix_2 <- saem(mmkin_sforb_1, quiet = TRUE, + no_random_effect = c("parent_free_0"), + transformations = "saemix") expect_equal( log(endpoints(dfop_saemix_1)$distimes[1:2]), log(endpoints(sforb_saemix_1)$distimes[1:2]), tolerance = 0.03) @@ -140,7 +142,7 @@ test_that("We can also use mkin solution methods for saem", { ) 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") + solution_type = "analytical", no_random_effect = "parent_0") 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 diff --git a/tests/testthat/test_tests.R b/tests/testthat/test_tests.R index 39649223..caa0f70a 100644 --- a/tests/testthat/test_tests.R +++ b/tests/testthat/test_tests.R @@ -8,6 +8,7 @@ test_that("The lack-of-fit test works and can be reproduced using nls", { # This code is a slightly modified version of that given in Ritz and Streibig # (2008) Nonlinear Regression using R, p. 64 + f_2_anova <- lm(value ~ as.factor(time), data = subset(DFOP_par_c, name == "parent")) Q <- as.numeric(- 2 * (logLik(f_2_nls) - logLik(f_2_anova))) df.Q <- df.residual(f_2_nls) - df.residual(f_2_anova) p_nls <- 1 - pchisq(Q, df.Q) -- cgit v1.2.1