From 56415f3c0638725aaa77db2c02eeaf9d4983452b Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 8 Nov 2022 08:46:57 +0100 Subject: Show DFOP g equivalent of SFORB in summary, test --- R/endpoints.R | 6 ++++-- R/summary.mkinfit.R | 13 +++++++------ log/test.log | 26 +++++++++++++------------- man/summary.mkinfit.Rd | 5 +++-- tests/testthat/test_SFORB.R | 15 +++++++++------ 5 files changed, 36 insertions(+), 29 deletions(-) diff --git a/R/endpoints.R b/R/endpoints.R index 227671b5..4aec8aa8 100644 --- a/R/endpoints.R +++ b/R/endpoints.R @@ -162,6 +162,7 @@ endpoints <- function(fit) { sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 - k_1output * k_21) b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp + g = (k_12 + k_21 - b1)/(b2 - b1) DT50_b1 = log(2)/b1 DT50_b2 = log(2)/b2 @@ -169,8 +170,7 @@ endpoints <- function(fit) { DT90_b2 = log(10)/b2 SFORB_fraction = function(t) { - ((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) + - ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t) + g * exp(-b1 * t) + (1 - g) * exp(-b2 * t) } f_50 <- function(log_t) (SFORB_fraction(exp(log_t)) - 0.5)^2 @@ -195,6 +195,8 @@ endpoints <- function(fit) { # Return the eigenvalues for comparison with DFOP rate constants ep$SFORB[[paste(obs_var, "b1", sep="_")]] = b1 ep$SFORB[[paste(obs_var, "b2", sep="_")]] = b2 + # Return g for comparison with DFOP + ep$SFORB[[paste(obs_var, "g", sep="_")]] = g ep$distimes[obs_var, c("DT50back")] = DT50_back ep$distimes[obs_var, c(paste("DT50", obs_var, "b1", sep = "_"))] = DT50_b1 diff --git a/R/summary.mkinfit.R b/R/summary.mkinfit.R index 4122873f..c25b836e 100644 --- a/R/summary.mkinfit.R +++ b/R/summary.mkinfit.R @@ -6,14 +6,14 @@ #' and optionally the data, consisting of observed, predicted and residual #' values. #' -#' @param object an object of class \code{\link{mkinfit}}. +#' @param object an object of class [mkinfit]. #' @param x an object of class \code{summary.mkinfit}. #' @param data logical, indicating whether the data should be included in the -#' summary. +#' summary. #' @param distimes logical, indicating whether DT50 and DT90 values should be -#' included. +#' included. #' @param alpha error level for confidence interval estimation from t -#' distribution +#' distribution #' @param digits Number of digits to use for printing #' @param \dots optional arguments passed to methods like \code{print}. #' @importFrom stats qt pt cov2cor @@ -37,7 +37,8 @@ #' \item{ff}{The estimated formation fractions derived from the fitted #' model.} #' \item{distimes}{The DT50 and DT90 values for each observed variable.} -#' \item{SFORB}{If applicable, eigenvalues of SFORB components of the model.} +#' \item{SFORB}{If applicable, eigenvalues and fractional eigenvector component +#' g of SFORB systems in the model.} #' The print method is called for its side effect, i.e. printing the summary. #' @author Johannes Ranke #' @references FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence @@ -264,7 +265,7 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . printSFORB <- !is.null(x$SFORB) if(printSFORB){ - cat("\nEstimated Eigenvalues of SFORB model(s):\n") + cat("\nEstimated Eigenvalues and DFOP g parameter of SFORB model(s):\n") print(x$SFORB, digits=digits,...) } diff --git a/log/test.log b/log/test.log index b82e290f..874d6e24 100644 --- a/log/test.log +++ b/log/test.log @@ -1,11 +1,11 @@ ℹ Testing mkin ✔ | F W S OK | Context ✔ | 5 | AIC calculation -✔ | 5 | Analytical solutions for coupled models [3.2s] +✔ | 5 | Analytical solutions for coupled models [3.3s] ✔ | 5 | Calculation of Akaike weights ✔ | 3 | Export dataset for reading into CAKE ✔ | 12 | Confidence intervals and p-values [1.0s] -✔ | 1 12 | Dimethenamid data from 2018 [32.4s] +✔ | 1 12 | Dimethenamid data from 2018 [31.5s] ──────────────────────────────────────────────────────────────────────────────── Skip (test_dmta.R:98:3): Different backends get consistent results for SFO-SFO3+, dimethenamid data Reason: Fitting this ODE model with saemix takes about 15 minutes on my system @@ -14,44 +14,44 @@ Reason: Fitting this ODE model with saemix takes about 15 minutes on my system ✔ | 5 | Time step normalisation ✔ | 4 | Calculation of FOCUS chi2 error levels [0.6s] ✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.8s] -✔ | 4 | Test fitting the decline of metabolites from their maximum [0.3s] +✔ | 4 | Test fitting the decline of metabolites from their maximum [0.4s] ✔ | 1 | Fitting the logistic model [0.2s] -✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [24.5s] +✔ | 10 | Batch fitting and diagnosing hierarchical kinetic models [24.0s] ✔ | 1 12 | Nonlinear mixed-effects models [0.3s] ──────────────────────────────────────────────────────────────────────────────── Skip (test_mixed.R:74:3): saemix results are reproducible for biphasic fits Reason: Fitting with saemix takes around 10 minutes when using deSolve ──────────────────────────────────────────────────────────────────────────────── ✔ | 3 | Test dataset classes mkinds and mkindsg -✔ | 10 | Special cases of mkinfit calls [0.6s] +✔ | 10 | Special cases of mkinfit calls [0.5s] ✔ | 3 | mkinfit features [0.7s] ✔ | 8 | mkinmod model generation and printing [0.2s] ✔ | 3 | Model predictions with mkinpredict [0.3s] ✔ | 7 | Multistart method for saem.mmkin models [36.3s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.5s] -✔ | 9 | Nonlinear mixed-effects models with nlme [9.1s] -✔ | 16 | Plotting [10.3s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [2.4s] +✔ | 9 | Nonlinear mixed-effects models with nlme [8.8s] +✔ | 16 | Plotting [10.1s] ✔ | 4 | Residuals extracted from mkinfit models -✔ | 1 37 | saemix parent models [72.3s] +✔ | 1 37 | saemix parent models [71.5s] ──────────────────────────────────────────────────────────────────────────────── Skip (test_saemix_parent.R:153:3): We can also use mkin solution methods for saem Reason: This still takes almost 2.5 minutes although we do not solve ODEs ──────────────────────────────────────────────────────────────────────────────── ✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s] ✔ | 11 | Processing of residue series -✔ | 7 | Fitting the SFORB model [3.8s] +✔ | 10 | Fitting the SFORB model [3.8s] ✔ | 1 | Summaries of old mkinfit objects ✔ | 5 | Summary [0.2s] ✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s] -✔ | 9 | Hypothesis tests [8.1s] +✔ | 9 | Hypothesis tests [8.0s] ✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 219.0 s +Duration: 216.0 s ── Skipped tests ────────────────────────────────────────────────────────────── • Fitting this ODE model with saemix takes about 15 minutes on my system (1) • Fitting with saemix takes around 10 minutes when using deSolve (1) • This still takes almost 2.5 minutes although we do not solve ODEs (1) -[ FAIL 0 | WARN 0 | SKIP 3 | PASS 265 ] +[ FAIL 0 | WARN 0 | SKIP 3 | PASS 268 ] diff --git a/man/summary.mkinfit.Rd b/man/summary.mkinfit.Rd index b6c1fb87..e315b7ab 100644 --- a/man/summary.mkinfit.Rd +++ b/man/summary.mkinfit.Rd @@ -10,7 +10,7 @@ \method{print}{summary.mkinfit}(x, digits = max(3, getOption("digits") - 3), ...) } \arguments{ -\item{object}{an object of class \code{\link{mkinfit}}.} +\item{object}{an object of class \link{mkinfit}.} \item{data}{logical, indicating whether the data should be included in the summary.} @@ -48,7 +48,8 @@ parameters, for use as starting parameters for related models.} \item{ff}{The estimated formation fractions derived from the fitted model.} \item{distimes}{The DT50 and DT90 values for each observed variable.} -\item{SFORB}{If applicable, eigenvalues of SFORB components of the model.} +\item{SFORB}{If applicable, eigenvalues and fractional eigenvector component +g of SFORB systems in the model.} The print method is called for its side effect, i.e. printing the summary. } \description{ diff --git a/tests/testthat/test_SFORB.R b/tests/testthat/test_SFORB.R index 91c8f2fb..88dcd761 100644 --- a/tests/testthat/test_SFORB.R +++ b/tests/testthat/test_SFORB.R @@ -1,15 +1,20 @@ context("Fitting the SFORB model") -# We do not want the warnings due to non-normality of residuals here -warn_option <- options(warn=-1) - test_that("Fitting the SFORB model is equivalent to fitting DFOP", { f_sforb <- mkinfit("SFORB", FOCUS_2006_C, quiet = TRUE) f_dfop <- mkinfit("DFOP", FOCUS_2006_C, quiet = TRUE) expect_equivalent(endpoints(f_sforb)$distimes, endpoints(f_dfop)$distimes, tolerance = 1e-6) + s_sforb_parms <- summary(f_sforb)$SFORB + expect_equivalent( + exp(f_dfop$par["log_k1"]), s_sforb_parms["parent_b1"]) + expect_equivalent( + exp(f_dfop$par["log_k2"]), s_sforb_parms["parent_b2"]) + expect_equivalent( + plogis(f_dfop$par["g_qlogis"]), s_sforb_parms["parent_g"]) + s_sforb <- capture_output(print(summary(f_sforb))) - expect_match(s_sforb, "Estimated Eigenvalues of SFORB model\\(s\\):") + expect_match(s_sforb, "Estimated Eigenvalues and DFOP g parameter of SFORB model\\(s\\):") expect_match(s_sforb, "parent_b1 parent_b2") expect_match(s_sforb, "0.45956 *0.01785") @@ -35,5 +40,3 @@ test_that("Fitting the SFORB model is equivalent to fitting DFOP", { expect_equivalent(endpoints(f_sforb_sfo_eigen)$distimes, endpoints(f_dfop_sfo)$distimes, tolerance = 1e-6) }) - -options(warn = warn_option$warn) -- cgit v1.2.1