aboutsummaryrefslogtreecommitdiff
path: root/tests/testthat
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2019-10-28 13:23:40 +0100
committerJohannes Ranke <jranke@uni-bremen.de>2019-10-28 13:23:40 +0100
commitbd761d879a95872f82e6d8f893634a61e122a938 (patch)
tree0b66fa5064eb5bd9e0014de18f433c2692fcbc00 /tests/testthat
parent6e5630a0df7e857697ff2ce4730a5f7f45b67377 (diff)
Fix the cutoff for likelihood based intervals
The cutoff now matches what is given by Venzon and Moolgavkar (1988). Also, confidence intervals closely match intervals obtained with stats4::confint in the test case where an stats4::mle object is created from the likelihood function in one test case. Static documentation rebuilt by pkgdown
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