From 80ae27dcc4186f03e17164be74158454651474e7 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Wed, 16 Nov 2022 14:50:12 +0100
Subject: Read in all data per default

---
 R/read_spreadsheet.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

(limited to 'R')

diff --git a/R/read_spreadsheet.R b/R/read_spreadsheet.R
index a20af6db..7ad09c3e 100644
--- a/R/read_spreadsheet.R
+++ b/R/read_spreadsheet.R
@@ -37,7 +37,7 @@
 #' and moisture normalisation factors in the sheet 'Datasets'?
 #' @export
 read_spreadsheet <- function(path, valid_datasets = "all",
-  parent_only = TRUE, normalize = TRUE)
+  parent_only = FALSE, normalize = TRUE)
 {
   if (!requireNamespace("readxl", quietly = TRUE))
     stop("Please install the readxl package to use this function")
-- 
cgit v1.2.1


From 9db8338d3a9240ed4685fcdd7aab9692031d5a04 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Fri, 18 Nov 2022 05:56:25 +0100
Subject: Improve logLik.mkinfit to attach nobs attribute

The lack of that attribute lead to a failure to calculate the BIC in
test_AIC.R on R-devel from yesterday.
---
 R/logLik.mkinfit.R | 1 +
 1 file changed, 1 insertion(+)

(limited to 'R')

diff --git a/R/logLik.mkinfit.R b/R/logLik.mkinfit.R
index 7cc10234..abf8517c 100644
--- a/R/logLik.mkinfit.R
+++ b/R/logLik.mkinfit.R
@@ -37,6 +37,7 @@ logLik.mkinfit <- function(object, ...) {
   val <- object$logLik
   # Number of estimated parameters
   attr(val, "df") <- length(object$bparms.optim) + length(object$errparms)
+  attr(val, "nobs") <- nobs(object)
   class(val) <- "logLik"
   return(val)
 }
-- 
cgit v1.2.1


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


From 5364f037a72863ef5ba81e14ba4417f68fd389f9 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Fri, 18 Nov 2022 19:14:47 +0100
Subject: Make mixed model test data permanent to ensure reproducibility

To ensure that tests on different platforms work on the same data, the
mixed modelling test data previosly generated in
tests/testthat/setup_script.R were generated once using the script in
inst/dataset/generation/ds_mixed.R, and are now distributed with the
package.
---
 R/ds_mixed.R | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)
 create mode 100644 R/ds_mixed.R

(limited to 'R')

diff --git a/R/ds_mixed.R b/R/ds_mixed.R
new file mode 100644
index 00000000..c5055712
--- /dev/null
+++ b/R/ds_mixed.R
@@ -0,0 +1,17 @@
+#'  Synthetic data for hierarchical kinetic degradation models
+#'
+#' The R code used to create this data object is installed with this package in
+#' the 'dataset_generation' directory.
+#'
+#' @name ds_mixed
+#' @aliases ds_sfo ds_fomc ds_dfop ds_hs ds_dfop_sfo
+#' @examples
+#' \dontrun{
+#'   sfo_mmkin <- mmkin("SFO", ds_sfo, quiet = TRUE, error_model = "tc", cores = 15)
+#'   sfo_saem <- saem(sfo_mmkin, no_random_effect = "parent_0")
+#'   plot(sfo_saem)
+#' }
+#'
+#' # This is the code used to generate the datasets
+#' cat(readLines(system.file("dataset_generation/ds_mixed.R", package = "mkin")), sep = "\n")
+NULL
-- 
cgit v1.2.1


From 61b9a4046582da5cf88bd3c04d0d6ca77adc3405 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Fri, 25 Nov 2022 19:10:28 +0100
Subject: mhmkin: Easy specification of ill-defined parms

The argument 'no_random_effect' now accepts an illparms.mhmkin object
---
 R/mhmkin.R | 38 ++++++++++++++++++++------------------
 1 file changed, 20 insertions(+), 18 deletions(-)

(limited to 'R')

diff --git a/R/mhmkin.R b/R/mhmkin.R
index 1f29dc40..6a61e8bb 100644
--- a/R/mhmkin.R
+++ b/R/mhmkin.R
@@ -12,13 +12,13 @@
 #' degradation models to the same data
 #' @param backend The backend to be used for fitting. Currently, only saemix is
 #' supported
-#' @param no_random_effect Default is NULL and will be passed to [saem]. If
-#' you specify "auto", random effects are only included if the number
-#' of datasets in which the parameter passed the t-test is at least 'auto_ranef_threshold'.
-#' Beware that while this may make for convenient model reduction or even
-#' numerical stability of the algorithm, it will likely lead to
-#' underparameterised models.
-#' @param auto_ranef_threshold See 'no_random_effect.
+#' @param no_random_effect Default is NULL and will be passed to [saem]. If a
+#' character vector is supplied, it will be passed to all calls to [saem],
+#' regardless if the corresponding parameter is in the model. Alternatively,
+#' an object of class [illparms.mhmkin] can be specified. This has to have
+#' the same dimensions as the return object of the current call. In this way,
+#' ill-defined parameters found in corresponding simpler versions of the
+#' degradation model can be specified.
 #' @param algorithm The algorithm to be used for fitting (currently not used)
 #' @param \dots Further arguments that will be passed to the nonlinear mixed-effects
 #' model fitting function.
