aboutsummaryrefslogtreecommitdiff
path: root/tests/testthat
diff options
context:
space:
mode:
Diffstat (limited to 'tests/testthat')
-rw-r--r--tests/testthat/FOCUS_2006_D.csf2
-rw-r--r--tests/testthat/setup_script.R38
-rw-r--r--tests/testthat/test_confidence.R47
-rw-r--r--tests/testthat/test_residuals.R8
-rw-r--r--tests/testthat/test_tests.R27
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)
+})

Contact - Imprint