Some improvements for borderline cases

- fit_with_errors for saem()
- test_log_parms for mean_degparms() and saem()

3 files changed, 61 insertions, 9 deletions

diff --git a/R/nlme.R b/R/nlme.R
index 9215aab0..d235a094 100644
--- a/R/nlme.R
+++ b/R/nlme.R
@@ -36,7 +36,7 @@
 #' nlme_f <- nlme_function(f)
 #' # These assignments are necessary for these objects to be
 #' # visible to nlme and augPred when evaluation is done by
-#' # pkgdown to generated the html docs.
+#' # pkgdown to generate the html docs.
 #' assign("nlme_f", nlme_f, globalenv())
 #' assign("grouped_data", grouped_data, globalenv())
 #'
@@ -130,13 +130,44 @@ nlme_function <- function(object) {
 #'   fixed and random effects, in the format required by the start argument of
 #'   nlme for the case of a single grouping variable ds.
 #' @param random Should a list with fixed and random effects be returned?
+#' @param test_log_parms If TRUE, log parameters are only considered in
+#'   the mean calculations if their untransformed counterparts (most likely
+#'   rate constants) pass the t-test for significant difference from zero.
+#' @param conf.level Possibility to adjust the required confidence level
+#'   for parameter that are tested if requested by 'test_log_parms'.
 #' @export
-mean_degparms <- function(object, random = FALSE) {
+mean_degparms <- function(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6)
+{
   if (nrow(object) > 1) stop("Only row objects allowed")
   parm_mat_trans <- sapply(object, parms, transformed = TRUE)
+
+  if (test_log_parms) {
+      parm_mat_dim <- dim(parm_mat_trans)
+      parm_mat_dimnames <- dimnames(parm_mat_trans)
+
+      log_parm_trans_names <- grep("^log_", rownames(parm_mat_trans), value = TRUE)
+      log_parm_names <- gsub("^log_", "", log_parm_trans_names)
+
+      t_test_back_OK <- matrix(
+        sapply(object, function(o) {
+          suppressWarnings(summary(o)$bpar[log_parm_names, "Pr(>t)"] < (1 - conf.level))
+        }), nrow = length(log_parm_names))
+      rownames(t_test_back_OK) <- log_parm_trans_names
+
+      parm_mat_trans_OK <- parm_mat_trans
+      for (trans_parm in log_parm_trans_names) {
+        parm_mat_trans_OK[trans_parm, ] <- ifelse(t_test_back_OK[trans_parm, ],
+          parm_mat_trans[trans_parm, ], NA)
+      }
+    } else {
+    parm_mat_trans_OK <- parm_mat_trans
+  }
+
   mean_degparm_names <- setdiff(rownames(parm_mat_trans), names(object[[1]]$errparms))
   degparm_mat_trans <- parm_mat_trans[mean_degparm_names, , drop = FALSE]
-  fixed <- apply(degparm_mat_trans, 1, mean)
+  degparm_mat_trans_OK <- parm_mat_trans_OK[mean_degparm_names, , drop = FALSE]
+
+  fixed <- apply(degparm_mat_trans_OK, 1, mean, na.rm = TRUE)
   if (random) {
     random <- t(apply(degparm_mat_trans[mean_degparm_names, , drop = FALSE], 2, function(column) column - fixed))
     # If we only have one parameter, apply returns a vector so we get a single row
diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R
index ff1f2fff..306600c6 100644
--- a/R/nlme.mmkin.R
+++ b/R/nlme.mmkin.R
@@ -24,7 +24,7 @@ get_deg_func <- function() {
 #' This functions sets up a nonlinear mixed effects model for an mmkin row
 #' object. An mmkin row object is essentially a list of mkinfit objects that
 #' have been obtained by fitting the same model to a list of datasets.
-#' 
+#'
 #' Note that the convergence of the nlme algorithms depends on the quality
 #' of the data. In degradation kinetics, we often only have few datasets
 #' (e.g. data for few soils) and complicated degradation models, which may
diff --git a/R/saem.R b/R/saem.R
index fd2a77b4..460edede 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -24,8 +24,16 @@ utils::globalVariables(c("predicted", "std"))
 #'   SFO or DFOP is used for the parent and there is either no metabolite or one.
 #' @param degparms_start Parameter values given as a named numeric vector will
 #'   be used to override the starting values obtained from the 'mmkin' object.
+#' @param test_log_parms If TRUE, an attempt is made to use more robust starting
+#'   values for population parameters fitted as log parameters in mkin (like
+#'   rate constants) by only considering rate constants that pass the t-test
+#'   when calculating mean degradation parameters using [mean_degparms].
+#' @param conf.level Possibility to adjust the required confidence level
+#'   for parameter that are tested if requested by 'test_log_parms'.
 #' @param solution_type Possibility to specify the solution type in case the
 #'   automatic choice is not desired
+#' @param fail_with_errors Should a failure to compute standard errors
+#'   from the inverse of the Fisher Information Matrix be a failure?
 #' @param quiet Should we suppress the messages saemix prints at the beginning
 #'   and the end of the optimisation process?
 #' @param control Passed to [saemix::saemix]
@@ -51,7 +59,7 @@ utils::globalVariables(c("predicted", "std"))
 #' # The returned saem.mmkin object contains an SaemixObject, therefore we can use
 #' # functions from saemix
 #' library(saemix)
-#' compare.saemix(list(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so))
+#' compare.saemix(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so)
 #' plot(f_saem_fomc$so, plot.type = "convergence")
 #' plot(f_saem_fomc$so, plot.type = "individual.fit")
 #' plot(f_saem_fomc$so, plot.type = "npde")
@@ -59,7 +67,7 @@ utils::globalVariables(c("predicted", "std"))
 #'
 #' f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc")
 #' f_saem_fomc_tc <- saem(f_mmkin_parent_tc["FOMC", ])
-#' compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so))
+#' compare.saemix(f_saem_fomc$so, f_saem_fomc_tc$so)
 #'
 #' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
 #'   A1 = mkinsub("SFO"))
@@ -104,19 +112,32 @@ saem <- function(object, ...) UseMethod("saem")
 saem.mmkin <- function(object,
   transformations = c("mkin", "saemix"),
   degparms_start = numeric(),
+  test_log_parms = FALSE,
+  conf.level = 0.6,
   solution_type = "auto",
   control = list(displayProgress = FALSE, print = FALSE,
     save = FALSE, save.graphs = FALSE),
+  fail_with_errors = TRUE,
   verbose = FALSE, quiet = FALSE, ...)
 {
   transformations <- match.arg(transformations)
   m_saemix <- saemix_model(object, verbose = verbose,
-    degparms_start = degparms_start, solution_type = solution_type,
+    degparms_start = degparms_start,
+    test_log_parms = test_log_parms, conf.level = conf.level,
+    solution_type = solution_type,
     transformations = transformations, ...)
   d_saemix <- saemix_data(object, verbose = verbose)
 
   fit_time <- system.time({
     utils::capture.output(f_saemix <- saemix::saemix(m_saemix, d_saemix, control), split = !quiet)
+    FIM_failed <- NULL
+    if (any(is.na(f_saemix@results@se.fixed))) FIM_failed <- c(FIM_failed, "fixed effects")
+    if (any(is.na(c(f_saemix@results@se.omega, f_saemix@results@se.respar)))) {
+      FIM_failed <- c(FIM_failed, "random effects and residual error parameters")
+    }
+    if (!is.null(FIM_failed) & fail_with_errors) {
+      stop("Could not invert FIM for ", paste(FIM_failed, collapse = " and "))
+    }
   })
 
   transparms_optim <- f_saemix@results@fixed.effects
@@ -203,13 +224,13 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
 #' @return An [saemix::SaemixModel] object.
 #' @export
 saemix_model <- function(object, solution_type = "auto", transformations = c("mkin", "saemix"),
-  degparms_start = numeric(), verbose = FALSE, ...)
+  degparms_start = numeric(), test_log_parms = FALSE, verbose = FALSE, ...)
 {
   if (nrow(object) > 1) stop("Only row objects allowed")
 
   mkin_model <- object[[1]]$mkinmod
 
-  degparms_optim <-  mean_degparms(object)
+  degparms_optim <-  mean_degparms(object, test_log_parms = test_log_parms)
   if (transformations == "saemix") {
     degparms_optim <- backtransform_odeparms(degparms_optim,
       object[[1]]$mkinmod,