aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2020-05-11 13:43:40 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2020-05-11 13:43:40 +0200
commitb36ae3d710858ee3ff2907eb2d780e0dff48a4f3 (patch)
treeb9e075d38233106465481c25b0a777ef043fb1c7
parent576fbc9d86f4db3d1be2fbd4e97b3fcd58f43c2b (diff)
Analytical solutions for all SFO variants
-rw-r--r--R/create_deg_func.R62
-rw-r--r--R/mkinfit.R1
-rw-r--r--test.log43
-rw-r--r--tests/testthat/DFOP_FOCUS_C_messages.txt2
-rw-r--r--tests/testthat/setup_script.R22
-rw-r--r--tests/testthat/test_analytical.R45
-rw-r--r--tests/testthat/test_from_max_mean.R14
-rw-r--r--tests/testthat/test_mkinfit_errors.R12
8 files changed, 128 insertions, 73 deletions
diff --git a/R/create_deg_func.R b/R/create_deg_func.R
index 6c0ae40b..11559799 100644
--- a/R/create_deg_func.R
+++ b/R/create_deg_func.R
@@ -31,18 +31,27 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
predicted_text <- character(0)
+ if (parent_type == "SFO") {
+ if (min_ff) {
+ targets <- c(spec[[1]]$to, if (spec[[1]]$sink) "sink" else NULL)
+ k_parent <- paste(paste0("k_", parent, "_", targets), collapse = " + ")
+ } else {
+ k_parent <- paste0("k_", parent)
+ }
+ }
+
predicted_text[parent] <- paste0(parent_type, ".solution(t, odeini['", parent,
if (parent_type == "SFORB") "_free", "'], ",
switch(parent_type,
- SFO = paste0("k_", parent, if (min_ff) "_sink" else "", ")"),
- FOMC = "alpha, beta)",
- IORE = paste0("k__iore_", parent, if (min_ff) "_sink" else "", ", N_", parent, ")"),
- DFOP = "k1, k2, g)",
- SFORB = paste0("k_", parent, "_free_bound, k_", parent, "_bound_free, k_", parent, "_free", if (min_ff) "_sink" else "", ")"),
- HS = "k1, k2, tb)",
- logistic = "kmax, k0, r)"
- )
- )
+ SFO = k_parent,
+ FOMC = "alpha, beta",
+ IORE = paste0("k__iore_", parent, if (min_ff) "_sink" else "", ", N_", parent),
+ DFOP = "k1, k2, g",
+ SFORB = paste0("k_", parent, "_free_bound, k_", parent, "_bound_free, k_", parent, "_free", if (min_ff) "_sink" else ""),
+ HS = "k1, k2, tb",
+ logistic = "kmax, k0, r"
+ ),
+ ")")
if (length(obs_vars) >= 2) {
supported <- FALSE # except for the following cases
@@ -51,20 +60,43 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
n10 <- paste0("odeini['", parent, "']")
n20 <- paste0("odeini['", n2, "']")
+ # sfo_sfo
+ if (all(spec[[1]]$sink == FALSE, length(obs_vars) == 2,
+ spec[[1]]$type == "SFO", spec[[2]]$type == "SFO")) {
+ supported <- TRUE
+ k1 <- k_parent
+ k2 <- paste0("k_", n2, if(min_ff) "_sink" else "")
+ predicted_text[n2] <- paste0(
+ "(((", k2, "-", k1, ")*", n20, "-", k1, "*", n10, ")*exp(-", k2, "*t)+",
+ k1, "*", n10, "*exp(-", k1, "*t))/(", k2, "-", k1, ")")
+ }
+
+ # sfo_f12_sfo
if (all(use_of_ff == "max", spec[[1]]$sink == TRUE, length(obs_vars) == 2,
spec[[1]]$type == "SFO", spec[[2]]$type == "SFO")) {
supported <- TRUE
- k1 <- paste0("k_", n1)
+ k1 <- k_parent
k2 <- paste0("k_", n2)
f12 <- paste0("f_", n1, "_to_", n2)
- if (parent_type == "SFO") {
- predicted_text[n2] <- paste0(
- "(((", k2, "-", k1, ")*", n20, "-", f12, "*", k1, "*", n10, ")*exp(-", k2, "*t)+",
- f12, "*", k1, "*", n10, "*exp(-", k1, "*t))/(", k2, "-", k1, ")")
- }
+ predicted_text[n2] <- paste0(
+ "(((", k2, "-", k1, ")*", n20, "-", f12, "*", k1, "*", n10, ")*exp(-", k2, "*t)+",
+ f12, "*", k1, "*", n10, "*exp(-", k1, "*t))/(", k2, "-", k1, ")")
+ }
+
+ # sfo_k120_sfo
+ if (all(use_of_ff == "min", spec[[1]]$sink == TRUE, length(obs_vars) == 2,
+ spec[[1]]$type == "SFO", spec[[2]]$type == "SFO")) {
+ supported <- TRUE
+ k12 <- paste0("k_", n1, "_", n2)
+ k10 <- paste0("k_", n1, "_sink")
+ k2 <- paste0("k_", n2, "_sink")
+ predicted_text[n2] <- paste0(
+ "(((", k2, "-", k12, "-", k10, ")*", n20, "-", k12, "*", n10, ")*exp(-", k2, "*t)+",
+ k12, "*", n10, "*exp(-(", k_parent, ")*t))/(", k2, "-(", k_parent, "))")
}
}
+
if (supported) {
deg_func <- function(observed, odeini, odeparms) {}
diff --git a/R/mkinfit.R b/R/mkinfit.R
index 9c3f2baa..683a6b3d 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -297,6 +297,7 @@ mkinfit <- function(mkinmod, observed,
# This is only used for simple decline models
if (length(obs_vars) > 1)
stop("Decline from maximum is only implemented for models with a single observed variable")
+ observed$name <- as.character(observed$name)
means <- aggregate(value ~ time, data = observed, mean, na.rm=TRUE)
t_of_max <- means[which.max(means$value), "time"]
diff --git a/test.log b/test.log
index 74e869da..f2f9ff85 100644
--- a/test.log
+++ b/test.log
@@ -4,44 +4,33 @@ Testing mkin
✔ | 2 | Export dataset for reading into CAKE
✔ | 13 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [1.1 s]
✔ | 4 | Calculation of FOCUS chi2 error levels [0.4 s]
-✔ | 7 | Fitting the SFORB model [3.2 s]
-✔ | 1 | Analytical solutions for coupled models [0.2 s]
+✔ | 7 | Fitting the SFORB model [3.4 s]
+✔ | 4 | Analytical solutions for coupled models [1.7 s]
✔ | 5 | Calculation of Akaike weights
-✔ | 10 | Confidence intervals and p-values [0.9 s]
-✔ | 14 | Error model fitting [3.7 s]
-✔ | 6 | Test fitting the decline of metabolites from their maximum [0.2 s]
+✔ | 10 | Confidence intervals and p-values [1.0 s]
+✔ | 14 | Error model fitting [4.0 s]
+✔ | 4 | Test fitting the decline of metabolites from their maximum [0.2 s]
✔ | 1 | Fitting the logistic model [0.2 s]
✔ | 1 | Test dataset class mkinds used in gmkin
✔ | 12 | Special cases of mkinfit calls [0.6 s]
✔ | 8 | mkinmod model generation and printing [0.2 s]
-✔ | 3 | Model predictions with mkinpredict [0.3 s]
-✔ | 14 2 | Evaluations according to 2015 NAFTA guidance [0.9 s]
-────────────────────────────────────────────────────────────────────────────────
-test_nafta.R:25: skip: Test data from Appendix B are correctly evaluated
-Reason: getRversion() > 4 is TRUE
-
-test_nafta.R:50: skip: Test data from Appendix D are correctly evaluated
-Reason: getRversion() > 4 is TRUE
-────────────────────────────────────────────────────────────────────────────────
-✔ | 9 | Nonlinear mixed-effects models [7.0 s]
-✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2 s]
+✔ | 3 | Model predictions with mkinpredict [0.4 s]
+✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.4 s]
+✔ | 9 | Nonlinear mixed-effects models [7.6 s]
+✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.4 s]
✔ | 3 | Summary
-✔ | 0 1 | Plotting [0.9 s]
-────────────────────────────────────────────────────────────────────────────────
-test_plots_summary_twa.R:98: skip: Plotting mkinfit and mmkin objects is reproducible
-Reason: getRversion() > 4 is TRUE
-────────────────────────────────────────────────────────────────────────────────
+✔ | 14 | Plotting [1.4 s]
✔ | 4 | AIC calculation
✔ | 4 | Residuals extracted from mkinfit models
-✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.3 s]
+✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4 s]
✔ | 1 | Summaries of old mkinfit objects
-✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2 s]
-✔ | 9 | Hypothesis tests [7.2 s]
+✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3 s]
+✔ | 9 | Hypothesis tests [6.5 s]
══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 32.9 s
+Duration: 36.4 s
-OK: 141
+OK: 158
Failed: 0
Warnings: 0
-Skipped: 3
+Skipped: 0
diff --git a/tests/testthat/DFOP_FOCUS_C_messages.txt b/tests/testthat/DFOP_FOCUS_C_messages.txt
index 7284dd19..6aa73e01 100644
--- a/tests/testthat/DFOP_FOCUS_C_messages.txt
+++ b/tests/testthat/DFOP_FOCUS_C_messages.txt
@@ -81,7 +81,7 @@ Sum of squared residuals at call 55: 4.364078
85.016328 -0.776316 -4.027611 1.248897
Sum of squared residuals at call 56: 4.364078
85.016327 -0.776316 -4.027611 1.248897
-Sum of squared residuals at call 57: 4.364078
+Sum of squared residuals at call 57: 4.364077
85.016327 -0.776316 -4.027611 1.248897
85.016327 -0.776316 -4.027611 1.248897
85.008939 -0.777792 -4.026307 1.247720
diff --git a/tests/testthat/setup_script.R b/tests/testthat/setup_script.R
index 86d5089f..58e328cd 100644
--- a/tests/testthat/setup_script.R
+++ b/tests/testthat/setup_script.R
@@ -39,25 +39,27 @@ fits <- mmkin(models,
quiet = TRUE, cores = n_cores)
# One metabolite
-SFO_SFO <- mkinmod(parent = list(type = "SFO", to = "m1"),
- m1 = list(type = "SFO"),
+SFO_SFO <- mkinmod(parent = mkinsub("SFO", to = "m1"),
+ m1 = mkinsub("SFO"),
use_of_ff = "min", quiet = TRUE)
-SFO_SFO.ff <- mkinmod(parent = list(type = "SFO", to = "m1"),
- m1 = list(type = "SFO"),
+SFO_SFO.ff <- mkinmod(parent = mkinsub("SFO", to = "m1"),
+ m1 = mkinsub("SFO"),
use_of_ff = "max", quiet = TRUE)
SFO_SFO.ff.nosink <- mkinmod(
parent = mkinsub("SFO", "m1", sink = FALSE),
m1 = mkinsub("SFO"), quiet = TRUE, use_of_ff = "max")
+FOMC_SFO <- mkinmod(parent = mkinsub("FOMC", to = "m1"),
+ m1 = mkinsub("SFO"), quiet = TRUE)
-f_sfo_sfo_desolve <- mkinfit(SFO_SFO,
- subset(FOCUS_2006_D, value != 0),
+# Avoid warning when fitting a dataset where zero value is removed
+FOCUS_D <- subset(FOCUS_2006_D, value != 0)
+
+f_sfo_sfo_desolve <- mkinfit(SFO_SFO, FOCUS_D,
solution_type = "deSolve", quiet = TRUE)
-f_sfo_sfo_eigen <- mkinfit(SFO_SFO,
- subset(FOCUS_2006_D, value != 0),
+f_sfo_sfo_eigen <- mkinfit(SFO_SFO, FOCUS_D,
solution_type = "eigen", quiet = TRUE)
-f_sfo_sfo.ff <- mkinfit(SFO_SFO.ff,
- subset(FOCUS_2006_D, value != 0),
+f_sfo_sfo.ff <- mkinfit(SFO_SFO.ff, FOCUS_D,
quiet = TRUE)
SFO_lin_a <- synthetic_data_for_UBA_2014[[1]]$data
diff --git a/tests/testthat/test_analytical.R b/tests/testthat/test_analytical.R
index 3d30e042..5972a18a 100644
--- a/tests/testthat/test_analytical.R
+++ b/tests/testthat/test_analytical.R
@@ -1,11 +1,46 @@
context("Analytical solutions for coupled models")
-test_that("The analytical solution of SFO-SFO is correct", {
- f_sfo_sfo.ff.analytical <- mkinfit(SFO_SFO.ff,
- subset(FOCUS_2006_D, value != 0),
- quiet = TRUE)
+test_that("The analytical solutions of SFO-SFO are correct", {
+ # No sink, no formation fractions
+ SFO_SFO_nosink <- mkinmod(
+ parent = mkinsub("SFO", to = "m1", sink = FALSE),
+ m1 = mkinsub("SFO"),
+ use_of_ff = "min", quiet = TRUE)
+ f_sfo_sfo_nosink <- mkinfit(SFO_SFO_nosink, FOCUS_D, quiet = TRUE)
+ f_sfo_sfo_nosink_deSolve <- mkinfit(SFO_SFO_nosink, FOCUS_D,
+ solution_type = "deSolve", quiet = TRUE)
+ expect_equal(
+ parms(f_sfo_sfo_nosink),
+ parms(f_sfo_sfo_nosink_deSolve)
+ )
+
+ # No sink, with formation fractions
+ SFO_SFO.ff_nosink <- mkinmod(
+ parent = mkinsub("SFO", to = "m1", sink = FALSE),
+ m1 = mkinsub("SFO"),
+ use_of_ff = "max", quiet = TRUE)
+ f_sfo_sfo_nosink <- mkinfit(SFO_SFO.ff_nosink, FOCUS_D, quiet = TRUE)
+ f_sfo_sfo_nosink_deSolve <- mkinfit(SFO_SFO.ff_nosink, FOCUS_D,
+ solution_type = "deSolve", quiet = TRUE)
+ expect_equal(
+ parms(f_sfo_sfo_nosink),
+ parms(f_sfo_sfo_nosink_deSolve)
+ )
+
+ # Without formation fraction
+ f_sfo_sfo_analytical <- mkinfit(SFO_SFO, FOCUS_D,
+ solution_type = "analytical", quiet = TRUE)
+ expect_equal(
+ parms(f_sfo_sfo_analytical),
+ parms(f_sfo_sfo_desolve)
+ )
+
+ # With formation fraction
+ f_sfo_sfo.ff_desolve <- mkinfit(SFO_SFO.ff, FOCUS_D,
+ solution_type = "deSolve", quiet = TRUE)
expect_equal(
parms(f_sfo_sfo.ff),
- parms(f_sfo_sfo.ff.analytical)
+ parms(f_sfo_sfo.ff_desolve)
)
+
})
diff --git a/tests/testthat/test_from_max_mean.R b/tests/testthat/test_from_max_mean.R
index 002b39e2..5a66c1ef 100644
--- a/tests/testthat/test_from_max_mean.R
+++ b/tests/testthat/test_from_max_mean.R
@@ -6,15 +6,15 @@ test_that("Fitting from maximum mean value works", {
expect_error(mkinfit(SFO_SFO, FOCUS_2006_D, from_max_mean = TRUE),
"only implemented for models with a single observed variable"),
"Observations with value of zero were removed")
-
+
# We can either explicitly create a model for m1, or subset the data
SFO_m1 <- mkinmod(m1 = mkinsub("SFO"))
- f.1 <- expect_warning(mkinfit(SFO_m1, FOCUS_2006_D, from_max_mean = TRUE, quiet = TRUE))
- expect_equivalent(endpoints(f.1)$distimes["m1", ], c(170.8, 567.5),
- scale = 1, tolerance = 0.1)
+ f.1 <- mkinfit(SFO_m1, FOCUS_D, from_max_mean = TRUE, quiet = TRUE)
+ expect_equivalent(endpoints(f.1)$distimes["m1", ], c(170.8, 567.5),
+ scale = 1, tolerance = 0.1)
- f.2 <- expect_warning(mkinfit("SFO", subset(FOCUS_2006_D, name == "m1"),
- from_max_mean = TRUE, quiet = TRUE))
- expect_equivalent(endpoints(f.2)$distimes["m1", ], c(170.8, 567.5),
+ f.2 <- mkinfit("SFO", subset(FOCUS_D, name == "m1"),
+ from_max_mean = TRUE, quiet = TRUE)
+ expect_equivalent(endpoints(f.2)$distimes["m1", ], c(170.8, 567.5),
scale = 1, tolerance = 0.1)
})
diff --git a/tests/testthat/test_mkinfit_errors.R b/tests/testthat/test_mkinfit_errors.R
index c1e9da1d..85ee574e 100644
--- a/tests/testthat/test_mkinfit_errors.R
+++ b/tests/testthat/test_mkinfit_errors.R
@@ -4,28 +4,24 @@ test_that("mkinfit stops to prevent and/or explain user errors", {
expect_error(mkinfit("foo", FOCUS_2006_A))
expect_error(mkinfit(3, FOCUS_2006_A))
- # We remove zero observations from FOCUS_2006_D beforehand in
- # order to avoid another expect_warning in the code
- FOCUS_2006_D <- subset(FOCUS_2006_D, value != 0)
-
# We get a warning if we use transform_fractions = FALSE with formation fractions
# and an error if any pathway to sink is turned off as well
expect_warning(
expect_error(
- mkinfit(SFO_SFO.ff.nosink, FOCUS_2006_D, transform_fractions = FALSE, quiet = TRUE),
+ mkinfit(SFO_SFO.ff.nosink, FOCUS_D, transform_fractions = FALSE, quiet = TRUE),
"turn off pathways to sink"
),
"sum of formation fractions may exceed one")
- expect_error(mkinfit(SFO_SFO.ff, FOCUS_2006_D, transform_fractions = TRUE,
+ expect_error(mkinfit(SFO_SFO.ff, FOCUS_D, transform_fractions = TRUE,
parms.ini = c(f_parent_to_m1 = 0.5), fixed_parms = "f_parent_to_m1", quiet = TRUE),
"not supported")
- expect_error(mkinfit(SFO_SFO.ff, FOCUS_2006_D,
+ expect_error(mkinfit(SFO_SFO.ff, FOCUS_D,
parms.ini = c(f_parent_to_m1 = 1.1), quiet = TRUE),
"sum up to more than 1")
- expect_error(mkinfit(SFO_SFO, FOCUS_2006_D, solution_type = "analytical"), "not implemented")
+ expect_error(mkinfit(FOMC_SFO, FOCUS_D, solution_type = "analytical"), "not implemented")
expect_error(mkinfit("FOMC", FOCUS_2006_A, solution_type = "eigen"), "coefficient matrix not present")
})

Contact - Imprint