Complete starting values in summary for saem.mmkin fits

Also update tests to the changes in mhmkin (see NEWS)

9 files changed, 55 insertions, 45 deletions

diff --git a/NEWS.md b/NEWS.md
index f07c6983..35f499d5 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,8 @@
 # mkin 1.2.2
 
-- 'R/mhmkin.R': Allow an 'illparms.mhmkin' object as value of the argument 'no_random_effects', making it possible to exclude random effects that were ill-defined in simpler variants of the set of degradation models.
+- 'R/mhmkin.R': Allow an 'illparms.mhmkin' object as value of the argument 'no_random_effects', making it possible to exclude random effects that were ill-defined in simpler variants of the set of degradation models. Remove the possibility to exclude random effects based on separate fits, as it did not work well.
+
+- 'R/summary.saem.mmkin.R': List all initial parameter values in the summary, including random effects and error model parameters
 
 # mkin 1.2.1 (2022-11-19)
 
diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R
index 2754e9f0..7337b0f3 100644
--- a/R/summary.saem.mmkin.R
+++ b/R/summary.saem.mmkin.R
@@ -225,13 +225,22 @@ print.summary.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3)
     obs = "Variance unique to each observed variable",
     tc = "Two-component variance function"), "\n")
 
-  cat("\nMean of starting values for individual parameters:\n")
+  cat("\nStarting values for degradation parameters:\n")
   print(x$mean_dp_start, digits = digits)
 
   cat("\nFixed degradation parameter values:\n")
   if(length(x$fixed$value) == 0) cat("None\n")
   else print(x$fixed, digits = digits)
 
+  cat("\nStarting values for random effects (square root of initial entries in omega):\n")
+  print(sqrt(x$so@model@omega.init), digits = digits)
+
+  cat("\nStarting values for error model parameters:\n")
+  errparms <- x$so@model@error.init
+  names(errparms) <- x$so@model@name.sigma
+  errparms <- errparms[x$so@model@indx.res]
+  print(errparms, digits = digits)
+
   cat("\nResults:\n\n")
   cat("Likelihood computed by importance sampling\n")
   print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik,
diff --git a/log/build.log b/log/build.log
index 6be01938..dbe0cd5b 100644
--- a/log/build.log
+++ b/log/build.log
@@ -5,5 +5,5 @@
 * creating vignettes ... OK
 * checking for LF line-endings in source and make files and shell scripts
 * checking for empty or unneeded directories
-* building ‘mkin_1.2.1.tar.gz’
+* building ‘mkin_1.2.2.tar.gz’
 
diff --git a/log/test.log b/log/test.log
index 0b9aa7eb..e17ecc1f 100644
--- a/log/test.log
+++ b/log/test.log
@@ -5,54 +5,54 @@
 ✔ |         5 | Calculation of Akaike weights
 ✔ |         3 | Export dataset for reading into CAKE
 ✔ |        12 | Confidence intervals and p-values [1.1s]
-✔ |     1  12 | Dimethenamid data from 2018 [33.3s]
+✔ |     1  12 | Dimethenamid data from 2018 [32.2s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip ('test_dmta.R:98'): 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.0s]
+✔ |        14 | Error model fitting [4.9s]
 ✔ |         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.4s]
-✔ |         1 | Fitting the logistic model [0.3s]
-✔ |        10 | Batch fitting and diagnosing hierarchical kinetic models [25.5s]
-✔ |     1  11 | Nonlinear mixed-effects models [13.2s]
+✔ |         4 | Test fitting the decline of metabolites from their maximum [0.3s]
+✔ |         1 | Fitting the logistic model [0.2s]
+✔ |         8 | Batch fitting and diagnosing hierarchical kinetic models [14.5s]
+✔ |     1  11 | Nonlinear mixed-effects models [13.1s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip ('test_mixed.R:78'): 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.5s]
+✔ |        10 | Special cases of mkinfit calls [0.4s]
 ✔ |         3 | mkinfit features [0.7s]
 ✔ |         8 | mkinmod model generation and printing [0.2s]
 ✔ |         3 | Model predictions with mkinpredict [0.3s]
-✔ |        12 | Multistart method for saem.mmkin models [51.1s]
-✔ |        16 | Evaluations according to 2015 NAFTA guidance [2.3s]
-✔ |         9 | Nonlinear mixed-effects models with nlme [9.2s]
-✔ |        15 | Plotting [10.3s]
+✔ |        12 | Multistart method for saem.mmkin models [50.1s]
+✔ |        16 | Evaluations according to 2015 NAFTA guidance [2.2s]
+✔ |         9 | Nonlinear mixed-effects models with nlme [8.7s]
+✔ |        15 | Plotting [10.2s]
 ✔ |         4 | Residuals extracted from mkinfit models
