Fix fitting HS with saemix transformations

4 files changed, 39 insertions, 27 deletions

diff --git a/NEWS.md b/NEWS.md
index 55462fd3..64700704 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,7 @@
 # mkin 1.1.2
 
+- 'R/saem.R': Implement and test saemix transformations for HS.
+
 - 'R/convergence.R': New generic to show convergence information with a method for 'mmin' objects.
 
 - 'R/illparms.R': New generic to show ill-defined parameters with methods for 'mkinfit' and 'mmkin' objects.
diff --git a/R/saem.R b/R/saem.R
index 36997ad7..d2240ea5 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -252,7 +252,7 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk
   }
   degparms_fixed <- object[[1]]$bparms.fixed
 
-  # Transformations are done in the degradation function
+  # Transformations are done in the degradation function by default
   transform.par = rep(0, length(degparms_optim))
 
   odeini_optim_parm_names <- grep('_0$', names(degparms_optim), value = TRUE)
@@ -339,13 +339,24 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk
           }
         }
         if (parent_type == "HS") {
-          model_function <- function(psi, id, xidep) {
-            tb <- exp(psi[id, 4])
-            t <- xidep[, "time"]
-            k1 = exp(psi[id, 2])
-            psi[id, 1] * ifelse(t <= tb,
-              exp(- k1 * t),
-              exp(- k1 * tb) * exp(- exp(psi[id, 3]) * (t - tb)))
+          if (transformations == "mkin") {
+            model_function <- function(psi, id, xidep) {
+              tb <- exp(psi[id, 4])
+              t <- xidep[, "time"]
+              k1 <- exp(psi[id, 2])
+              psi[id, 1] * ifelse(t <= tb,
+                exp(- k1 * t),
+                exp(- k1 * tb) * exp(- exp(psi[id, 3]) * (t - tb)))
+            }
+          } else {
+            model_function <- function(psi, id, xidep) {
+              tb <- exp(psi[id, 4])
+              t <- xidep[, "time"]
+              psi[id, 1] * ifelse(t <= tb,
+                exp(- psi[id, 2] * t),
+                exp(- psi[id, 2] * tb) * exp(- psi[id, 3] * (t - tb)))
+            }
+            transform.par = c(0, 1, 1, 1)
           }
         }
       }
diff --git a/log/test.log b/log/test.log
index 266a566f..776b6f49 100644
--- a/log/test.log
+++ b/log/test.log
@@ -1,13 +1,11 @@
-ℹ Loading mkin
-Loading required package: parallel
 ℹ 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 [31.6s]
+✔ |     1  12 | Dimethenamid data from 2018 [32.0s]
 ────────────────────────────────────────────────────────────────────────────────
 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
@@ -16,7 +14,7 @@ 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]
 ✔ |     1  12 | Nonlinear mixed-effects models [0.2s]
 ────────────────────────────────────────────────────────────────────────────────
@@ -28,13 +26,13 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve
 ✔ |         3 | mkinfit features [0.7s]
 ✔ |         8 | mkinmod model generation and printing [0.2s]
 ✔ |         3 | Model predictions with mkinpredict [0.4s]
-✔ |        16 | Evaluations according to 2015 NAFTA guidance [1.7s]
-✔ |         9 | Nonlinear mixed-effects models with nlme [8.3s]
-✔ |        16 | Plotting [10.6s]
+✔ |        16 | Evaluations according to 2015 NAFTA guidance [1.8s]
+✔ |         9 | Nonlinear mixed-effects models with nlme [9.1s]
+✔ |        16 | Plotting [10.7s]
 ✔ |         4 | Residuals extracted from mkinfit models
-✔ |        25 | saemix parent models [172.9s]
+✔ |        27 | saemix parent models [178.6s]
 ✔ |         2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s]
-✔ |         7 | Fitting the SFORB model [3.7s]
+✔ |         7 | Fitting the SFORB model [3.6s]
 ✔ |         1 | Summaries of old mkinfit objects
 ✔ |         5 | Summary [0.2s]
 ✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.1s]
@@ -42,10 +40,10 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve
 ✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 255.9 s
+Duration: 263.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)
 
-[ FAIL 0 | WARN 0 | SKIP 2 | PASS 225 ]
+[ FAIL 0 | WARN 0 | SKIP 2 | PASS 227 ]
diff --git a/tests/testthat/test_saemix_parent.R b/tests/testthat/test_saemix_parent.R
index 731228d9..7f550999 100644
--- a/tests/testthat/test_saemix_parent.R
+++ b/tests/testthat/test_saemix_parent.R
@@ -10,6 +10,7 @@ test_that("Parent fits using saemix are correctly implemented", {
 
   sfo_saem_2 <- saem(mmkin_sfo_1, quiet = TRUE, transformations = "mkin")
   sfo_saem_3 <- expect_error(saem(mmkin_sfo_2, quiet = TRUE), "at least two parameters")
+  expect_equal(endpoints(sfo_saem_1), endpoints(sfo_saem_2), tol = 0.01)
   s_sfo_s1 <- summary(sfo_saem_1)
   s_sfo_s2 <- summary(sfo_saem_2)
 
@@ -76,20 +77,20 @@ test_that("Parent fits using saemix are correctly implemented", {
 
   mmkin_hs_1 <- mmkin("HS", ds_hs, quiet = TRUE, error_model = "const", cores = n_cores)
   hs_saem_1 <- saem(mmkin_hs_1, quiet = TRUE)
+  hs_saem_2 <- saem(mmkin_hs_1, quiet = TRUE, transformations = "mkin")
+  expect_equal(endpoints(hs_saem_1), endpoints(hs_saem_2), tol = 0.01)
   ci_hs_s1 <- summary(hs_saem_1)$confint_back
 
   hs_pop <- as.numeric(hs_pop)
-  # expect_true(all(ci_hs_s1[, "lower"] < hs_pop)) # k1 is overestimated
+  #expect_true(all(ci_hs_s1[, "lower"] < hs_pop)) # k1 is overestimated
   expect_true(all(ci_hs_s1[, "upper"] > hs_pop))
 
   mmkin_hs_2 <- update(mmkin_hs_1, state.ini = 100, fixed_initials = "parent")
-  hs_saem_2 <- saem(mmkin_hs_2, quiet = TRUE)
-  ci_hs_s2 <- summary(hs_saem_2)$confint_back
+  hs_saem_3 <- saem(mmkin_hs_2, quiet = TRUE)
+  ci_hs_s3 <- summary(hs_saem_3)$confint_back
 
-  #expect_true(all(ci_hs_s2[, "lower"] < hs_pop[2:4])) # k1 again overestimated
-  expect_true(all(ci_hs_s2[, "upper"] > hs_pop[2:4]))
-
-  # HS would likely benefit from implemenation of transformations = "saemix"
+  #expect_true(all(ci_hs_s3[, "lower"] < hs_pop[2:4])) # k1 again overestimated
+  expect_true(all(ci_hs_s3[, "upper"] > hs_pop[2:4]))
 })
 
 test_that("We can also use mkin solution methods for saem", {