Show DFOP g equivalent of SFORB in summary, test

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)