Use all cores per default, confint tolerance

Also, use more intelligent starting values for the variance of the
random effects for saemix. While this does not appear to speed up
the convergence, it shows where this variance is greatly reduced
by using mixed-effects models as opposed to the separate independent
fits.

5 files changed, 78 insertions, 50 deletions

diff --git a/R/confint.mkinfit.R b/R/confint.mkinfit.R
index 78dda95d..53eb45ee 100644
--- a/R/confint.mkinfit.R
+++ b/R/confint.mkinfit.R
@@ -27,7 +27,10 @@
 #' @param backtransform If we approximate the likelihood in terms of the
 #'   transformed parameters, should we backtransform the parameters with
 #'   their confidence intervals?
-#' @param cores The number of cores to be used for multicore processing. 
+#' @param rel_tol If the method is 'profile', what should be the accuracy
+#'   of the lower and upper bounds, relative to the estimate obtained from
+#'   the quadratic method?
+#' @param cores The number of cores to be used for multicore processing.
 #'   On Windows machines, cores > 1 is currently not supported.
 #' @param quiet Should we suppress the message "Profiling the likelihood"
 #' @return A matrix with columns giving lower and upper confidence limits for
@@ -121,7 +124,7 @@ confint.mkinfit <- function(object, parm,
   level = 0.95, alpha = 1 - level, cutoff,
   method = c("quadratic", "profile"),
   transformed = TRUE, backtransform = TRUE,
-  cores = round(detectCores()/2), quiet = FALSE, ...)
+  cores = parallel::detectCores(), rel_tol = 0.01, quiet = FALSE, ...)
 {
   tparms <- parms(object, transformed = TRUE)
   bparms <- parms(object, transformed = FALSE)
@@ -140,50 +143,50 @@ confint.mkinfit <- function(object, parm,
 
   a <- c(alpha / 2, 1 - (alpha / 2))
 
-  if (method == "quadratic") {
+  quantiles <- qt(a, object$df.residual)
 
-    quantiles <- qt(a, object$df.residual)
-
-    covar_pnames <- if (missing(parm)) {
-      if (transformed) tpnames else bpnames
-    } else {
-      parm
-    }
+  covar_pnames <- if (missing(parm)) {
+    if (transformed) tpnames else bpnames
+  } else {
+    parm
+  }
 
-    return_parms <- if (backtransform) bparms[return_pnames]
-      else tparms[return_pnames]
+  return_parms <- if (backtransform) bparms[return_pnames]
+    else tparms[return_pnames]
 
-    covar_parms <- if (transformed) tparms[covar_pnames]
-      else bparms[covar_pnames]
+  covar_parms <- if (transformed) tparms[covar_pnames]
+    else bparms[covar_pnames]
 
-    if (transformed) {
-      covar <- try(solve(object$hessian), silent = TRUE)
-    } else {
-      covar <- try(solve(object$hessian_notrans), silent = TRUE)
-    }
+  if (transformed) {
+    covar <- try(solve(object$hessian), silent = TRUE)
+  } else {
+    covar <- try(solve(object$hessian_notrans), silent = TRUE)
+  }
 
-    # If inverting the covariance matrix failed or produced NA values
-    if (!is.numeric(covar) | is.na(covar[1])) {
-      ses <- lci <- uci <- rep(NA, p)
-    } else {
-      ses     <- sqrt(diag(covar))[covar_pnames]
-      lci    <- covar_parms + quantiles[1] * ses
-      uci    <- covar_parms + quantiles[2] * ses
-      if (transformed & backtransform) {
-        lci_back <- backtransform_odeparms(lci,
-          object$mkinmod, object$transform_rates, object$transform_fractions)
-        uci_back <- backtransform_odeparms(uci,
-          object$mkinmod, object$transform_rates, object$transform_fractions)
-
-        return_errparm_names <- intersect(names(object$errparms), return_pnames)
-        lci <- c(lci_back, lci[return_errparm_names])
-        uci <- c(uci_back, uci[return_errparm_names])
-      }
+  # If inverting the covariance matrix failed or produced NA values
+  if (!is.numeric(covar) | is.na(covar[1])) {
+    ses <- lci <- uci <- rep(NA, p)
+  } else {
+    ses     <- sqrt(diag(covar))[covar_pnames]
+    lci    <- covar_parms + quantiles[1] * ses
+    uci    <- covar_parms + quantiles[2] * ses
+    if (transformed & backtransform) {
+      lci_back <- backtransform_odeparms(lci,
+        object$mkinmod, object$transform_rates, object$transform_fractions)
+      uci_back <- backtransform_odeparms(uci,
+        object$mkinmod, object$transform_rates, object$transform_fractions)
+
+      return_errparm_names <- intersect(names(object$errparms), return_pnames)
+      lci <- c(lci_back, lci[return_errparm_names])
+      uci <- c(uci_back, uci[return_errparm_names])
     }
-    ci <- cbind(lower = lci, upper = uci)
   }
+  ci <- cbind(lower = lci, upper = uci)
 
   if (method == "profile") {
+
+    ci_quadratic <- ci
+
     if (!quiet) message("Profiling the likelihood")
 
     lci <- uci <- rep(NA, p)
@@ -215,9 +218,14 @@ confint.mkinfit <- function(object, parm,
         (cutoff - (object$logLik - profile_ll(x)))^2
       }
 
-      lci_pname <- optimize(cost, lower = 0, upper = all_parms[pname])$minimum
+      lower_quadratic <- ci_quadratic["lower"][pname]
+      upper_quadratic <- ci_quadratic["upper"][pname]
+      ltol <- if (!is.na(lower_quadratic)) rel_tol * lower_quadratic else .Machine$double.eps^0.25
+      utol <- if (!is.na(upper_quadratic)) rel_tol * upper_quadratic else .Machine$double.eps^0.25
+      lci_pname <- optimize(cost, lower = 0, upper = all_parms[pname], tol = ltol)$minimum
       uci_pname <- optimize(cost, lower = all_parms[pname],
-        upper = ifelse(grepl("^f_|^g$", pname), 1, 15 * all_parms[pname]))$minimum
+        upper = ifelse(grepl("^f_|^g$", pname), 1, 15 * all_parms[pname]),
+        tol = utol)$minimum
       return(c(lci_pname, uci_pname))
     }
     ci <- t(parallel::mcmapply(get_ci, profile_pnames, mc.cores = cores))
diff --git a/R/mmkin.R b/R/mmkin.R
index dbb61b78..37c4e87d 100644
--- a/R/mmkin.R
+++ b/R/mmkin.R
@@ -12,7 +12,8 @@
 #' @param cores The number of cores to be used for multicore processing. This
 #'   is only used when the \code{cluster} argument is \code{NULL}. On Windows
 #'   machines, cores > 1 is not supported, you need to use the \code{cluster}
-#'   argument to use multiple logical processors.
+#'   argument to use multiple logical processors. Per default, all cores
+#'   detected by [parallel::detectCores()] are used.
 #' @param cluster A cluster as returned by \code{\link{makeCluster}} to be used
 #'   for parallel execution.
 #' @param \dots Further arguments that will be passed to \code{\link{mkinfit}}.
@@ -62,7 +63,7 @@
 #'
 #' @export mmkin
 mmkin <- function(models = c("SFO", "FOMC", "DFOP"), datasets,
-                  cores = round(detectCores()/2), cluster = NULL, ...)
+  cores = detectCores(), cluster = NULL, ...)
 {
   parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE", "logistic")
   n.m <- length(models)
diff --git a/R/nlme.R b/R/nlme.R
index 3ee7b9fd..20987064 100644
--- a/R/nlme.R
+++ b/R/nlme.R
@@ -125,7 +125,7 @@ nlme_function <- function(object) {
 #' @return If random is FALSE (default), a named vector containing mean values
 #'   of the fitted degradation model parameters. If random is TRUE, a list with
 #'   fixed and random effects, in the format required by the start argument of
-#'   nlme for the case of a single grouping variable ds?
+#'   nlme for the case of a single grouping variable ds.
 #' @param random Should a list with fixed and random effects be returned?
 #' @export
 mean_degparms <- function(object, random = FALSE) {
diff --git a/R/parms.mkinfit.R b/R/parms.mkinfit.R
index aae6fa52..a1f2d209 100644
--- a/R/parms.mkinfit.R
+++ b/R/parms.mkinfit.R
@@ -21,11 +21,13 @@
 #' ds <- lapply(experimental_data_for_UBA_2019[6:10],
 #'  function(x) subset(x$data[c("name", "time", "value")]))
 #' names(ds) <- paste("Dataset", 6:10)
-#' fits <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE)
+#' \dontrun{
+#' fits <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE, cores = 1)
 #' parms(fits["SFO", ])
 #' parms(fits[, 2])
 #' parms(fits)
 #' parms(fits, transformed = TRUE)
+#' }
 #' @export
 parms <- function(object, ...)
 {
diff --git a/R/saemix.R b/R/saemix.R
index 69e5fc50..24c0f0d0 100644
--- a/R/saemix.R
+++ b/R/saemix.R
@@ -5,25 +5,37 @@
 #' list of mkinfit objects that have been obtained by fitting the same model to
 #' a list of datasets.
 #'
+#' Starting values for the fixed effects (population mean parameters, argument psi0 of
+#' [saemix::saemixModel()] are the mean values of the parameters found using
+#' mmkin. Starting variances of the random effects (argument omega.init) are the
+#' variances of the deviations of the parameters from these mean values.
+#'
 #' @param object An mmkin row object containing several fits of the same model to different datasets
+#' @param cores The number of cores to be used for multicore processing.
+#'   On Windows machines, cores > 1 is currently not supported.
 #' @rdname saemix
 #' @importFrom saemix saemixData saemixModel
+#' @importFrom stats var
 #' @examples
 #' ds <- lapply(experimental_data_for_UBA_2019[6:10],
 #'  function(x) subset(x$data[c("name", "time", "value")]))
 #' names(ds) <- paste("Dataset", 6:10)
 #' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
 #'   A1 = mkinsub("SFO"))
-#' f_mmkin <- mmkin(list("SFO-SFO" = sfo_sfo), ds, quiet = TRUE, cores = 5)
+#' \dontrun{
+#' f_mmkin <- mmkin(list("SFO-SFO" = sfo_sfo), ds, quiet = TRUE)
+#' library(saemix)
 #' m_saemix <- saemix_model(f_mmkin)
 #' d_saemix <- saemix_data(f_mmkin)
-#' saemix_options <- list(seed = 123456, save = FALSE, save.graphs = FALSE)
-#' \dontrun{
-#'   saemix(m_saemix, d_saemix, saemix_options)
+#' saemix_options <- list(seed = 123456,
+#'   save = FALSE, save.graphs = FALSE, displayProgress = FALSE,
+#'   nbiter.saemix = c(200, 80))
+#' f_saemix <- saemix(m_saemix, d_saemix, saemix_options)
+#' plot(f_saemix, plot.type = "convergence")
 #' }
 #' @return An [saemix::SaemixModel] object.
 #' @export
-saemix_model <- function(object) {
+saemix_model <- function(object, cores = parallel::detectCores()) {
   if (nrow(object) > 1) stop("Only row objects allowed")
 
   mkin_model <- object[[1]]$mkinmod
@@ -81,14 +93,19 @@ saemix_model <- function(object) {
           out_values <- out_wide[out_index]
         }
         return(out_values)
-      }, mc.cores = 15)
+      }, mc.cores = cores)
       res <- unlist(res_list)
       return(res)
   }
 
+  raneff_0 <- mean_degparms(object, random = TRUE)$random$ds
+  var_raneff_0 <- apply(raneff_0, 2, var)
+
   res <- saemixModel(model_function, psi0,
     "Mixed model generated from mmkin object",
-    transform.par = rep(0, length(degparms_optim)))
+    transform.par = rep(0, length(degparms_optim)),
+    omega.init = diag(var_raneff_0)
+  )
   return(res)
 }