From 20b9c584e7c43ecbb708459e531c24a1a4751e17 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Sat, 9 Nov 2019 01:05:51 +0100 Subject: Add a lack-of-fit test - Switch an example dataset in the test setup to a dataset with replicates, adapt tests - Skip the test for lrtest with an update specification as it does not only fail when pkgdown generates static help pages, but also in testthat --- tests/testthat/FOCUS_2006_D.csf | 2 +- tests/testthat/setup_script.R | 10 +++++----- tests/testthat/test_confidence.R | 7 ++++--- tests/testthat/test_tests.R | 22 ++++++++++++++++++++-- 4 files changed, 30 insertions(+), 11 deletions(-) (limited to 'tests') diff --git a/tests/testthat/FOCUS_2006_D.csf b/tests/testthat/FOCUS_2006_D.csf index 528e2b61..09940aa3 100644 --- a/tests/testthat/FOCUS_2006_D.csf +++ b/tests/testthat/FOCUS_2006_D.csf @@ -5,7 +5,7 @@ Description: MeasurementUnits: % AR TimeUnits: days Comments: Created using mkin::CAKE_export -Date: 2019-11-05 +Date: 2019-11-09 Optimiser: IRLS [Data] diff --git a/tests/testthat/setup_script.R b/tests/testthat/setup_script.R index 9becdd2a..e33f4af7 100644 --- a/tests/testthat/setup_script.R +++ b/tests/testthat/setup_script.R @@ -32,9 +32,6 @@ 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) -f_2_mkin <- mkinfit("DFOP", FOCUS_2006_C, quiet = TRUE) -f_2_nls <- nls(value ~ SSbiexp(time, A1, lrc1, A2, lrc2), data = FOCUS_2006_C) - # mmkin object of parent fits for tests models <- c("SFO", "FOMC", "DFOP", "HS") fits <- mmkin(models, @@ -62,11 +59,14 @@ f_sfo_sfo.ff <- mkinfit(SFO_SFO.ff, subset(FOCUS_2006_D, value != 0), quiet = TRUE) -# Two metabolites SFO_lin_a <- synthetic_data_for_UBA_2014[[1]]$data - 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 m_synth_SFO_lin <- mkinmod( parent = mkinsub("SFO", "M1"), M1 = mkinsub("SFO", "M2"), diff --git a/tests/testthat/test_confidence.R b/tests/testthat/test_confidence.R index a2bf1401..e85fdb7a 100644 --- a/tests/testthat/test_confidence.R +++ b/tests/testthat/test_confidence.R @@ -54,11 +54,12 @@ test_that("Quadratic confidence intervals for rate constants are comparable to v # Another case: se_mkin_2 <- summary(f_2_mkin)$par[1:4, "Std. Error"] se_nls_2 <- summary(f_2_nls)$coefficients[, "Std. Error"] - # Here we the ratio of standard errors can be explained by the same + # Here the ratio of standard errors can be explained by the same # principle up to about 3% + nobs_DFOP_par_c_parent <- nrow(subset(DFOP_par_c, name == "parent")) expect_equivalent( se_nls_2[c("lrc1", "lrc2")] / se_mkin_2[c("log_k1", "log_k2")], - rep(sqrt(nrow(FOCUS_2006_C) / (nrow(FOCUS_2006_C) - 4)), 2), + rep(sqrt(nobs_DFOP_par_c_parent / (nobs_DFOP_par_c_parent - 4)), 2), tolerance = 0.03) }) @@ -73,7 +74,7 @@ test_that("Likelihood profile based confidence intervals work", { } f_mle <- stats4::mle(f_nll, start = as.list(parms(f)), nobs = nrow(FOCUS_2006_C)) - ci_mkin_1_p_0.95 <- confint(f, method = "profile", level = 0.95, + ci_mkin_1_p_0.95 <- confint(f, method = "profile", level = 0.95, cores = n_cores, quiet = TRUE) # Magically, we get very similar boundaries as stats4::mle diff --git a/tests/testthat/test_tests.R b/tests/testthat/test_tests.R index 5a522f8e..bdc72f08 100644 --- a/tests/testthat/test_tests.R +++ b/tests/testthat/test_tests.R @@ -1,5 +1,20 @@ context("Hypothesis tests") +test_that("The lack-of-fit test works and can be reproduced using nls", { + + expect_error(loftest(f_1_mkin_trans), "Not defined for fits to data without replicates") + + loftest_mkin <- loftest(f_2_mkin) + + # This code is inspired by Ritz and Streibig (2008) Nonlinear Regression using R, p. 64 + 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) + + expect_equal(loftest_mkin[["2", "Pr(>Chisq)"]], p_nls, tolerance = 1e-5) + +}) + test_that("The likelihood ratio test works", { expect_error(lrtest(f_1_mkin_trans, f_2_mkin), "not been fitted to the same data") @@ -25,7 +40,7 @@ test_that("Updating fitted models works", { parent = mkinsub("DFOP", to = "A1"), A1 = mkinsub("SFO", to = "A2"), A2 = mkinsub("SFO"), - use_of_ff = "max" + use_of_ff = "max", quiet = TRUE ) f_soil_1_tc <- mkinfit(dfop_sfo_sfo, @@ -41,6 +56,9 @@ test_that("Updating fitted models works", { }) test_that("We can do a likelihood ratio test using an update specification", { + skip("This errors out if called by testthat while it works in a normal R session") test_2_mkin_k2 <- lrtest(f_2_mkin, fixed_parms = c(k2 = 0)) - expect_equivalent(test_2_mkin_k2[["2", "Pr(>Chisq)"]], 1.139e-6, tolerance = 1e-8) + expect_equivalent(test_2_mkin_k2[["2", "Pr(>Chisq)"]], 4.851e-8, tolerance = 1e-8) + test_2_mkin_tc <- lrtest(f_2_mkin, error_model = "tc") + expect_equivalent(test_2_mkin_tc[["2", "Pr(>Chisq)"]], 7.302e-5, tolerance = 1e-7) }) -- cgit v1.2.1