diff options
Diffstat (limited to 'tests')
| -rw-r--r-- | tests/testthat/FOCUS_2006_D.csf | 2 | ||||
| -rw-r--r-- | tests/testthat/setup_script.R | 38 | ||||
| -rw-r--r-- | tests/testthat/test_confidence.R | 47 | ||||
| -rw-r--r-- | tests/testthat/test_residuals.R | 8 | ||||
| -rw-r--r-- | tests/testthat/test_tests.R | 27 |
5 files changed, 82 insertions, 40 deletions
diff --git a/tests/testthat/FOCUS_2006_D.csf b/tests/testthat/FOCUS_2006_D.csf index c9da7d61..942d56e1 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-10-28 +Date: 2019-10-31 Optimiser: IRLS [Data] diff --git a/tests/testthat/setup_script.R b/tests/testthat/setup_script.R index cfa978fc..fc972a3d 100644 --- a/tests/testthat/setup_script.R +++ b/tests/testthat/setup_script.R @@ -14,6 +14,27 @@ if (Sys.getenv("TRAVIS") != "") n_cores = 1 # 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)) + +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, @@ -54,8 +75,8 @@ m_synth_DFOP_par <- mkinmod(parent = mkinsub("DFOP", c("M1", "M2")), M2 = mkinsub("SFO"), use_of_ff = "max", quiet = TRUE) -f_SFO_lin_mkin_OLS <- mkinfit(m_synth_SFO_lin, SFO_lin_a, quiet = TRUE) -f_SFO_lin_mkin_ML <- mkinfit(m_synth_SFO_lin, SFO_lin_a, quiet = TRUE, +fit_nw_1 <- mkinfit(m_synth_SFO_lin, SFO_lin_a, quiet = TRUE) +fit_nw_1_ML <- mkinfit(m_synth_SFO_lin, SFO_lin_a, quiet = TRUE, error_model = "const", error_model_algorithm = "direct") # We know direct optimization is OK and direct needs 4 sec versus 5.5 for threestep and 6 for IRLS @@ -69,5 +90,14 @@ fit_tc_1 <- mkinfit(m_synth_SFO_lin, SFO_lin_a, error_model = "tc", quiet = TRUE f_tc_2 <- mkinfit(m_synth_DFOP_par, DFOP_par_c, error_model = "tc", error_model_algorithm = "direct", quiet = TRUE) -#f_tc_2_ntf <- mkinfit(m_synth_DFOP_par, DFOP_par_c, error_model = "tc", -# transform_fractions = FALSE, error_model_algorithm = "direct", quiet = TRUE) +# Experimental data for UBA +dfop_sfo_sfo <- mkinmod( + parent = mkinsub("DFOP", to = "A1"), + A1 = mkinsub("SFO", to = "A2"), + A2 = mkinsub("SFO"), + use_of_ff = "max" +) + +f_soil_1_tc <- mkinfit(dfop_sfo_sfo, + experimental_data_for_UBA_2019[[1]]$data, + error_model = "tc", quiet = TRUE) diff --git a/tests/testthat/test_confidence.R b/tests/testthat/test_confidence.R index 5f76c344..2443fa66 100644 --- a/tests/testthat/test_confidence.R +++ b/tests/testthat/test_confidence.R @@ -1,44 +1,21 @@ -# 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)) - -f_1_mkin_OLS <- mkinfit("SFO", FOCUS_2006_A, quiet = TRUE) -f_1_mkin_OLS_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) - context("Confidence intervals and p-values") test_that("The confint method 'quadratic' is consistent with the summary", { expect_equivalent( - confint(f_SFO_lin_mkin_ML, method = "quadratic"), - summary(f_SFO_lin_mkin_ML)$bpar[, c("Lower", "Upper")]) + confint(fit_nw_1, method = "quadratic"), + summary(fit_nw_1)$bpar[, c("Lower", "Upper")]) expect_equivalent( - confint(f_SFO_lin_mkin_ML, method = "quadratic", backtransform = FALSE), - summary(f_SFO_lin_mkin_ML)$par[, c("Lower", "Upper")]) + confint(fit_nw_1, method = "quadratic", backtransform = FALSE), + summary(fit_nw_1)$par[, c("Lower", "Upper")]) - f_notrans <- mkinfit("SFO", FOCUS_2006_C, quiet = TRUE, transform_rates = FALSE) expect_equivalent( - confint(f_notrans, method = "quadratic", transformed = FALSE), - summary(f_notrans)$par[, c("Lower", "Upper")]) + confint(f_1_mkin_notrans, method = "quadratic", transformed = FALSE), + summary(f_1_mkin_notrans)$par[, c("Lower", "Upper")]) expect_equivalent( - confint(f_notrans, method = "quadratic", transformed = FALSE), - summary(f_notrans)$bpar[, c("Lower", "Upper")]) + confint(f_1_mkin_notrans, method = "quadratic", transformed = FALSE), + summary(f_1_mkin_notrans)$bpar[, c("Lower", "Upper")]) }) @@ -46,20 +23,20 @@ test_that("Quadratic confidence intervals for rate constants are comparable to v # Check fitted parameter values expect_equivalent( - (f_1_mkin_OLS$bparms.optim -coef(f_1_nls_notrans))/f_1_mkin_OLS$bparms.optim, + (f_1_mkin_trans$bparms.optim -coef(f_1_nls_notrans))/f_1_mkin_trans$bparms.optim, rep(0, 2), tolerance = 1e-6) expect_equivalent( - (f_1_mkin_OLS$par[1:2] - coef(f_1_nls_trans))/f_1_mkin_OLS$par[1:2], + (f_1_mkin_trans$par[1:2] - coef(f_1_nls_trans))/f_1_mkin_trans$par[1:2], rep(0, 2), tolerance = 1e-6) # Check the standard error for the transformed parameters se_nls <- summary(f_1_nls_trans)$coefficients[, "Std. Error"] # This is of similar magnitude as the standard error obtained with the mkin - se_mkin <- summary(f_1_mkin_OLS)$par[1:2, "Std. Error"] + se_mkin <- summary(f_1_mkin_trans)$par[1:2, "Std. Error"] se_nls_notrans <- summary(f_1_nls_notrans)$coefficients[, "Std. Error"] # This is also of similar magnitude as the standard error obtained with the mkin - se_mkin_notrans <- summary(f_1_mkin_OLS_notrans)$par[1:2, "Std. Error"] + se_mkin_notrans <- summary(f_1_mkin_notrans)$par[1:2, "Std. Error"] # The difference can partly be explained by the ratio between # the maximum likelihood estimate of the standard error sqrt(rss/n) diff --git a/tests/testthat/test_residuals.R b/tests/testthat/test_residuals.R new file mode 100644 index 00000000..0fe05b4f --- /dev/null +++ b/tests/testthat/test_residuals.R @@ -0,0 +1,8 @@ +context("Residuals extracted from mkinfit models") + +test_that("Residuals are correctly returned", { + f <- fits[["FOMC", "FOCUS_C"]] + expect_equal(residuals(f)[1:3], c(-0.7748906, 2.7090589, -1.9451989)) + + expect_equivalent(residuals(f_tc_2, standardized = TRUE)[1:3], c(0.52579103, 0.40714911, 1.66394233)) +}) diff --git a/tests/testthat/test_tests.R b/tests/testthat/test_tests.R new file mode 100644 index 00000000..523edc4a --- /dev/null +++ b/tests/testthat/test_tests.R @@ -0,0 +1,27 @@ +context("Hypothesis tests") + +test_that("The likelihood ratio test works", { + + expect_error(lrtest(fit_tc_1, f_tc_2), "not been fitted to the same data") + + res <- lrtest(fit_nw_1, fit_tc_1) + expect_equal(res[["2", "Pr(>Chisq)"]], 0.9999998) + +}) + +test_that("We can conveniently fix parameters using 'fixed_parms'", { + f_k2_fixed <- mkinfit("DFOP", FOCUS_2006_C, fixed_parms = c(k2 = 0.05), quiet = TRUE) + expect_equivalent(f_k2_fixed$bparms.ode["k2"], 0.05) +}) + +test_that("Updating fitted models works", { + f_dfop_const <- mkinfit("DFOP", FOCUS_2006_C, quiet = TRUE) + f_dfop_tc <- update(f_dfop_const, error_model = "tc") + + f_soil_1_nw <- update(f_soil_1_tc, error_model = "const") + f_soil_1_nw_A2 <- update(f_soil_1_nw, fixed_parms = c(k_A2 = 0)) + test_nw_tc <- lrtest(f_soil_1_nw, f_soil_1_tc) + expect_equivalent(test_nw_tc[["2", "Pr(>Chisq)"]], 2.113e-6) + test_nw_A2 <- lrtest(f_soil_1_nw, f_soil_1_nw_A2) + expect_equivalent(test_nw_A2[["2", "Pr(>Chisq)"]], 0.9999468) +}) |