diff options
| -rw-r--r-- | test.log | 41 | ||||
| -rw-r--r-- | tests/testthat/summary_saem_biphasic_s.txt | 3 | ||||
| -rw-r--r-- | tests/testthat/test-dmta.R | 64 | ||||
| -rw-r--r-- | tests/testthat/test_dmta.R | 120 | ||||
| -rw-r--r-- | tests/testthat/test_mixed.R | 93 | ||||
| -rw-r--r-- | tests/testthat/test_saemix_parent.R | 91 |
6 files changed, 236 insertions, 176 deletions
@@ -7,38 +7,41 @@ Loading required package: parallel ✔ | 5 | Calculation of Akaike weights ✔ | 2 | Export dataset for reading into CAKE ✔ | 12 | Confidence intervals and p-values [1.0s] -✔ | 14 | Error model fitting [4.8s] +⠋ | 1 | Dimethenamid data from 2018, parent fits +✔ | 24 | Dimethenamid data from 2018, parent fits [36.7s] +✔ | 14 | Error model fitting [6.8s] ✔ | 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.4s] -✔ | 1 | Fitting the logistic model [0.2s] -✔ | 1 35 | Nonlinear mixed-effects models [26.8s] +✔ | 4 | Calculation of FOCUS chi2 error levels [0.8s] +✔ | 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.5s] +✔ | 1 | Fitting the logistic model [0.3s] +✔ | 1 12 | Nonlinear mixed-effects models [0.2s] ──────────────────────────────────────────────────────────────────────────────── -Skip (test_mixed.R:161:3): saem results are reproducible for biphasic fits +Skip (test_mixed.R:67:3): saemix results are reproducible for biphasic fits Reason: Fitting with saemix takes around 10 minutes when using deSolve ──────────────────────────────────────────────────────────────────────────────── ✔ | 2 | Test dataset classes mkinds and mkindsg -✔ | 10 | Special cases of mkinfit calls [0.4s] -✔ | 1 | mkinfit features [0.3s] +✔ | 10 | Special cases of mkinfit calls [0.6s] +✔ | 1 | mkinfit features [0.5s] ✔ | 8 | mkinmod model generation and printing [0.2s] -✔ | 3 | Model predictions with mkinpredict [0.4s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.5s] -✔ | 9 | Nonlinear mixed-effects models with nlme [8.3s] -✔ | 16 | Plotting [1.3s] +✔ | 3 | Model predictions with mkinpredict [0.3s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.1s] +✔ | 9 | Nonlinear mixed-effects models with nlme [9.4s] +✔ | 16 | Plotting [1.5s] ✔ | 4 | Residuals extracted from mkinfit models -✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5s] -✔ | 7 | Fitting the SFORB model [3.9s] +✔ | 23 | saemix_parent [29.2s] +✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.7s] +✔ | 7 | Fitting the SFORB model [4.4s] ✔ | 1 | Summaries of old mkinfit objects ✔ | 4 | Summary [0.1s] -✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3s] -✔ | 9 | Hypothesis tests [8.6s] +✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.5s] +✔ | 9 | Hypothesis tests [9.5s] ✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 69.5 s +Duration: 115.5 s ── Skipped tests ────────────────────────────────────────────────────────────── • Fitting with saemix takes around 10 minutes when using deSolve (1) -[ FAIL 0 | WARN 0 | SKIP 1 | PASS 206 ] +[ FAIL 0 | WARN 0 | SKIP 1 | PASS 230 ] diff --git a/tests/testthat/summary_saem_biphasic_s.txt b/tests/testthat/summary_saem_biphasic_s.txt index 995e81c8..4569099f 100644 --- a/tests/testthat/summary_saem_biphasic_s.txt +++ b/tests/testthat/summary_saem_biphasic_s.txt @@ -17,7 +17,8 @@ Data: Model predictions using solution type analytical -Fitted in test time 0 s using 300, 100 iterations +Fitted in test time 0 s +Using 300, 100 iterations and 4 chains Variance model: Constant variance diff --git a/tests/testthat/test-dmta.R b/tests/testthat/test-dmta.R deleted file mode 100644 index 12bbcb8e..00000000 --- a/tests/testthat/test-dmta.R +++ /dev/null @@ -1,64 +0,0 @@ -local_edition(3) - -# Data -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"]] <- dmta_ds[["Elliot 2"]] <- NULL - -# mkin -nlm_dfop <- mmkin("DFOP", dmta_ds) -nlm_dfop_tc <- mmkin("DFOP", dmta_ds, error_model = "tc") -parms(nlm_dfop_tc) - -# nlme -nlme_dfop_tc <- nlme(nlm_dfop_tc) -summary(nlme_dfop_tc) -intervals(nlme_dfop_tc) - -# saemix -saem_saemix_dfop_tc <- saem(nlm_dfop_tc) -saem_saemix_dfop_tc$so <- saemix::llgq.saemix(saem_saemix_dfop_tc$so) -summary(saem_saemix_dfop_tc) -intervals(saem_saemix_dfop_tc) -AIC(saem_saemix_dfop_tc$so) -AIC(saem_saemix_dfop_tc$so, "gq") -AIC(saem_saemix_dfop_tc$so, "lin") -saemix::plot(saem_saemix_dfop_tc$so, plot.type = "likelihood") -saemix::plot(saem_saemix_dfop_tc$so, plot.type = "convergence") - -saem_saemix_dfop_tc_1k <- saem(nlm_dfop_tc, nbiter.saemix = c(1000, 100)) -AIC(saem_saemix_dfop_tc_1k$so) -saemix::plot(saem_saemix_dfop_tc_1k$so, plot.type = "convergence") -saemix::plot(saem_saemix_dfop_tc_1k$so, plot.type = "likelihood") -intervals(saem_saemix_dfop_tc_1k) - -saem_saemix_dfop_tc_1.5k <- saem(nlm_dfop_tc, nbiter.saemix = c(1500, 100)) -saem_saemix_dfop_tc_1.5k$so <- saemix::llgq.saemix(saem_saemix_dfop_tc_1.5k$so) -saemix::plot(saem_saemix_dfop_tc_1.5k$so, plot.type = "convergence") -AIC(saem_saemix_dfop_tc_1.5k$so) -AIC(saem_saemix_dfop_tc_1.5k$so, "gq") -intervals(saem_saemix_dfop_tc_1.5k) - -# nlmixr saem -saem_nlmixr_dfop_tc <- nlmixr(nlm_dfop_tc, est = "saem", - control = nlmixr::saemControl(nBurn = 300, nEm = 100, nmc = 9, print = 0)) -intervals(saem_nlmixr_dfop_tc) -summary(saem_nlmixr_dfop_tc) -AIC(saem_nlmixr_dfop_tc$nm) - -saem_nlmixr_dfop_tc_1k <- nlmixr(nlm_dfop_tc, est = "saem", - control = nlmixr::saemControl(nBurn = 1000, nEm = 300, nmc = 9, print = 0)) -intervals(saem_nlmixr_dfop_tc_1k) -summary(saem_nlmixr_dfop_tc_1k) -AIC(saem_nlmixr_dfop_tc_1k$nm) - -focei_nlmixr_dfop_tc <- nlmixr(nlm_dfop_tc, est = "focei") -intervals(focei_nlmixr_dfop_tc) - -AIC(saem_nlmixr_dfop_tc$nm) diff --git a/tests/testthat/test_dmta.R b/tests/testthat/test_dmta.R new file mode 100644 index 00000000..3437966f --- /dev/null +++ b/tests/testthat/test_dmta.R @@ -0,0 +1,120 @@ +context("Dimethenamid data from 2018, parent fits") + +# Data +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"]] <- dmta_ds[["Elliot 2"]] <- NULL + +# mkin +dmta_dfop <- mmkin("DFOP", dmta_ds, quiet = TRUE) +dmta_dfop_tc <- mmkin("DFOP", dmta_ds, error_model = "tc", quiet = TRUE) + +test_that("Different backends get consistent results for DFOP tc, dimethenamid data", { + + # nlme + expect_warning( + nlme_dfop_tc <- nlme(dmta_dfop_tc), + "Iteration 3, .* false convergence") + ints_nlme <- intervals(nlme_dfop_tc) + + # saemix + saem_saemix_dfop_tc <- saem(dmta_dfop_tc) + ints_saemix <- intervals(saem_saemix_dfop_tc) + + # saemix mkin transformations + saem_saemix_dfop_tc_mkin <- saem(dmta_dfop_tc, transformations = "mkin") + ints_saemix_mkin <- intervals(saem_saemix_dfop_tc_mkin) + + # nlmixr saem + saem_nlmixr_dfop_tc <- nlmixr(dmta_dfop_tc, est = "saem", + control = nlmixr::saemControl(nBurn = 300, nEm = 100, nmc = 9, print = 0)) + ints_nlmixr_saem <- intervals(saem_nlmixr_dfop_tc) + + # nlmixr focei + # We get three warnings about nudged etas, the initial optimization and + # gradient problems with initial estimate and covariance + # We need to capture output, otherwise it pops up in testthat output + expect_warning(tmp <- capture_output(focei_nlmixr_dfop_tc <- nlmixr( + dmta_dfop_tc, est = "focei", + control = nlmixr::foceiControl(print = 0), all = TRUE))) + ints_nlmixr_focei <- intervals(focei_nlmixr_dfop_tc) + + # Fixed effects + ## saemix vs. nlme + expect_true(all(ints_saemix$fixed[, "est."] > + backtransform_odeparms(ints_nlme$fixed[, "lower"], dmta_dfop$mkinmod))) + expect_true(all(ints_saemix$fixed[, "est."] < + backtransform_odeparms(ints_nlme$fixed[, "upper"], dmta_dfop$mkinmod))) + + ## saemix mkin vs. nlme + expect_true(all(ints_saemix_mkin$fixed[, "est."] > + backtransform_odeparms(ints_nlme$fixed[, "lower"], dmta_dfop$mkinmod))) + expect_true(all(ints_saemix_mkin$fixed[, "est."] < + backtransform_odeparms(ints_nlme$fixed[, "upper"], dmta_dfop$mkinmod))) + + ## nlmixr saem vs. nlme + expect_true(all(ints_nlmixr_saem$fixed[, "est."] > + backtransform_odeparms(ints_nlme$fixed[, "lower"], dmta_dfop$mkinmod))) + expect_true(all(ints_nlmixr_saem$fixed[, "est."] < + backtransform_odeparms(ints_nlme$fixed[, "upper"], dmta_dfop$mkinmod))) + + ## nlmixr focei vs. nlme + expect_true(all(ints_nlmixr_focei$fixed[, "est."] > + backtransform_odeparms(ints_nlme$fixed[, "lower"], dmta_dfop$mkinmod))) + expect_true(all(ints_nlmixr_focei$fixed[, "est."] < + backtransform_odeparms(ints_nlme$fixed[, "upper"], dmta_dfop$mkinmod))) + + # Random effects + ## for saemix with saemix transformations, the comparison would be complicated... + ## saemix mkin vs. nlme + expect_true(all(ints_saemix$random[, "est."] > + backtransform_odeparms(ints_nlme$reStruct$ds[, "lower"], dmta_dfop$mkinmod))) + expect_true(all(ints_saemix$fixed[, "est."] < + backtransform_odeparms(ints_nlme$fixed[, "upper"], dmta_dfop$mkinmod))) + + ## nlmixr saem vs. nlme + expect_true(all(ints_nlmixr_saem$random[, "est."] > + backtransform_odeparms(ints_nlme$reStruct$ds[, "lower"], dmta_dfop$mkinmod))) + expect_true(all(ints_nlmixr_saem$random[, "est."] < + backtransform_odeparms(ints_nlme$reStruct$ds[, "upper"], dmta_dfop$mkinmod))) + + ## nlmixr focei vs. nlme + expect_true(all(ints_nlmixr_focei$random[, "est."] > + backtransform_odeparms(ints_nlme$reStruct$ds[, "lower"], dmta_dfop$mkinmod))) + expect_true(all(ints_nlmixr_focei$random[, "est."] < + backtransform_odeparms(ints_nlme$reStruct$ds[, "upper"], dmta_dfop$mkinmod))) + + # Variance function + # saemix vs. nlme + expect_true(all(ints_saemix[[3]][, "est."] > + ints_nlme$varStruct[, "lower"])) + expect_true(all(ints_saemix[[3]][, "est."] < + ints_nlme$varStruct[, "upper"])) + + # saemix with mkin transformations vs. nlme + expect_true(all(ints_saemix_mkin[[3]][, "est."] > + ints_nlme$varStruct[, "lower"])) + expect_true(all(ints_saemix_mkin[[3]][, "est."] < + ints_nlme$varStruct[, "upper"])) + + # nlmixr saem vs. nlme + expect_true(all(ints_nlmixr_saem[[3]][, "est."] > + ints_nlme$varStruct[, "lower"])) + expect_true(all(ints_nlmixr_saem[[3]][, "est."] < + ints_nlme$varStruct[, "upper"])) + + # nlmixr focei vs. nlme + # We only test for the proportional part (rsd_high), as the + # constant part (sigma_low) obtained with nlmixr/FOCEI is below the lower + # bound of the confidence interval obtained with nlme + expect_true(ints_nlmixr_focei[[3]]["rsd_high", "est."] > + ints_nlme$varStruct["prop", "lower"]) + expect_true(ints_nlmixr_focei[[3]]["rsd_high", "est."] < + ints_nlme$varStruct["prop", "upper"]) +}) diff --git a/tests/testthat/test_mixed.R b/tests/testthat/test_mixed.R index 40bd3fdf..dbcc66ce 100644 --- a/tests/testthat/test_mixed.R +++ b/tests/testthat/test_mixed.R @@ -1,96 +1,5 @@ context("Nonlinear mixed-effects models") -test_that("Parent fits using saemix are correctly implemented", { - - expect_error(saem(fits), "Only row objects") - # Some fits were done in the setup script - mmkin_sfo_2 <- update(mmkin_sfo_1, fixed_initials = c(parent = 100)) - expect_error(update(mmkin_sfo_1, models = c("SFOOO")), "Please supply models.*") - - 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") - s_sfo_s1 <- summary(sfo_saem_1) - s_sfo_s2 <- summary(sfo_saem_2) - - sfo_nlme_1 <- expect_warning(nlme(mmkin_sfo_1), "not converge") - s_sfo_n <- summary(sfo_nlme_1) - - # Compare with input - expect_equal(round(s_sfo_s2$confint_ranef["SD.log_k_parent", "est."], 1), 0.3) - # k_parent is a bit different from input 0.03 here - expect_equal(round(s_sfo_s1$confint_back["k_parent", "est."], 3), 0.035) - expect_equal(round(s_sfo_s2$confint_back["k_parent", "est."], 3), 0.035) - - # But the result is pretty unanimous between methods - expect_equal(round(s_sfo_s1$confint_back["k_parent", "est."], 3), - round(s_sfo_s2$confint_back["k_parent", "est."], 3)) - expect_equal(round(s_sfo_s1$confint_back["k_parent", "est."], 3), - round(s_sfo_n$confint_back["k_parent", "est."], 3)) - - mmkin_fomc_1 <- mmkin("FOMC", ds_fomc, quiet = TRUE, error_model = "tc", cores = n_cores) - fomc_saem_1 <- saem(mmkin_fomc_1, quiet = TRUE) - ci_fomc_s1 <- summary(fomc_saem_1)$confint_back - - fomc_pop <- as.numeric(fomc_pop) - expect_true(all(ci_fomc_s1[, "lower"] < fomc_pop)) - expect_true(all(ci_fomc_s1[, "upper"] > fomc_pop)) - - mmkin_fomc_2 <- update(mmkin_fomc_1, state.ini = 100, fixed_initials = "parent") - fomc_saem_2 <- saem(mmkin_fomc_2, quiet = TRUE, transformations = "mkin") - ci_fomc_s2 <- summary(fomc_saem_2)$confint_back - - expect_true(all(ci_fomc_s2[, "lower"] < fomc_pop[2:3])) - expect_true(all(ci_fomc_s2[, "upper"] > fomc_pop[2:3])) - - 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)) - 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 < 15% 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.15)) - - # 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)) - - mmkin_hs_1 <- mmkin("HS", ds_hs, quiet = TRUE, error_model = "const", cores = n_cores) - hs_saem_1 <- saem(mmkin_hs_1, quiet = TRUE) - 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_2 <- saem(mmkin_hs_2, quiet = TRUE) - ci_hs_s2 <- summary(hs_saem_2)$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" -}) - test_that("Print methods work", { expect_known_output(print(fits[, 2:3], digits = 2), "print_mmkin_parent.txt") expect_known_output(print(mmkin_biphasic_mixed, digits = 2), "print_mmkin_biphasic_mixed.txt") @@ -122,7 +31,7 @@ test_that("nlme results are reproducible to some degree", { expect_true(all(ci_dfop_sfo_n[, "upper"] > dfop_sfo_pop)) }) -test_that("saem results are reproducible for biphasic fits", { +test_that("saemix results are reproducible for biphasic fits", { test_summary <- summary(saem_biphasic_s) test_summary$saemixversion <- "Dummy 0.0 for testing" diff --git a/tests/testthat/test_saemix_parent.R b/tests/testthat/test_saemix_parent.R new file mode 100644 index 00000000..2f05c175 --- /dev/null +++ b/tests/testthat/test_saemix_parent.R @@ -0,0 +1,91 @@ +test_that("Parent fits using saemix are correctly implemented", { + + expect_error(saem(fits), "Only row objects") + # Some fits were done in the setup script + mmkin_sfo_2 <- update(mmkin_sfo_1, fixed_initials = c(parent = 100)) + expect_error(update(mmkin_sfo_1, models = c("SFOOO")), "Please supply models.*") + + 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") + s_sfo_s1 <- summary(sfo_saem_1) + s_sfo_s2 <- summary(sfo_saem_2) + + sfo_nlme_1 <- expect_warning(nlme(mmkin_sfo_1), "not converge") + s_sfo_n <- summary(sfo_nlme_1) + + # Compare with input + expect_equal(round(s_sfo_s2$confint_ranef["SD.log_k_parent", "est."], 1), 0.3) + # k_parent is a bit different from input 0.03 here + expect_equal(round(s_sfo_s1$confint_back["k_parent", "est."], 3), 0.035) + expect_equal(round(s_sfo_s2$confint_back["k_parent", "est."], 3), 0.035) + + # But the result is pretty unanimous between methods + expect_equal(round(s_sfo_s1$confint_back["k_parent", "est."], 3), + round(s_sfo_s2$confint_back["k_parent", "est."], 3)) + expect_equal(round(s_sfo_s1$confint_back["k_parent", "est."], 3), + round(s_sfo_n$confint_back["k_parent", "est."], 3)) + + mmkin_fomc_1 <- mmkin("FOMC", ds_fomc, quiet = TRUE, error_model = "tc", cores = n_cores) + fomc_saem_1 <- saem(mmkin_fomc_1, quiet = TRUE) + ci_fomc_s1 <- summary(fomc_saem_1)$confint_back + + fomc_pop <- as.numeric(fomc_pop) + expect_true(all(ci_fomc_s1[, "lower"] < fomc_pop)) + expect_true(all(ci_fomc_s1[, "upper"] > fomc_pop)) + + mmkin_fomc_2 <- update(mmkin_fomc_1, state.ini = 100, fixed_initials = "parent") + fomc_saem_2 <- saem(mmkin_fomc_2, quiet = TRUE, transformations = "mkin") + ci_fomc_s2 <- summary(fomc_saem_2)$confint_back + + expect_true(all(ci_fomc_s2[, "lower"] < fomc_pop[2:3])) + expect_true(all(ci_fomc_s2[, "upper"] > fomc_pop[2:3])) + + 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)) + 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 < 15% 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.15)) + + # 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)) + + mmkin_hs_1 <- mmkin("HS", ds_hs, quiet = TRUE, error_model = "const", cores = n_cores) + hs_saem_1 <- saem(mmkin_hs_1, quiet = TRUE) + 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_2 <- saem(mmkin_hs_2, quiet = TRUE) + ci_hs_s2 <- summary(hs_saem_2)$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" +}) + |