From ae4ca17b89047052b35acee8e636ff8f31636c13 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 22 Apr 2020 16:09:53 +0200 Subject: Support SFORB with formation fractions --- NEWS.md | 6 ++++- R/mkinfit.R | 9 +++---- R/mkinmod.R | 46 +++++++++++++----------------------- man/mkinmod.Rd | 4 ++-- test.log | 18 +++++++------- tests/testthat/FOCUS_2006_D.csf | 2 +- tests/testthat/test_SFORB.R | 18 ++++++++++++++ tests/testthat/test_mkinfit_errors.R | 3 ++- tests/testthat/test_mkinmod.R | 4 ---- 9 files changed, 59 insertions(+), 51 deletions(-) diff --git a/NEWS.md b/NEWS.md index 790e2e03..ba1daa87 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,8 @@ -# mkin 0.9.49.11 (unreleased) +# mkin 0.9.50.1 (unreleased) + +- Support SFORB with formation fractions + +# mkin 0.9.49.11 (2020-04-20) - Increase a test tolerance to make it pass on all CRAN check machines diff --git a/R/mkinfit.R b/R/mkinfit.R index 1c409569..5c092612 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -378,8 +378,9 @@ mkinfit <- function(mkinmod, observed, if (parmname == "r") parms.ini[parmname] = 0.2 } # Default values for formation fractions in case they are present - for (box in mod_vars) { - f_names <- mkinmod$parms[grep(paste0("^f_", box), mkinmod$parms)] + for (obs_var in obs_vars) { + origin <- mkinmod$map[[obs_var]][[1]] + f_names <- mkinmod$parms[grep(paste0("^f_", origin), mkinmod$parms)] if (length(f_names) > 0) { # We need to differentiate between default and specified fractions # and set the unspecified to 1 - sum(specified)/n_unspecified @@ -388,9 +389,9 @@ mkinfit <- function(mkinmod, observed, sum_f_specified = sum(parms.ini[f_specified_names]) if (sum_f_specified > 1) { stop("Starting values for the formation fractions originating from ", - box, " sum up to more than 1.") + origin, " sum up to more than 1.") } - if (mkinmod$spec[[box]]$sink) n_unspecified = length(f_default_names) + 1 + if (mkinmod$spec[[obs_var]]$sink) n_unspecified = length(f_default_names) + 1 else { n_unspecified = length(f_default_names) } diff --git a/R/mkinmod.R b/R/mkinmod.R index 62f16e73..4587e210 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -182,7 +182,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb "are only supported with formation fractions (use_of_ff = 'max')") } if(spec[[varname]]$sink) { - # If sink is required, add first-order/IORE sink term + # If sink is requested, add first-order/IORE sink term k_compound_sink <- paste("k", box_1, "sink", sep = "_") if(spec[[varname]]$type == "IORE") { k_compound_sink <- paste("k__iore", box_1, "sink", sep = "_") @@ -197,7 +197,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb } else { # otherwise no decline term needed here decline_term = "0" } - } else { + } else { # Maximum use of formation fractions k_compound <- paste("k", box_1, sep = "_") if(spec[[varname]]$type == "IORE") { k_compound <- paste("k__iore", box_1, sep = "_") @@ -236,26 +236,13 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb diffs[[box_1]] <- paste(diffs[[box_1]], "-", decline_term)#}}} if(spec[[varname]]$type == "SFORB") { # {{{ Add SFORB reversible binding terms box_2 = map[[varname]][[2]] - if (use_of_ff == "min") { # Minimum use of formation fractions - k_free_bound <- paste("k", varname, "free", "bound", sep = "_") - k_bound_free <- paste("k", varname, "bound", "free", sep = "_") - parms <- c(parms, k_free_bound, k_bound_free) - reversible_binding_term_1 <- paste("-", k_free_bound, "*", box_1, "+", - k_bound_free, "*", box_2) - reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-", - k_bound_free, "*", box_2) - } else { # Use formation fractions also for the free compartment - stop("The maximum use of formation fractions is not supported for SFORB models") - # The problems were: Calculation of dissipation times did not work in this case - # and the coefficient matrix is not generated correctly by the code present - # in this file in this case - #f_free_bound <- paste("f", varname, "free", "bound", sep = "_") - #k_bound_free <- paste("k", varname, "bound", "free", sep = "_") - #parms <- c(parms, f_free_bound, k_bound_free) - #reversible_binding_term_1 <- paste("+", k_bound_free, "*", box_2) - #reversible_binding_term_2 <- paste("+", f_free_bound, "*", k_compound, "*", box_1, "-", - # k_bound_free, "*", box_2) - } + k_free_bound <- paste("k", varname, "free", "bound", sep = "_") + k_bound_free <- paste("k", varname, "bound", "free", sep = "_") + parms <- c(parms, k_free_bound, k_bound_free) + reversible_binding_term_1 <- paste("-", k_free_bound, "*", box_1, "+", + k_bound_free, "*", box_2) + reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-", + k_bound_free, "*", box_2) diffs[[box_1]] <- paste(diffs[[box_1]], reversible_binding_term_1) diffs[[box_2]] <- paste(diffs[[box_2]], reversible_binding_term_2) } #}}} @@ -286,7 +273,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb k_from_to, "*", origin_box) } else { # Do not introduce a formation fraction if this is the only target - if (spec[[origin_box]]$sink == FALSE && n_targets == 1) { + if (spec[[varname]]$sink == FALSE && n_targets == 1) { diffs[[target_box]] <- paste(diffs[[target_box]], "+", decline_term) } else { @@ -302,7 +289,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb model <- list(diffs = diffs, parms = parms, map = map, spec = spec, use_of_ff = use_of_ff) - # Create coefficient matrix if appropriate#{{{ + # Create coefficient matrix if possible #{{{ if (mat) { boxes <- names(diffs) n <- length(boxes) @@ -321,12 +308,12 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb } else { # off-diagonal elements k.candidate = paste("k", from, to, sep = "_") - if (sub("_free$", "", from) == sub("_bound$", "", to)) { + if (sub("_free$", "", from) == sub("_bound$", "", to)) { k.candidate = paste("k", sub("_free$", "_free_bound", from), sep = "_") - } - if (sub("_bound$", "", from) == sub("_free$", "", to)) { + } + if (sub("_bound$", "", from) == sub("_free$", "", to)) { k.candidate = paste("k", sub("_bound$", "_bound_free", from), sep = "_") - } + } k.effective = intersect(model$parms, k.candidate) m[to, from] = ifelse(length(k.effective) > 0, k.effective, "0") @@ -350,7 +337,8 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb } else { # off-diagonal elements f.candidate = paste("f", from, "to", to, sep = "_") k.candidate = paste("k", from, to, sep = "_") - # SFORB with maximum use of formation fractions not implemented, see above + k.candidate = sub("free.*bound", "free_bound", k.candidate) + k.candidate = sub("bound.*free", "bound_free", k.candidate) m[to, from] = ifelse(f.candidate %in% model$parms, paste(f.candidate, " * k_", from, sep = ""), ifelse(k.candidate %in% model$parms, k.candidate, "0")) diff --git a/man/mkinmod.Rd b/man/mkinmod.Rd index d2b851b6..020917b9 100644 --- a/man/mkinmod.Rd +++ b/man/mkinmod.Rd @@ -60,8 +60,8 @@ A list of class \code{mkinmod} for use with \code{\link{mkinfit}}, The coefficient matrix, if the system of differential equations can be represented by one. } - \item{ll}{ - The likelihood function, taking the parameter vector as the first argument. + \item{cf}{ + If generated, the compiled function as returned by cfunction. } } \description{ diff --git a/test.log b/test.log index 684f6729..a7879c33 100644 --- a/test.log +++ b/test.log @@ -3,16 +3,16 @@ Testing mkin ✔ | OK F W S | Context ✔ | 2 | Export dataset for reading into CAKE ✔ | 13 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [3.5 s] -✔ | 4 | Calculation of FOCUS chi2 error levels [2.2 s] -✔ | 4 | Fitting the SFORB model [1.7 s] +✔ | 4 | Calculation of FOCUS chi2 error levels [2.1 s] +✔ | 6 | Fitting the SFORB model [8.8 s] ✔ | 5 | Calculation of Akaike weights ✔ | 10 | Confidence intervals and p-values [9.4 s] -✔ | 14 | Error model fitting [38.0 s] +✔ | 14 | Error model fitting [37.2 s] ✔ | 6 | Test fitting the decline of metabolites from their maximum [0.8 s] ✔ | 1 | Fitting the logistic model [0.9 s] ✔ | 1 | Test dataset class mkinds used in gmkin ✔ | 12 | Special cases of mkinfit calls [2.3 s] -✔ | 9 | mkinmod model generation and printing [0.2 s] +✔ | 8 | mkinmod model generation and printing [0.2 s] ✔ | 3 | Model predictions with mkinpredict [0.3 s] ✔ | 16 | Evaluations according to 2015 NAFTA guidance [4.1 s] ✔ | 9 | Nonlinear mixed-effects models [11.9 s] @@ -21,15 +21,15 @@ Testing mkin ✔ | 14 | Plotting [4.9 s] ✔ | 4 | AIC calculation ✔ | 4 | Residuals extracted from mkinfit models -✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [5.3 s] +✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [5.4 s] ✔ | 1 | Summaries of old mkinfit objects -✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [7.1 s] -✔ | 9 | Hypothesis tests [36.8 s] +✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [6.9 s] +✔ | 9 | Hypothesis tests [36.7 s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 131.8 s +Duration: 137.8 s -OK: 154 +OK: 155 Failed: 0 Warnings: 0 Skipped: 0 diff --git a/tests/testthat/FOCUS_2006_D.csf b/tests/testthat/FOCUS_2006_D.csf index 81a4ed74..8faa19db 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: 2020-04-20 +Date: 2020-04-22 Optimiser: IRLS [Data] diff --git a/tests/testthat/test_SFORB.R b/tests/testthat/test_SFORB.R index b168a4ee..ad9881a8 100644 --- a/tests/testthat/test_SFORB.R +++ b/tests/testthat/test_SFORB.R @@ -9,4 +9,22 @@ test_that("Fitting the SFORB model is equivalent to fitting DFOP", { expect_match(s_sforb, "Estimated Eigenvalues of SFORB model\\(s\\):") expect_match(s_sforb, "parent_b1 parent_b2") expect_match(s_sforb, "0.45956 *0.01785") + + DFOP_SFO <- mkinmod(parent = mkinsub("DFOP", "M1"), + M1 = mkinsub("SFO"), + use_of_ff = "max", quiet = TRUE) + SFORB_SFO <- mkinmod(parent = mkinsub("SFORB", "M1"), + M1 = mkinsub("SFO"), + use_of_ff = "max", quiet = TRUE) + + SFORB_SFO$coefmat + + f_dfop_sfo <- mkinfit(DFOP_SFO, DFOP_par_c, quiet = TRUE) + f_sforb_sfo <- mkinfit(SFORB_SFO, DFOP_par_c, quiet = TRUE) + f_sforb_sfo_eigen <- mkinfit(SFORB_SFO, DFOP_par_c, solution_type = "eigen", quiet = TRUE) + + expect_equivalent(endpoints(f_sforb_sfo)$distimes, endpoints(f_dfop_sfo)$distimes, + tolerance = 1e-6) + expect_equivalent(endpoints(f_sforb_sfo_eigen)$distimes, endpoints(f_dfop_sfo)$distimes, + tolerance = 1e-6) }) diff --git a/tests/testthat/test_mkinfit_errors.R b/tests/testthat/test_mkinfit_errors.R index 7987d291..940e0b2e 100644 --- a/tests/testthat/test_mkinfit_errors.R +++ b/tests/testthat/test_mkinfit_errors.R @@ -7,6 +7,7 @@ test_that("mkinfit stops to prevent and/or explain user errors", { # 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( @@ -14,7 +15,7 @@ test_that("mkinfit stops to prevent and/or explain user errors", { mkinfit(SFO_SFO.ff.nosink, FOCUS_2006_D, transform_fractions = FALSE, quiet = TRUE), "turn off pathways to sink" ), - "sum of formation fractions") + "sum of formation fractions may exceed one") expect_error(mkinfit(SFO_SFO.ff, FOCUS_2006_D, transform_fractions = TRUE, parms.ini = c(f_parent_to_m1 = 0.5), fixed_parms = "f_parent_to_m1", quiet = TRUE), diff --git a/tests/testthat/test_mkinmod.R b/tests/testthat/test_mkinmod.R index 1bdf094b..3178186a 100644 --- a/tests/testthat/test_mkinmod.R +++ b/tests/testthat/test_mkinmod.R @@ -20,10 +20,6 @@ test_that("mkinmod stops to prevent and/or explain user errors", { expect_error(mkinmod(parent = mkinsub("IORE", "m1"), m1 = mkinsub("SFO"), use_of_ff = "min"), "only supported with formation fractions") - - expect_error(mkinmod(parent = mkinsub("SFORB", "m1"), - m1 = mkinsub("SFO"), use_of_ff = "max"), - "not supported") }) test_that("Printing mkinmod models is reproducible", { -- cgit v1.2.1