Pass error.init to saemix_model, show in parplot

Due to an oversight, error.init was not really passed to saemix_model in
saem.mmkin. The new initial values were reverted to c(1, 1), in order to
avoid changing the test results. Initial values for error model
parameters are now shown in parplot.multistart.

9 files changed, 51 insertions, 29 deletions

diff --git a/R/parplot.R b/R/parplot.R
index 627a4eb8..98579779 100644
--- a/R/parplot.R
+++ b/R/parplot.R
@@ -4,6 +4,10 @@
 #' either by the parameters of the run with the highest likelihood,
 #' or by their medians as proposed in the paper by Duchesne et al. (2021).
 #'
+#' Starting values of degradation model parameters and error model parameters
+#' are shown as green circles. The results obtained in the original run
+#' are shown as red circles.
+#'
 #' @param object The [multistart] object
 #' @param llmin The minimum likelihood of objects to be shown
 #' @param scale By default, scale parameters using the best available fit.
@@ -32,7 +36,7 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best"
 
   orig <- attr(object, "orig")
   orig_parms <- parms(orig)
-  start_parms <- orig$mean_dp_start
+  start_degparms <- orig$mean_dp_start
   all_parms <- parms(object)
 
   if (inherits(object, "multistart.saem.mmkin")) {
@@ -49,18 +53,18 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best"
 
   par(las = 1)
   if (orig$transformations == "mkin") {
-    degparm_names_transformed <- names(start_parms)
+    degparm_names_transformed <- names(start_degparms)
     degparm_index <- which(names(orig_parms) %in% degparm_names_transformed)
     orig_parms[degparm_names_transformed] <- backtransform_odeparms(
       orig_parms[degparm_names_transformed],
       orig$mmkin[[1]]$mkinmod,
       transform_rates = orig$mmkin[[1]]$transform_rates,
       transform_fractions = orig$mmkin[[1]]$transform_fractions)
-    start_parms <- backtransform_odeparms(start_parms,
+    start_degparms <- backtransform_odeparms(start_degparms,
       orig$mmkin[[1]]$mkinmod,
       transform_rates = orig$mmkin[[1]]$transform_rates,
       transform_fractions = orig$mmkin[[1]]$transform_fractions)
-    degparm_names <- names(start_parms)
+    degparm_names <- names(start_degparms)
 
     names(orig_parms) <- c(degparm_names, names(orig_parms[-degparm_index]))
 
@@ -72,6 +76,13 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best"
     colnames(selected_parms)[1:length(degparm_names)] <- degparm_names
   }
 
+  start_errparms <- orig$so@model@error.init
+  names(start_errparms) <- orig$so@model@name.sigma
+
+  start_omegaparms <- orig$so@model@omega.init
+
+  start_parms <- c(start_degparms, start_errparms)
+
   scale <- match.arg(scale)
   parm_scale <- switch(scale,
     best = selected_parms[which.best(object[selected]), ],
diff --git a/R/saem.R b/R/saem.R
index 696ea0ee..c77ff70f 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -149,7 +149,7 @@ saem.mmkin <- function(object,
   covariates = NULL,
   covariate_models = NULL,
   no_random_effect = NULL,
-  error.init = c(3, 0.1),
+  error.init = c(1, 1),
   nbiter.saemix = c(300, 100),
   control = list(displayProgress = FALSE, print = FALSE,
     nbiter.saemix = nbiter.saemix,
@@ -708,6 +708,7 @@ saemix_model <- function(object, solution_type = "auto",
     covariance.model = covariance.model,
     covariate.model = covariate.model,
     omega.init = omega.init,
+    error.init = error.init,
     ...
   )
   attr(res, "mean_dp_start") <- degparms_optim
diff --git a/log/build.log b/log/build.log
index a56a64df..c4f9b8a2 100644
--- a/log/build.log
+++ b/log/build.log
@@ -6,3 +6,4 @@
 * checking for LF line-endings in source and make files and shell scripts
 * checking for empty or unneeded directories
 * building ‘mkin_1.2.0.tar.gz’
+
diff --git a/log/test.log b/log/test.log
index b305bf58..af8e52fd 100644
--- a/log/test.log
+++ b/log/test.log
@@ -1,53 +1,53 @@
 ℹ Testing mkin
 ✔ | F W S  OK | Context
 ✔ |         5 | AIC calculation
-✔ |         5 | Analytical solutions for coupled models [3.5s]
+✔ |         5 | Analytical solutions for coupled models [3.2s]
 ✔ |         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 [34.5s]
+✔ |        12 | Confidence intervals and p-values [1.0s]
+✔ |     1  12 | Dimethenamid data from 2018 [31.4s]
 ────────────────────────────────────────────────────────────────────────────────
 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.3s]
+✔ |        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.2s]
-✔ |        10 | Batch fitting and diagnosing hierarchical kinetic models [25.6s]
-✔ |     1  12 | Nonlinear mixed-effects models [0.4s]
+✔ |        10 | Batch fitting and diagnosing hierarchical kinetic models [23.9s]
+✔ |     1  12 | Nonlinear mixed-effects models [0.3s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip ('test_mixed.R:74'): 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]
-✔ |         9 | Multistart method for saem.mmkin models [40.2s]
-✔ |        16 | Evaluations according to 2015 NAFTA guidance [2.6s]
-✔ |         9 | Nonlinear mixed-effects models with nlme [9.5s]
-✔ |        16 | Plotting [10.5s]
+✔ |         3 | Model predictions with mkinpredict [0.4s]
+✔ |         9 | Multistart method for saem.mmkin models [37.0s]
+✔ |        16 | Evaluations according to 2015 NAFTA guidance [2.5s]
+✔ |         9 | Nonlinear mixed-effects models with nlme [8.8s]
+✔ |        16 | Plotting [10.1s]
 ✔ |         4 | Residuals extracted from mkinfit models
-✔ |     1  36 | saemix parent models [69.8s]
+✔ |     1  36 | saemix parent models [66.0s]
 ────────────────────────────────────────────────────────────────────────────────
 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.6s]
+✔ |         2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s]
 ✔ |        11 | Processing of residue series
-✔ |        10 | Fitting the SFORB model [3.9s]
+✔ |        10 | Fitting the SFORB model [3.7s]
 ✔ |         1 | Summaries of old mkinfit objects
 ✔ |         5 | Summary [0.2s]
-✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3s]
-✔ |         9 | Hypothesis tests [8.5s]
-✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [2.3s]
+✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s]
+✔ |         9 | Hypothesis tests [8.0s]
+✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 225.9 s
+Duration: 211.3 s
 
 ── Skipped tests  ──────────────────────────────────────────────────────────────
 • Fitting this ODE model with saemix takes about 15 minutes on my system (1)