@@ -51,7 +51,7 @@ mhmkin.mmkin <- function(objects, ...) {
 #' @export
 #' @rdname mhmkin
 mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
-  no_random_effect = NULL, auto_ranef_threshold = 3,
+  no_random_effect = NULL,
   ...,
   cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(), cluster = NULL)
 {
@@ -102,22 +102,24 @@ mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
     deg_index <- w[1]
     error_index <- w[2]
     mmkin_row <- objects[[error_index]][deg_index, ]
-    if (identical(no_random_effect, "auto")) {
-      ip <- illparms(mmkin_row)
-      ipt <- table(unlist(lapply(as.list(ip), strsplit, ", ")))
-      no_ranef <- names(ipt)[which((length(ds) - ipt) <= auto_ranef_threshold)]
-      transparms <- rownames(mmkin_row[[1]]$start_transformed)
-      bparms <- rownames(mmkin_row[[1]]$start)
-      names(transparms) <- bparms
-      no_ranef_trans <- transparms[no_ranef]
+    if (is(no_random_effect, "illparms.mhmkin")) {
+      if (identical(dim(no_random_effect), dim(fit_indices))) {
+        no_ranef_split <- strsplit(no_random_effect[[fit_index]], ", ")
+        no_ranef <- sapply(no_ranef_split, function(x) {
+          gsub("sd\\((.*)\\)", "\\1", x)
+        })
+      } else {
+        stop("Dimensions of illparms.mhmkin object given as 'no_random_effect' are not suitable")
+      }
     } else {
-      no_ranef_trans <- no_random_effect
+      no_ranef <- no_random_effect
     }
     res <- try(do.call(backend_function,
-        args = c(list(mmkin_row), dot_args, list(no_random_effect = no_ranef_trans))))
+        args = c(list(mmkin_row), dot_args, list(no_random_effect = no_ranef))))
     return(res)
   }
 
