Error on unsupported saemix transformations

Also, support FOMC with saemix transformations

4 files changed, 26 insertions, 13 deletions

diff --git a/NEWS.md b/NEWS.md
index 64700704..93f688ff 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,6 @@
 # mkin 1.1.2
 
-- 'R/saem.R': Implement and test saemix transformations for HS.
+- 'R/saem.R': Implement and test saemix transformations for FOMC and HS. Also, error out if saemix transformations are requested but not supported.
 
 - 'R/convergence.R': New generic to show convergence information with a method for 'mmin' objects.
 
diff --git a/R/saem.R b/R/saem.R
index d2240ea5..fe968346 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -253,6 +253,7 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk
   degparms_fixed <- object[[1]]$bparms.fixed
 
   # Transformations are done in the degradation function by default
+  # (transformations = "mkin")
   transform.par = rep(0, length(degparms_optim))
 
   odeini_optim_parm_names <- grep('_0$', names(degparms_optim), value = TRUE)
@@ -276,6 +277,9 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk
     if (length(mkin_model$spec) == 1) {
       parent_type <- mkin_model$spec[[1]]$type
       if (length(odeini_fixed) == 1) {
+        if (transformations == "saemix") {
+          stop("saemix transformations are not supported for parent fits with fixed initial parent value")
+        }
         if (parent_type == "SFO") {
           stop("saemix needs at least two parameters to work on.")
         }
@@ -316,8 +320,15 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk
           }
         }
         if (parent_type == "FOMC") {
-          model_function <- function(psi, id, xidep) {
-            psi[id, 1] / (xidep[, "time"]/exp(psi[id, 3]) + 1)^exp(psi[id, 2])
+          if (transformations == "mkin") {
+            model_function <- function(psi, id, xidep) {
+              psi[id, 1] / (xidep[, "time"]/exp(psi[id, 3]) + 1)^exp(psi[id, 2])
+            }
+          } else {
+            model_function <- function(psi, id, xidep) {
+              psi[id, 1] / (xidep[, "time"]/psi[id, 3] + 1)^psi[id, 2]
+            }
+            transform.par = c(0, 1, 1)
           }
         }
         if (parent_type == "DFOP") {
diff --git a/log/test.log b/log/test.log
index 776b6f49..67ba8841 100644
--- a/log/test.log
+++ b/log/test.log
@@ -5,12 +5,12 @@
 ✔ |         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.0s]
+✔ |     1  12 | Dimethenamid data from 2018 [32.9s]
 ────────────────────────────────────────────────────────────────────────────────
 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 | Error model fitting [5.1s]
+✔ |        14 | Error model fitting [5.0s]
 ✔ |         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]
@@ -27,23 +27,23 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve
 ✔ |         8 | mkinmod model generation and printing [0.2s]
 ✔ |         3 | Model predictions with mkinpredict [0.4s]
 ✔ |        16 | Evaluations according to 2015 NAFTA guidance [1.8s]
-✔ |         9 | Nonlinear mixed-effects models with nlme [9.1s]
-✔ |        16 | Plotting [10.7s]
+✔ |         9 | Nonlinear mixed-effects models with nlme [8.3s]
+✔ |        16 | Plotting [10.6s]
 ✔ |         4 | Residuals extracted from mkinfit models
-✔ |        27 | saemix parent models [178.6s]
+✔ |        28 | saemix parent models [182.2s]
 ✔ |         2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s]
-✔ |         7 | Fitting the SFORB model [3.6s]
+✔ |         7 | Fitting the SFORB model [3.9s]
 ✔ |         1 | Summaries of old mkinfit objects
 ✔ |         5 | Summary [0.2s]
 ✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.1s]
-✔ |         9 | Hypothesis tests [7.8s]
+✔ |         9 | Hypothesis tests [7.7s]
 ✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 263.0 s
+Duration: 266.8 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)
 
-[ FAIL 0 | WARN 0 | SKIP 2 | PASS 227 ]
+[ FAIL 0 | WARN 0 | SKIP 2 | PASS 228 ]
diff --git a/tests/testthat/test_saemix_parent.R b/tests/testthat/test_saemix_parent.R
index 7f550999..5d9a01de 100644
--- a/tests/testthat/test_saemix_parent.R
+++ b/tests/testthat/test_saemix_parent.R
@@ -30,12 +30,14 @@ test_that("Parent fits using saemix are correctly implemented", {
     round(s_sfo_n$confint_back["k_parent", "est."], 3))
 
   mmkin_fomc_1 <- mmkin("FOMC", ds_fomc, quiet = TRUE, error_model = "tc", cores = n_cores)
-  fomc_saem_1 <- saem(mmkin_fomc_1, quiet = TRUE)
+  fomc_saem_1 <- saem(mmkin_fomc_1, quiet = TRUE, transformations = "saemix")
+  fomc_saem_2 <- saem(mmkin_fomc_1, quiet = TRUE, transformations = "mkin")
   ci_fomc_s1 <- summary(fomc_saem_1)$confint_back
 
   fomc_pop <- as.numeric(fomc_pop)
   expect_true(all(ci_fomc_s1[, "lower"] < fomc_pop))
   expect_true(all(ci_fomc_s1[, "upper"] > fomc_pop))
+  expect_equal(endpoints(fomc_saem_1), endpoints(fomc_saem_2), tol = 0.01)
 
   mmkin_fomc_2 <- update(mmkin_fomc_1, state.ini = 100, fixed_initials = "parent")
   fomc_saem_2 <- saem(mmkin_fomc_2, quiet = TRUE, transformations = "mkin")