From df0cff4b829f1abf62f037591a24a8019174dd0a Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Fri, 18 Nov 2022 08:37:40 +0100
Subject: 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.
---
 R/parplot.R | 19 +++++++++++++++----
 R/saem.R    |  3 ++-
 2 files changed, 17 insertions(+), 5 deletions(-)

(limited to 'R')

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
-- 
cgit v1.2.1