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/test_confidence.R103
2 files changed, 45 insertions, 60 deletions
diff --git a/tests/testthat/FOCUS_2006_D.csf b/tests/testthat/FOCUS_2006_D.csf
index 22a4b125..c9da7d61 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-26
+Date: 2019-10-28
Optimiser: IRLS
[Data]
diff --git a/tests/testthat/test_confidence.R b/tests/testthat/test_confidence.R
index 13652436..5f76c344 100644
--- a/tests/testthat/test_confidence.R
+++ b/tests/testthat/test_confidence.R
@@ -1,3 +1,25 @@
+# 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", {
@@ -20,24 +42,7 @@ test_that("The confint method 'quadratic' is consistent with the summary", {
})
-test_that("Quadratic confidence intervals for rate constants are comparable to confint.nls", {
-
- 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)
+test_that("Quadratic confidence intervals for rate constants are comparable to values in summary.nls", {
# Check fitted parameter values
expect_equivalent(
@@ -70,8 +75,6 @@ test_that("Quadratic confidence intervals for rate constants are comparable to c
expect_equivalent(se_nls_notrans[2] / se_mkin_notrans[2], sqrt(8/5), tolerance = 0.01)
# Another case:
- 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)
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
@@ -82,42 +85,24 @@ test_that("Quadratic confidence intervals for rate constants are comparable to c
tolerance = 0.03)
})
-#test_that("Likelihood profile based confidence intervals work", {
-
-# f <- mkinfit("SFO", FOCUS_2006_C, quiet = TRUE)
-# f <- fits[["SFO", "FOCUS_C"]]
-# f_notrans <- mkinfit("SFO", FOCUS_2006_C, quiet = TRUE, transform_rates = FALSE)
-# pf <- parms(f)
-# f_nll <- function(parent_0, k_parent_sink, sigma) {
-# - f$ll(c(parent_0 = as.numeric(parent_0),
-# k_parent_sink = as.numeric(k_parent_sink),
-# sigma = as.numeric(sigma)))
-# }
-# f_nll(parent_0 = 100, k_parent_sink = 0.3, sigma = 4.7)
-# f_nll(pf[1], pf[2], pf[3])
-# f_mle <- stats4::mle(f_nll, start = as.list(parms(f)), nobs = nrow(FOCUS_2006_C))
-# f_mle <- stats4::mle(f_nll, start = as.list(parms(f)), nobs = 17L)
-# logLik(f_mle)
-# stats4::confint(f_mle, nobs = nrow(FOCUS_2006_C))
-#
-# ci_mkin_1 <- confint(f,
-# method = "quadratic", backtransform = FALSE,
-# transformed = TRUE)
-# ci_maxLik_1 <- maxLik::confint.maxLik(f_maxLik)
-#
-# f_tc_2_maxLik <- maxLik::maxLik(f_tc_2$ll,
-# start = f_tc_2$par)
-# ci_tc_2 <- confint(f_tc_2, method = "quadratic",
-# backtransform = FALSE, transformed = TRUE,
-# distribution = "normal")
-# ci_tc_2_maxLik <- confint(f_tc_2_maxLik)
-# rel_diff_ci_tc_2 <- (ci_tc_2 - ci_tc_2_maxLik)/ci_tc_2_maxLik
-# # The ilr transformed variable appear to misbehave in the numeric differentiation
-# expect_equivalent(
-# rel_diff_ci_tc_2[c(2, 3, 4, 6, 7, 9, 10), ], rep(0, 14),
-# tolerance = 0.05)
-#
-# ci_tc_2[, 1]
-# ci_tc_2_maxLik[, 1]
-
-#})
+test_that("Likelihood profile based confidence intervals work", {
+ f <- fits[["SFO", "FOCUS_C"]]
+
+ # negative log-likelihood for use with mle
+ f_nll <- function(parent_0, k_parent_sink, sigma) {
+ - f$ll(c(parent_0 = as.numeric(parent_0),
+ k_parent_sink = as.numeric(k_parent_sink),
+ sigma = as.numeric(sigma)))
+ }
+ 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, quiet = TRUE)
+
+ # Magically, we get very similar boundaries as stats4::mle
+ # (we need to capture the output to avoid printing this while testing as
+ # stats4::confint uses cat() for its message, instead of message(), so
+ # suppressMessage() has no effect)
+ msg <- capture.output(ci_mle_1_0.95 <- stats4::confint(f_mle, level = 0.95))
+ rel_diff_ci <- (ci_mle_1_0.95 - ci_mkin_1_p_0.95)/ci_mle_1_0.95
+ expect_equivalent(as.numeric(rel_diff_ci), rep(0, 6), tolerance = 1e-4)
+})

Contact - Imprint