From df0cff4b829f1abf62f037591a24a8019174dd0a Mon Sep 17 00:00:00 2001 From: Johannes Ranke 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