From b36ae3d710858ee3ff2907eb2d780e0dff48a4f3 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Mon, 11 May 2020 13:43:40 +0200 Subject: Analytical solutions for all SFO variants --- R/create_deg_func.R | 62 ++++++++++++++++++++++++-------- R/mkinfit.R | 1 + test.log | 43 +++++++++------------- tests/testthat/DFOP_FOCUS_C_messages.txt | 2 +- tests/testthat/setup_script.R | 22 ++++++------ tests/testthat/test_analytical.R | 45 ++++++++++++++++++++--- tests/testthat/test_from_max_mean.R | 14 ++++---- tests/testthat/test_mkinfit_errors.R | 12 +++---- 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") }) -- cgit v1.2.1