Add likelihood ratio test and other methods, fixes

The likelihood ratio test method is lrtest, in addition,
methods for update and residuals were added.

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)
+})