+
   fit_time <- system.time({
     if (is.null(cluster)) {
       results <- parallel::mclapply(as.list(1:n.fits), fit_function,
-- 
cgit v1.2.1


From aaa4cab7e0c7212f91147a9789af54b97fe342ca Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Tue, 29 Nov 2022 20:23:17 +0100
Subject: Complete starting values in summary for saem.mmkin fits

Also update tests to the changes in mhmkin (see NEWS)
---
 R/summary.saem.mmkin.R | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

(limited to 'R')

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


From 74e44dfed5af6e6fd421abe82d3e3f190771f85a Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Thu, 1 Dec 2022 11:20:00 +0100
Subject: Possibility to manually specify no_random_effects in mhmkin

---
 R/mhmkin.R | 91 ++++++++++++++++++++++++++++++++++++++++++++++++--------------
 1 file changed, 71 insertions(+), 20 deletions(-)

(limited to 'R')

diff --git a/R/mhmkin.R b/R/mhmkin.R
index 6a61e8bb..a9798ff4 100644
--- a/R/mhmkin.R
+++ b/R/mhmkin.R
@@ -14,11 +14,12 @@
 #' supported
 #' @param no_random_effect Default is NULL and will be passed to [saem]. If a
 #' character vector is supplied, it will be passed to all calls to [saem],
-#' regardless if the corresponding parameter is in the model. Alternatively,
-#' an object of class [illparms.mhmkin] can be specified. This has to have
-#' the same dimensions as the return object of the current call. In this way,
-#' ill-defined parameters found in corresponding simpler versions of the
-#' degradation model can be specified.
+#' which will exclude random effects for all matching parameters. Alternatively,
+#' a list of character vectors or an object of class [illparms.mhmkin] can be
+#' specified. They have to have the same dimensions that the return object of
+#' the current call will have, i.e. the number of rows must match the number
+#' of degradation models in the mmkin object(s), and the number of columns must
+#' match the number of error models used in the mmkin object(s).
 #' @param algorithm The algorithm to be used for fitting (currently not used)
 #' @param \dots Further arguments that will be passed to the nonlinear mixed-effects
 #' model fitting function.
@@ -50,6 +51,42 @@ mhmkin.mmkin <- function(objects, ...) {
 
 #' @export
 #' @rdname mhmkin
+#' @examples
+#' \dontrun{
+#' # We start with separate evaluations of all the first six datasets with two
+#' # degradation models and two error models
+#' f_sep_const <- mmkin(c("SFO", "FOMC"), ds_fomc[1:6], cores = 2, quiet = TRUE)
+#' f_sep_tc <- update(f_sep_const, error_model = "tc")
+#' # The mhmkin function sets up hierarchical degradation models aka
+#' # nonlinear mixed-effects models for all four combinations, specifying
+#' # uncorrelated random effects for all degradation parameters
+#' f_saem_1 <- mhmkin(list(f_sep_const, f_sep_tc), cores = 2)
+#' status(f_saem_1)
+#' # The 'illparms' function shows that in all hierarchical fits, at least
+#' # one random effect is ill-defined (the confidence interval for the
+#' # random effect expressed as standard deviation includes zero)
+#' illparms(f_saem_1)
+#' # Therefore we repeat the fits, excluding the ill-defined random effects
+#' f_saem_2 <- update(f_saem_1, no_random_effect = illparms(f_saem_1))
+#' status(f_saem_2)
+#' illparms(f_saem_2)
+#' # Model comparisons show that FOMC with two-component error is preferable,
+#' # and confirms our reduction of the default parameter model
+#' anova(f_saem_1)
+#' anova(f_saem_2)
+#' # The convergence plot for the selected model looks fine
+#' saemix::plot(f_saem_2[["FOMC", "tc"]]$so, plot.type = "convergence")
+#' # The plot of predictions versus data shows that we have a pretty data-rich
+#' # situation with homogeneous distribution of residuals, because we used the
+#' # same degradation model, error model and parameter distribution model that
+#' # was used in the data generation.
+#' plot(f_saem_2[["FOMC", "tc"]])
+#' # We can specify the same parameter model reductions manually
+#' no_ranef <- list("parent_0", "log_beta", "parent_0", c("parent_0", "log_beta"))
+#' dim(no_ranef) <- c(2, 2)
+#' f_saem_2m <- update(f_saem_1, no_random_effect = no_ranef)
+#' anova(f_saem_2m)
+#' }
 mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
   no_random_effect = NULL,
   ...,
@@ -97,25 +134,38 @@ mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
   dimnames(fit_indices) <- list(degradation = names(deg_models),
                                 error = error_models)
 
-  fit_function <- function(fit_index) {
-    w <- which(fit_indices == fit_index, arr.ind = TRUE)
-    deg_index <- w[1]
-    error_index <- w[2]
-    mmkin_row <- objects[[error_index]][deg_index, ]
+  if (is.null(no_random_effect) || length(dim(no_random_effect)) == 1) {
+    no_ranef <- rep(list(no_random_effect), n.fits)
+    dim(no_ranef) <- dim(fit_indices)
+  } else {
+    if (!identical(dim(no_random_effect), dim(fit_indices))) {
+      stop("Dimensions of argument 'no_random_effect' are not suitable")
+    }
     if (is(no_random_effect, "illparms.mhmkin")) {
-      if (identical(dim(no_random_effect), dim(fit_indices))) {
-        no_ranef_split <- strsplit(no_random_effect[[fit_index]], ", ")
-        no_ranef <- sapply(no_ranef_split, function(x) {
-          gsub("sd\\((.*)\\)", "\\1", x)
+      no_ranef_dim <- dim(no_random_effect)
+      no_ranef <- lapply(no_random_effect, function(x) {
+        no_ranef_split <- strsplit(x, ", ")
+        ret <- sapply(no_ranef_split, function(y) {
+          gsub("sd\\((.*)\\)", "\\1", y)
         })
-      } else {
-        stop("Dimensions of illparms.mhmkin object given as 'no_random_effect' are not suitable")
-      }
+        return(ret)
+      })
+      dim(no_ranef) <- no_ranef_dim
     } else {
       no_ranef <- no_random_effect
     }
+  }
+
+  fit_function <- function(fit_index) {
+    w <- which(fit_indices == fit_index, arr.ind = TRUE)
+    deg_index <- w[1]
+    error_index <- w[2]
+    mmkin_row <- objects[[error_index]][deg_index, ]
     res <- try(do.call(backend_function,
-        args = c(list(mmkin_row), dot_args, list(no_random_effect = no_ranef))))
+        args = c(
+          list(mmkin_row),
+          dot_args,
+          list(no_random_effect = no_ranef[[deg_index, error_index]]))))
     return(res)
   }
 
@@ -145,15 +195,16 @@ mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
 #' @param j Column index selecting the fits to specific datasets
 #' @param drop If FALSE, the method always returns an mhmkin object, otherwise
 #'   either a list of fit objects or a single fit object.
-#' @return An object of class \code{\link{mhmkin}}.
+#' @return An object inheriting from \code{\link{mhmkin}}.
 #' @rdname mhmkin
 #' @export
 `[.mhmkin` <- function(x, i, j, ..., drop = FALSE) {
+  original_class <- class(x)
   class(x) <- NULL
   x_sub <- x[i, j, drop = drop]
 
   if (!drop) {
-    class(x_sub) <- "mhmkin"
+    class(x_sub) <- original_class
   }
   return(x_sub)
 }
-- 
cgit v1.2.1


From cf73354f985be17352a4d538c6f9ea69f6be9aa9 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Fri, 2 Dec 2022 13:28:26 +0100
Subject: Avoid redundant warnings in summaries

---
 R/summary.saem.mmkin.R | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

(limited to 'R')

diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R
index 7337b0f3..46ab548b 100644
--- a/R/summary.saem.mmkin.R
+++ b/R/summary.saem.mmkin.R
@@ -136,10 +136,11 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes =
   }
 
   #  Correlation of fixed effects (inspired by summary.nlme)
-  varFix <- try(vcov(object$so)[1:n_fixed, 1:n_fixed])
-  if (inherits(varFix, "try-error")) {
+  cov_so <- try(solve(object$so@results@fim), silent = TRUE)
+  if (inherits(cov_so, "try-error")) {
     object$corFixed <- NA
   } else {
+    varFix <- cov_so[1:n_fixed, 1:n_fixed]
     stdFix <- sqrt(diag(varFix))
     object$corFixed <- array(
       t(varFix/stdFix)/stdFix,
-- 
cgit v1.2.1


From 478c6d5eec4c84b22b43adcbdf36888b302ead00 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Tue, 6 Dec 2022 10:33:24 +0100
Subject: Some parplot improvements

llquant argument, improved legend text, tests
---
 R/parplot.R | 16 ++++++++++++----
 1 file changed, 12 insertions(+), 4 deletions(-)

(limited to 'R')

diff --git a/R/parplot.R b/R/parplot.R
index 98579779..63306ac2 100644
--- a/R/parplot.R
+++ b/R/parplot.R
@@ -10,7 +10,10 @@
 #'
 #' @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.
+#' @param llquant Fractional value for selecting only the fits with higher
+#' likelihoods. Overrides 'llmin'.
+#' @param scale By default, scale parameters using the best
+#' available fit.
 #' If 'median', parameters are scaled using the median parameters from all fits.
 #' @param main Title of the plot
 #' @param lpos Positioning of the legend.
@@ -28,7 +31,8 @@ parplot <- function(object, ...) {
 
 #' @rdname parplot
 #' @export
-parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best", "median"),
+parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA,
+  scale = c("best", "median"),
   lpos = "bottomleft", main = "", ...)
 {
   oldpar <- par(no.readonly = TRUE)
@@ -48,6 +52,10 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best"
     stop("parplot is only implemented for multistart.saem.mmkin objects")
   }
   ll <- sapply(object, llfunc)
+  if (!is.na(llquant[1])) {
+    if (llmin != -Inf) warning("Overriding 'llmin' because 'llquant' was specified")
+    llmin <- quantile(ll, 1 - llquant)
+  }
   selected <- which(ll > llmin)
   selected_parms <- all_parms[selected, ]
 
@@ -110,7 +118,7 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best"
   legend(lpos, inset = c(0.05, 0.05), bty = "n",
     pch = 1, col = 3:1, lty = c(NA, NA, 1),
     legend = c(
-      "Starting parameters",
-      "Original run",
+      "Original start",
+      "Original results",
       "Multistart runs"))
 }
-- 
cgit v1.2.1


From 97f71fc3d086bd447ab3e4d19abf32bb3114085b Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Wed, 7 Dec 2022 14:41:52 +0100
Subject: Check slopes in saemix covariate models

---
 R/illparms.R | 14 +++++++++++++-
 1 file changed, 13 insertions(+), 1 deletion(-)

(limited to 'R')

diff --git a/R/illparms.R b/R/illparms.R
index 01e75cf1..eef4bd33 100644
--- a/R/illparms.R
+++ b/R/illparms.R
@@ -20,6 +20,9 @@
 #' @param random For hierarchical fits, should random effects be tested?
 #' @param errmod For hierarchical fits, should error model parameters be
 #' tested?
+#' @param slopes For hierarchical [saem] fits using saemix as backend,
+#' should slope parameters in the covariate model(starting with 'beta_') be
+#' tested?
 #' @return For [mkinfit] or [saem] objects, a character vector of parameter
 #' names. For [mmkin] or [mhmkin] objects, a matrix like object of class
 #' 'illparms.mmkin' or 'illparms.mhmkin'.
@@ -92,7 +95,7 @@ print.illparms.mmkin <- function(x, ...) {
 
 #' @rdname illparms
 #' @export
-illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...) {
+illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod = TRUE, slopes = TRUE, ...) {
   if (inherits(object$so, "try-error")) {
     ill_parms <- NA
   } else {
@@ -106,6 +109,15 @@ illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod
       ill_parms_errmod <- ints$errmod[, "lower"] < 0 & ints$errmod[, "est."] > 0
       ill_parms <- c(ill_parms, names(which(ill_parms_errmod)))
     }
+    if (slopes) {
+      if (is.null(object$so)) stop("Slope testing is only implemented for the saemix backend")
+      slope_names <- grep("^beta_", object$so@model@name.fixed, value = TRUE)
+      ci <- object$so@results@conf.int
+      rownames(ci) <- ci$name
+      slope_ci <- ci[slope_names, ]
+      ill_parms_slopes <- slope_ci[, "lower"] < 0 & slope_ci[, "estimate"] > 0
+      ill_parms <- c(ill_parms, slope_names[ill_parms_slopes])
+    }
   }
   class(ill_parms) <- "illparms.saem.mmkin"
   return(ill_parms)
-- 
cgit v1.2.1


From 904ba9668eb76eaae4960e2188134e8c88da07ee Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Wed, 7 Dec 2022 16:19:54 +0100
Subject: Fix parplot for the case of failed multistart runs

---
 R/parms.R   | 2 +-
 R/parplot.R | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

(limited to 'R')

diff --git a/R/parms.R b/R/parms.R
index bd4e479b..bb04a570 100644
--- a/R/parms.R
+++ b/R/parms.R
@@ -77,6 +77,6 @@ parms.multistart <- function(object, exclude_failed = TRUE, ...) {
   successful <- which(!is.na(res[, 1]))
   first_success <- successful[1]
   colnames(res) <- names(parms(object[[first_success]]))
-  if (exclude_failed) res <- res[successful, ]
+  if (exclude_failed[1]) res <- res[successful, ]
   return(res)
 }
diff --git a/R/parplot.R b/R/parplot.R
index 63306ac2..e9c18947 100644
--- a/R/parplot.R
+++ b/R/parplot.R
@@ -41,7 +41,7 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA,
   orig <- attr(object, "orig")
   orig_parms <- parms(orig)
   start_degparms <- orig$mean_dp_start
-  all_parms <- parms(object)
+  all_parms <- parms(object, exclude_failed = FALSE)
 
   if (inherits(object, "multistart.saem.mmkin")) {
     llfunc <- function(object) {
-- 
cgit v1.2.1


From a54bd290bc3884d0000c52c1b29bc557825d9eae Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Thu, 15 Dec 2022 14:50:28 +0100
Subject: List random effects correlations in output if any

Update docs
---
 R/intervals.R          |  8 ++++++--
 R/summary.saem.mmkin.R | 16 ++++++++++++++--
 2 files changed, 20 insertions(+), 4 deletions(-)

(limited to 'R')

diff --git a/R/intervals.R b/R/intervals.R
index 705ef6eb..fcdbaea9 100644
--- a/R/intervals.R
+++ b/R/intervals.R
@@ -78,8 +78,12 @@ intervals.saem.mmkin <- function(object, level = 0.95, backtransform = TRUE, ...
 
   # Random effects
   sdnames <- intersect(rownames(conf.int), paste("SD", pnames, sep = "."))
-  ranef_ret <- as.matrix(conf.int[sdnames, c("lower", "est.", "upper")])
-  rownames(ranef_ret) <- paste0(gsub("SD\\.", "sd(", sdnames), ")")
+  corrnames <- grep("^Corr.", rownames(conf.int), value = TRUE)
+  ranef_ret <- as.matrix(conf.int[c(sdnames, corrnames), c("lower", "est.", "upper")])
+  sdnames_ret <- paste0(gsub("SD\\.", "sd(", sdnames), ")")
+  corrnames_ret <- gsub("Corr\\.(.*)\\.(.*)", "corr(\\1,\\2)", corrnames)
+  rownames(ranef_ret) <- c(sdnames_ret, corrnames_ret)
+
   attr(ranef_ret, "label") <- "Random effects:"
 
 
diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R
index 46ab548b..49b02a50 100644
--- a/R/summary.saem.mmkin.R
+++ b/R/summary.saem.mmkin.R
@@ -75,10 +75,21 @@
 #' f_saem_dfop_sfo <- saem(f_mmkin_dfop_sfo)
 #' print(f_saem_dfop_sfo)
 #' illparms(f_saem_dfop_sfo)
-#' f_saem_dfop_sfo_2 <- update(f_saem_dfop_sfo, covariance.model = diag(c(0, 0, 1, 1, 1, 0)))
+#' f_saem_dfop_sfo_2 <- update(f_saem_dfop_sfo,
+#'   no_random_effect = c("parent_0", "log_k_m1"))
 #' illparms(f_saem_dfop_sfo_2)
 #' intervals(f_saem_dfop_sfo_2)
 #' summary(f_saem_dfop_sfo_2, data = TRUE)
+#' # Add a correlation between random effects of g and k2
+#' cov_model_3 <- f_saem_dfop_sfo_2$so@model@covariance.model
+#' cov_model_3["log_k2", "g_qlogis"] <- 1
+#' cov_model_3["g_qlogis", "log_k2"] <- 1
+#' f_saem_dfop_sfo_3 <- update(f_saem_dfop_sfo,
+#'   covariance.model = cov_model_3)
+#' intervals(f_saem_dfop_sfo_3)
+#' # The correlation does not improve the fit judged by AIC and BIC, although
+#' # the likelihood is higher with the additional parameter
+#' anova(f_saem_dfop_sfo, f_saem_dfop_sfo_2, f_saem_dfop_sfo_3)
 #' }
 #'
 #' @export
@@ -150,7 +161,8 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes =
 
   # Random effects
   sdnames <- intersect(rownames(conf.int), paste0("SD.", pnames))
-  confint_ranef <- as.matrix(conf.int[sdnames, c("estimate", "lower", "upper")])
+  corrnames <- grep("^Corr.", rownames(conf.int), value = TRUE)
+  confint_ranef <- as.matrix(conf.int[c(sdnames, corrnames), c("estimate", "lower", "upper")])
   colnames(confint_ranef)[1] <- "est."
 
   # Error model
-- 
cgit v1.2.1


From 886c9ef013124aa954d960c655b349b5340ff154 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Mon, 19 Dec 2022 12:31:56 +0100
Subject: Rename template folder, create format

Instead of rmarkdown::pdf_document, mkin::hierarchical_kinetics is used
as a document format in the template. In this way, the template file can
be freed from some R code and yaml options that the average user does
not have to be aware of.
---
 R/hierarchical_kinetics.R | 39 +++++++++++++++++++++++++++++++++++++++
 R/parplot.R               |  2 +-
 2 files changed, 40 insertions(+), 1 deletion(-)
 create mode 100644 R/hierarchical_kinetics.R

(limited to 'R')

diff --git a/R/hierarchical_kinetics.R b/R/hierarchical_kinetics.R
new file mode 100644
index 00000000..f7ffb333
--- /dev/null
+++ b/R/hierarchical_kinetics.R
@@ -0,0 +1,39 @@
+#' Hierarchical kinetics template
+#'
+#' R markdown format for setting up hierarchical kinetics based on a template
+#' provided with the mkin package.
+#'
+#' @inheritParams rmarkdown::pdf_document
+#' @param ... Arguments to \code{rmarkdown::pdf_document}
+#'
+#' @return R Markdown output format to pass to
+#'   \code{\link[rmarkdown:render]{render}}
+#'
+#' @examples
+#'
+#' \dontrun{
+#' library(rmarkdown)
+#' draft("New analysis.rmd", template = "hierarchical_kinetics", package = "mkin")
+#' }
+#'
+#' @export
+hierarchical_kinetics <- function(..., keep_tex = FALSE) {
+
+  if (getRversion() < "4.1.0")
+    stop("You need R with version > 4.1.0 to compile this document")
+
+  if (!requireNamespace("knitr")) stop("Please install the knitr package to use this template")
+  if (!requireNamespace("rmarkdown")) stop("Please install the rmarkdown package to use this template")
+  knitr::opts_chunk$set(echo = FALSE, cache = TRUE, comment = "", tidy = FALSE)
+  knitr::opts_chunk$set(fig.align = "center", fig.pos = "H")
+  options(knitr.kable.NA = "")
+
+  fmt <- rmarkdown::pdf_document(...,
+    keep_tex = keep_tex,
+    toc = TRUE,
+    includes = rmarkdown::includes(in_header = "header.tex"),
+    extra_dependencies = c("float", "listing", "framed")
+  )
+
+  return(fmt)
+}
diff --git a/R/parplot.R b/R/parplot.R
index e9c18947..3da4b51a 100644
--- a/R/parplot.R
+++ b/R/parplot.R
@@ -23,7 +23,7 @@
 #' of the in vitro erythropoiesis. BMC Bioinformatics. 2021 Oct 4;22(1):478.
 #' doi: 10.1186/s12859-021-04373-4.
 #' @seealso [multistart]
-#' @importFrom stats median
+#' @importFrom stats median quantile
 #' @export
 parplot <- function(object, ...) {
   UseMethod("parplot")
-- 
cgit v1.2.1


From a89a835e0fd6fe55e4c7353e99cc442c4d309bc1 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Mon, 2 Jan 2023 10:25:47 +0100
Subject: Echo R code per default in markdown template

---
 R/hierarchical_kinetics.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

(limited to 'R')

diff --git a/R/hierarchical_kinetics.R b/R/hierarchical_kinetics.R
index f7ffb333..70e0100a 100644
--- a/R/hierarchical_kinetics.R
+++ b/R/hierarchical_kinetics.R
@@ -24,7 +24,7 @@ hierarchical_kinetics <- function(..., keep_tex = FALSE) {
 
   if (!requireNamespace("knitr")) stop("Please install the knitr package to use this template")
   if (!requireNamespace("rmarkdown")) stop("Please install the rmarkdown package to use this template")
-  knitr::opts_chunk$set(echo = FALSE, cache = TRUE, comment = "", tidy = FALSE)
+  knitr::opts_chunk$set(echo = FALSE, cache = TRUE, comment = "", tidy = FALSE, echo = TRUE)
   knitr::opts_chunk$set(fig.align = "center", fig.pos = "H")
   options(knitr.kable.NA = "")
 
-- 
cgit v1.2.1


From 435b7bc8fcbb55b5f487cf5d92d60833daf75ae7 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Mon, 2 Jan 2023 10:26:40 +0100
Subject: Fix no_random_effect with character vector

---
 R/mhmkin.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

(limited to 'R')

diff --git a/R/mhmkin.R b/R/mhmkin.R
index a9798ff4..6265a59e 100644
--- a/R/mhmkin.R
+++ b/R/mhmkin.R
@@ -134,7 +134,7 @@ mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
   dimnames(fit_indices) <- list(degradation = names(deg_models),
                                 error = error_models)
 
-  if (is.null(no_random_effect) || length(dim(no_random_effect)) == 1) {
+  if (is.null(no_random_effect) || is.null(dim(no_random_effect))) {
     no_ranef <- rep(list(no_random_effect), n.fits)
     dim(no_ranef) <- dim(fit_indices)
   } else {
-- 
cgit v1.2.1


From ed2ab7bba07e3cb036782227fe27943f4b5583fa Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Tue, 3 Jan 2023 09:31:48 +0100
Subject: Improved skeleton for hierarchical fits

Now with working pathway fits using SFORB-SFO2 (only two parallel
metabolites instead of three) as the data for compound Ia was not
sufficient for a reliable fit.
---
 R/hierarchical_kinetics.R | 1 +
 1 file changed, 1 insertion(+)

(limited to 'R')

diff --git a/R/hierarchical_kinetics.R b/R/hierarchical_kinetics.R
index 70e0100a..a90dd619 100644
--- a/R/hierarchical_kinetics.R
+++ b/R/hierarchical_kinetics.R
@@ -31,6 +31,7 @@ hierarchical_kinetics <- function(..., keep_tex = FALSE) {
   fmt <- rmarkdown::pdf_document(...,
     keep_tex = keep_tex,
     toc = TRUE,
+    toc_depth = 3,
     includes = rmarkdown::includes(in_header = "header.tex"),
     extra_dependencies = c("float", "listing", "framed")
   )
-- 
cgit v1.2.1


From a583cd7e3eecef4c70ac09304495c765a84e76ce Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Wed, 4 Jan 2023 05:34:42 +0100
Subject: Update documentation of 'use_of_ff' argument

---
 R/mkinmod.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

(limited to 'R')

diff --git a/R/mkinmod.R b/R/mkinmod.R
index d8740aed..47307ab7 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -26,7 +26,7 @@
 #'   Additionally, [mkinsub()] has an argument \code{to}, specifying names of
 #'   variables to which a transfer is to be assumed in the model.
 #'   If the argument \code{use_of_ff} is set to "min"
-#'   (default) and the model for the compartment is "SFO" or "SFORB", an
+#'   and the model for the compartment is "SFO" or "SFORB", an
 #'   additional [mkinsub()] argument can be \code{sink = FALSE}, effectively
 #'   fixing the flux to sink to zero.
 #'   In print.mkinmod, this argument is currently not used.
-- 
cgit v1.2.1


From a5903e74d9cf54c764d5bbc48e461cecd5f56e72 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Thu, 5 Jan 2023 15:03:06 +0100
Subject: Don't preschedule multistart runs

Sometimes a lot of them fail, so we were wasting time
---
 R/multistart.R | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

(limited to 'R')

diff --git a/R/multistart.R b/R/multistart.R
index 61ef43dc..bdfbfe63 100644
--- a/R/multistart.R
+++ b/R/multistart.R
@@ -100,9 +100,12 @@ multistart.saem.mmkin <- function(object, n = 50, cores = 1,
   }
 
   if (is.null(cluster)) {
-    res <- parallel::mclapply(1:n, fit_function, mc.cores = cores)
+    res <- parallel::mclapply(1:n, fit_function,
+      mc.cores = cores, mc.preschedule = FALSE)
   } else {
-    res <- parallel::parLapply(cluster, 1:n, fit_function)
+    res <- parallel::parLapplyLB(cluster, 1:n, fit_function,
+      mc.preschedule = FALSE
+    )
   }
   attr(res, "orig") <- object
   attr(res, "start_parms") <- start_parms
-- 
cgit v1.2.1


From 24eb77216700cf8b2f2bde3abad84c1f83f9e32a Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Mon, 9 Jan 2023 06:22:04 +0100
Subject: Prebuilt PDF vignettes, summary_listing

---
 R/hierarchical_kinetics.R |  2 +-
 R/summary_listing.R       | 59 +++++++++++++++++++++++++++++++++++++++++++++++
 R/tex_listing.R           | 32 -------------------------
 3 files changed, 60 insertions(+), 33 deletions(-)
 create mode 100644 R/summary_listing.R
 delete mode 100644 R/tex_listing.R

(limited to 'R')

diff --git a/R/hierarchical_kinetics.R b/R/hierarchical_kinetics.R
index a90dd619..e545754a 100644
--- a/R/hierarchical_kinetics.R
+++ b/R/hierarchical_kinetics.R
@@ -13,7 +13,7 @@
 #'
 #' \dontrun{
 #' library(rmarkdown)
-#' draft("New analysis.rmd", template = "hierarchical_kinetics", package = "mkin")
+#' draft("example_analysis.rmd", template = "hierarchical_kinetics", package = "mkin")
 #' }
 #'
 #' @export
diff --git a/R/summary_listing.R b/R/summary_listing.R
new file mode 100644
index 00000000..38b52240
--- /dev/null
+++ b/R/summary_listing.R
@@ -0,0 +1,59 @@
+#' Display the output of a summary function according to the output format
+#'
+#' This function is intended for use in a R markdown code chunk with the chunk
+#' option `results = "asis"`.
+#'
+#' @param object The object for which the summary is to be listed
+#' @param caption An optional caption
+#' @param label An optional label, ignored in html output
+#' @param clearpage Should a new page be started after the listing? Ignored in html output
+#' @export
+summary_listing <- function(object, caption = NULL, label = NULL,
+  clearpage = TRUE) {
+  if (knitr::is_latex_output()) {
+    tex_listing(object = object, caption = caption, label = label,
+      clearpage = clearpage)
+  }
+  if (knitr::is_html_output()) {
+    html_listing(object = object, caption = caption)
+  }
+}
+
+#' @rdname summary_listing
+#' @export
+tex_listing <- function(object, caption = NULL, label = NULL,
+  clearpage = TRUE) {
+  cat("\n")
+  cat("\\begin{listing}", "\n")
+  if (!is.null(caption)) {
+    cat("\\caption{", caption, "}", "\n", sep = "")
+  }
+  if (!is.null(label)) {
+    cat("\\caption{", label, "}", "\n", sep = "")
+  }
+  cat("\\begin{snugshade}", "\n")
+  cat("\\scriptsize", "\n")
+  cat("\\begin{verbatim}", "\n")
+  cat(capture.output(suppressWarnings(summary(object))), sep = "\n")
+  cat("\n")
+  cat("\\end{verbatim}", "\n")
+  cat("\\end{snugshade}", "\n")
+  cat("\\end{listing}", "\n")
+  if (clearpage) {
+    cat("\\clearpage", "\n")
+  }
+}
+
+#' @rdname summary_listing
+#' @export
+html_listing <- function(object, caption = NULL) {
+  cat("\n")
+  if (!is.null(caption)) {
+    cat("<caption>", caption, "</caption>", "\n", sep = "")
+  }
+  cat("<pre><code>\n")
+  cat(capture.output(suppressWarnings(summary(object))), sep = "\n")
+  cat("\n")
+  cat("</pre></code>\n")
+}
+
diff --git a/R/tex_listing.R b/R/tex_listing.R
deleted file mode 100644
index 05f662e4..00000000
--- a/R/tex_listing.R
+++ /dev/null
@@ -1,32 +0,0 @@
-#' Wrap the output of a summary function in tex listing environment
-#'
-#' This function can be used in a R markdown code chunk with the chunk
-#' option `results = "asis"`.
-#'
-#' @param object The object for which the summary is to be listed
-#' @param caption An optional caption
-#' @param label An optional label
-#' @param clearpage Should a new page be started after the listing?
-#' @export
-tex_listing <- function(object, caption = NULL, label = NULL,
-  clearpage = TRUE) {
-  cat("\n")
-  cat("\\begin{listing}", "\n")
-  if (!is.null(caption)) {
-    cat("\\caption{", caption, "}", "\n", sep = "")
-  }
-  if (!is.null(label)) {
-    cat("\\caption{", label, "}", "\n", sep = "")
-  }
-  cat("\\begin{snugshade}", "\n")
-  cat("\\scriptsize", "\n")
-  cat("\\begin{verbatim}", "\n")
-  cat(capture.output(suppressWarnings(summary(object))), sep = "\n")
-  cat("\n")
-  cat("\\end{verbatim}", "\n")
-  cat("\\end{snugshade}", "\n")
-  cat("\\end{listing}", "\n")
-  if (clearpage) {
-    cat("\\clearpage", "\n")
-  }
-}
-- 
cgit v1.2.1