-✔ |     1  36 | saemix parent models [73.7s]
+✔ |     1  36 | saemix parent models [103.8s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip ('test_saemix_parent.R:143'): 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.5s]
+✔ |         2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s]
 ✔ |        11 | Processing of residue series
 ✔ |        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.4s]
-✔ |         9 | Hypothesis tests [8.8s]
+✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s]
+✔ |         9 | Hypothesis tests [8.1s]
 ✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 251.2 s
+Duration: 266.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 270 ]
+[ FAIL 0 | WARN 0 | SKIP 3 | PASS 268 ]
 Error while shutting down parallel: unable to terminate some child processes
diff --git a/tests/testthat/illparms_hfits_synth_no_ranef_auto.txt b/tests/testthat/illparms_hfits_synth_no_ranef_auto.txt
deleted file mode 100644
index a64ed222..00000000
--- a/tests/testthat/illparms_hfits_synth_no_ranef_auto.txt
+++ /dev/null
@@ -1,4 +0,0 @@
-           error
-degradation const tc          
-       SFO        sd(parent_0)
-       FOMC       b.1         
diff --git a/tests/testthat/print_hfits_synth_no_ranef_auto.txt b/tests/testthat/print_hfits_synth_no_ranef_auto.txt
deleted file mode 100644
index 9af1cbcd..00000000
--- a/tests/testthat/print_hfits_synth_no_ranef_auto.txt
+++ /dev/null
@@ -1,9 +0,0 @@
-<mhmkin> object
-Status of individual fits:
-
-           error
-degradation const tc
-       SFO  OK    OK
-       FOMC OK    OK
-
-OK: Fit terminated successfully
diff --git a/tests/testthat/summary_hfit_sfo_tc.txt b/tests/testthat/summary_hfit_sfo_tc.txt
index bb5bf6fb..0a61f75f 100644
--- a/tests/testthat/summary_hfit_sfo_tc.txt
+++ b/tests/testthat/summary_hfit_sfo_tc.txt
@@ -17,13 +17,22 @@ Using 300, 100 iterations and 9 chains
 
 Variance model: Two-component variance function 
 
-Mean of starting values for individual parameters:
+Starting values for degradation parameters:
     parent_0 log_k_parent 
          101           -3 
 
 Fixed degradation parameter values:
 None
 
+Starting values for random effects (square root of initial entries in omega):
+             parent_0 log_k_parent
+parent_0            4          0.0
+log_k_parent        0          0.4
+
+Starting values for error model parameters:
+a.1 b.1 
+  1   1 
+
 Results:
 
 Likelihood computed by importance sampling
diff --git a/tests/testthat/summary_saem_dfop_sfo_s.txt b/tests/testthat/summary_saem_dfop_sfo_s.txt
index 7c337843..6468ff17 100644
--- a/tests/testthat/summary_saem_dfop_sfo_s.txt
+++ b/tests/testthat/summary_saem_dfop_sfo_s.txt
@@ -22,7 +22,7 @@ Using 300, 100 iterations and 4 chains
 
 Variance model: Two-component variance function 
 
-Mean of starting values for individual parameters:
+Starting values for degradation parameters:
       parent_0           k_m1 f_parent_to_m1             k1             k2 
          1e+02          7e-03          5e-01          1e-01          2e-02 
              g 
@@ -31,6 +31,19 @@ Mean of starting values for individual parameters:
 Fixed degradation parameter values:
 None
 
+Starting values for random effects (square root of initial entries in omega):
+               parent_0 k_m1 f_parent_to_m1 k1 k2 g
+parent_0            101    0              0  0  0 0
+k_m1                  0    1              0  0  0 0
+f_parent_to_m1        0    0              1  0  0 0
+k1                    0    0              0  1  0 0
+k2                    0    0              0  0  1 0
+g                     0    0              0  0  0 1
+
+Starting values for error model parameters:
+a.1 b.1 
+  1   1 
+
 Results:
 
 Likelihood computed by importance sampling
diff --git a/tests/testthat/test_mhmkin.R b/tests/testthat/test_mhmkin.R
index 93333ac1..e2339f28 100644
--- a/tests/testthat/test_mhmkin.R
+++ b/tests/testthat/test_mhmkin.R
@@ -46,14 +46,4 @@ test_that("Multiple hierarchical kinetic models can be fitted and diagnosed", {
     print(fits_synth_const),
     "print_fits_synth_const.txt")
 
-  hfits_no_ranef_auto <- update(hfits, no_random_effect = "auto", auto_ranef_threshold = 2)
-
-  expect_known_output(
-    print(hfits_no_ranef_auto),
-    "print_hfits_synth_no_ranef_auto.txt")
-
-  expect_known_output(
-    print(illparms(hfits_no_ranef_auto)),
-    "illparms_hfits_synth_no_ranef_auto.txt")
-
 })