aboutsummaryrefslogtreecommitdiff
path: root/tests/testthat
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2019-10-31 01:55:01 +0100
committerJohannes Ranke <jranke@uni-bremen.de>2019-10-31 01:59:05 +0100
commit7091d3738e7e55acb20edb88772b228f6f5b6c98 (patch)
treeb6e31700074605c702662e5238162c57de330453 /tests/testthat
parent5e4ea59a41e00b05ea6664c08c7922e892e8ab77 (diff)
Add likelihood ratio test and other methods, fixes
The likelihood ratio test method is lrtest, in addition, methods for update and residuals were added.
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