diff --git a/man/parplot.Rd b/man/parplot.Rd
index 37c5841d..ac9e02cf 100644
--- a/man/parplot.Rd
+++ b/man/parplot.Rd
@@ -35,6 +35,11 @@ Produces a boxplot with all parameters from the multiple runs, scaled
 either by the parameters of the run with the highest likelihood,
 or by their medians as proposed in the paper by Duchesne et al. (2021).
 }
+\details{
+Starting values of degradation model parameters and error model parameters
+are shown as green circles. The results obtained in the original run
+are shown as red circles.
+}
 \references{
 Duchesne R, Guillemin A, Gandrillon O, Crauste F. Practical
 identifiability in the frame of nonlinear mixed effects models: the example
diff --git a/man/saem.Rd b/man/saem.Rd
index 11463351..984d341b 100644
--- a/man/saem.Rd
+++ b/man/saem.Rd
@@ -24,7 +24,7 @@ saem(object, ...)
   covariates = NULL,
   covariate_models = NULL,
   no_random_effect = NULL,
-  error.init = c(3, 0.1),
+  error.init = c(1, 1),
   nbiter.saemix = c(300, 100),
   control = list(displayProgress = FALSE, print = FALSE, nbiter.saemix = nbiter.saemix,
     save = FALSE, save.graphs = FALSE),
diff --git a/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg b/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg
index c0332fd5..7017908e 100644
--- a/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg
+++ b/tests/testthat/_snaps/multistart/parplot-for-biphasic-saemix-fit.svg
@@ -173,6 +173,8 @@
 <circle cx='228.40' cy='269.36' r='8.10' style='stroke-width: 0.75; stroke: #61D04F;' />
 <circle cx='270.11' cy='283.47' r='8.10' style='stroke-width: 0.75; stroke: #61D04F;' />
 <circle cx='311.83' cy='285.42' r='8.10' style='stroke-width: 0.75; stroke: #61D04F;' />
+<circle cx='353.54' cy='258.25' r='8.10' style='stroke-width: 0.75; stroke: #61D04F;' />
+<circle cx='395.26' cy='-596.36' r='8.10' style='stroke-width: 0.75; stroke: #61D04F;' />
 <circle cx='103.26' cy='280.60' r='5.40' style='stroke-width: 0.75; stroke: #DF536B;' />
 <circle cx='144.97' cy='281.94' r='5.40' style='stroke-width: 0.75; stroke: #DF536B;' />
 <circle cx='186.69' cy='281.42' r='5.40' style='stroke-width: 0.75; stroke: #DF536B;' />
diff --git a/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg b/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg
index f3373901..18eb7fcc 100644
--- a/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg
+++ b/tests/testthat/_snaps/multistart/parplot-for-sfo-fit.svg
@@ -89,6 +89,8 @@
 <g clip-path='url(#cpNTkuMDR8Njg5Ljc2fDU5LjA0fDUwMi41Ng==)'>
 <circle cx='140.80' cy='280.97' r='8.10' style='stroke-width: 0.75; stroke: #61D04F;' />
 <circle cx='257.60' cy='281.05' r='8.10' style='stroke-width: 0.75; stroke: #61D04F;' />
+<circle cx='374.40' cy='257.89' r='8.10' style='stroke-width: 0.75; stroke: #61D04F;' />
+<circle cx='491.20' cy='-615.24' r='8.10' style='stroke-width: 0.75; stroke: #61D04F;' />
 <circle cx='140.80' cy='280.75' r='5.40' style='stroke-width: 0.75; stroke: #DF536B;' />
 <circle cx='257.60' cy='280.67' r='5.40' style='stroke-width: 0.75; stroke: #DF536B;' />
 <circle cx='374.40' cy='280.50' r='5.40' style='stroke-width: 0.75; stroke: #DF536B;' />
diff --git a/tests/testthat/test_multistart.R b/tests/testthat/test_multistart.R
index 98d3fb6d..502cee98 100644
--- a/tests/testthat/test_multistart.R
+++ b/tests/testthat/test_multistart.R
@@ -3,16 +3,16 @@ context("Multistart method for saem.mmkin models")
 test_that("multistart works for saem.mmkin models", {
   skip_on_cran() # Save CRAN time
   set.seed(123456)
-  saem_sfo_s_multi <- multistart(sfo_saem_1_reduced, n = 8, cores = n_cores,
-    no_random_effect = "parent_0")
+  saem_sfo_s_multi <- multistart(sfo_saem_1_reduced, n = 8, cores = n_cores)
   anova_sfo <- anova(sfo_saem_1,
     sfo_saem_1_reduced,
     best(saem_sfo_s_multi),
     test = TRUE
   )
   # On winbuilder, sfo_saem_1 gives an AIC of 1310.8, while we get 1311.7
-  # locally on Linux and Windows. The other, well-determined fits
-  # both give 1309.7
+  # locally (using saemix 3.2, which likely makes the difference due to the
+  # error parameter patch) on Linux and Windows. The other, well-determined
+  # fits both give 1309.7.
   expect_equal(round(anova_sfo, 1)["sfo_saem_1_reduced", "AIC"], 1309.7)
   expect_equal(round(anova_sfo, 1)["best(saem_sfo_s_multi)", "AIC"], 1309.7)
   expect_true(anova_sfo[3, "Pr(>Chisq)"] > 0.2) # Local: 1, CRAN: 0.4