Batch processing for hierarchical fits

- 'R/mhmkin.R': New method for performing multiple hierarchical mkin fits in one function call, optionally in parallel.

- 'R/saem.R': 'logLik' and 'update' methods for 'saem.mmkin' objects.

- 'R/illparms.R': Add methods for 'saem.mmkin' and 'mhmkin' objects.

tests: Use 2 cores on travis, should work according to docs

3 files changed, 285 insertions, 3 deletions

diff --git a/R/illparms.R b/R/illparms.R
index f23f1cae..e4b28c56 100644
--- a/R/illparms.R
+++ b/R/illparms.R
@@ -1,12 +1,29 @@
 #' Method to get the names of ill-defined parameters
 #'
+#' The method for generalised nonlinear regression fits as obtained
+#' with [mkinfit] and [mmkin] checks if the degradation parameters
+#' pass the Wald test (in degradation kinetics often simply called t-test) for
+#' significant difference from zero. For this test, the parameterisation
+#' without parameter transformations is used.
+#'
+#' The method for hierarchical model fits, also known as nonlinear
+#' mixed-effects model fits as obtained with [saem] and [mhmkin]
+#' checks if any of the confidence intervals for the random
+#' effects expressed as standard deviations include zero, and if
+#' the confidence intervals for the error model parameters include
+#' zero.
+#'
 #' @param object The object to investigate
 #' @param x The object to be printed
 #' @param conf.level The confidence level for checking p values
 #' @param \dots For potential future extensions
-#' @return For [mkinfit] objects, a character vector of parameter names
-#' For [mmkin] objects, an object of class 'illparms.mmkin' with a
-#' suitable printing method.
+#' @param random For hierarchical fits, should random effects be tested?
+#' @param errmod For hierarchical fits, should error model parameters 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'. The latter objects have a suitable
+#' printing method.
 #' @export
 illparms <- function(object, ...)
 {
@@ -60,3 +77,50 @@ print.illparms.mmkin <- function(x, ...) {
   class(x) <- NULL
   print(x, quote = FALSE)
 }
+
+#' @rdname illparms
+#' @export
+illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...) {
+  if (inherits(object, "try-error")) {
+    return(NA)
+  } else {
+    ints <- intervals(object, conf.level = conf.level)
+    failed <- NULL
+    if (random) {
+      failed_random <- ints$random[, "lower"] < 0
+      failed <- c(failed, names(which(failed_random)))
+    }
+    if (errmod) {
+      failed_errmod <- ints$errmod[, "lower"] < 0 & ints$errmod[, "est."] > 0
+      failed <- c(failed, names(which(failed_errmod)))
+    }
+  }
+  return(failed)
+}
+
+#' @rdname illparms
+#' @export
+illparms.mhmkin <- function(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...) {
+  result <- lapply(object,
+    function(fit) {
+      if (inherits(fit, "try-error")) return("E")
+      ill <- illparms(fit, conf.level = conf.level, random = random, errmod = errmod)
+      if (length(ill) > 0) {
+        return(paste(ill, collapse = ", "))
+      } else {
+        return("")
+      }
+    })
+  result <- unlist(result)
+  dim(result) <- dim(object)
+  dimnames(result) <- dimnames(object)
+  class(result) <- "illparms.mhmkin"
+  return(result)
+}
+
+#' @rdname illparms
+#' @export
+print.illparms.mhmkin <- function(x, ...) {
+  class(x) <- NULL
+  print(x, quote = FALSE)
+}
diff --git a/R/mhmkin.R b/R/mhmkin.R
new file mode 100644
index 00000000..61d67204
--- /dev/null
+++ b/R/mhmkin.R
@@ -0,0 +1,184 @@
+#' Fit nonlinear mixed-effects models built from one or more kinetic
+#' degradation models and one or more error models
+#'
+#' The name of the methods expresses that (**m**ultiple) **h**ierarchichal
+#' (also known as multilevel) **m**ulticompartment **kin**etic models are
+#' fitted. Our kinetic models are nonlinear, so we can use various nonlinear
+#' mixed-effects model fitting functions.
+#'
+#' @param objects A list of [mmkin] objects containing fits of the same
+#' degradation models to the same data, but using different error models.
+#' @param backend The backend to be used for fitting. Currently, only saemix is
+#' supported
+#' @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.
+#' @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. Per default, all cores detected
+#' by [parallel::detectCores()] are used, except on Windows where the default
+#' is 1.
+#' @param cluster A cluster as returned by [makeCluster] to be used for
+#' parallel execution.
+#' @importFrom parallel mclapply parLapply detectCores
+#' @return A two-dimensional [array] of fit objects and/or try-errors that can
+#' be indexed using the degradation model names for the first index (row index)
+#' and the error model names for the second index (column index), with class
+#' attribute 'mhmkin'.
+#' @author Johannes Ranke
+#' @seealso \code{\link{[.mhmkin}} for subsetting [mhmkin] objects
+#' @export
+mhmkin <- function(objects, backend = "saemix", algorithm = "saem", ...) {
+  UseMethod("mhmkin")
+}
+
+#' @export
+#' @rdname mhmkin
+mhmkin.list <- function(objects, backend = "saemix",
+  ...,
+  cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(), cluster = NULL)
+{
+  call <- match.call()
+  dot_args <- list(...)
+  backend_function <- switch(backend,
+    saemix = "saem"
+  )
+
+  deg_models <- lapply(objects[[1]][, 1], function(x) x$mkinmod)
+  names(deg_models) <- dimnames(objects[[1]])$model
+  n.deg <- length(deg_models)
+
+  ds <- lapply(objects[[1]][1, ], function(x) x$data)
+
+  for (other in objects[-1]) {
+    # Check if the degradation models in all objects are the same
+    for (deg_model_name in names(deg_models)) {
+      if (!all.equal(other[[deg_model_name, 1]]$mkinmod$spec,
+            deg_models[[deg_model_name]]$spec))
+      {
+        stop("The mmkin objects have to be based on the same degradation models")
+      }
+    }
+    # Check if they have been fitted to the same dataset
+    other_object_ds <- lapply(other[1, ], function(x) x$data)
+    for (i in 1:length(ds)) {
+      if (!all.equal(ds[[i]][c("time", "variable", "observed")],
+          other_object_ds[[i]][c("time", "variable", "observed")]))
+      {
+        stop("The mmkin objects have to be fitted to the same datasets")
+      }
+    }
+  }
+
+  n.o <- length(objects)
+
+  error_models = sapply(objects, function(x) x[[1]]$err_mod)
+  n.e <- length(error_models)
+
+  n.fits <- n.deg * n.e
+  fit_indices <- matrix(1:n.fits, ncol = n.e)
+  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, ]
+    res <- try(do.call(backend_function, args = c(list(mmkin_row), dot_args)))
+    return(res)
+  }
+
+  fit_time <- system.time({
+    if (is.null(cluster)) {
+      results <- parallel::mclapply(as.list(1:n.fits), fit_function,
+        mc.cores = cores, mc.preschedule = FALSE)
+    } else {
+      results <- parallel::parLapply(cluster, as.list(1:n.fits), fit_function)
+    }
+  })
+
+  attributes(results) <- attributes(fit_indices)
+  attr(results, "call") <- call
+  attr(results, "time") <- fit_time
+  class(results) <- "mhmkin"
+  return(results)
+}
+
+#' Subsetting method for mhmkin objects
+#'
+#' @param x An [mhmkin] object.
+#' @param i Row index selecting the fits for specific models
+#' @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}}.
+#' @rdname mhmkin
+#' @export
+`[.mhmkin` <- function(x, i, j, ..., drop = FALSE) {
+  class(x) <- NULL
+  x_sub <- x[i, j, drop = drop]
+
+  if (!drop) {
+    class(x_sub) <- "mhmkin"
+  }
+  return(x_sub)
+}
+
+#' Print method for mhmkin objects
+#'
+#' @rdname mhmkin
+#' @export
+print.mhmkin <- function(x, ...) {
+  cat("<mhmkin> object\n")
+  cat("Status of individual fits:\n\n")
+  print(convergence(x))
+}
+
+#' @export
+convergence.mhmkin <- function(object, ...) {
+  all_summary_warnings <- character()
+
+  result <- lapply(object,
+    function(fit) {
+      if (inherits(fit, "try-error")) return("E")
+      else {
+        return("OK")
+      }
+  })
+  result <- unlist(result)
+  dim(result) <- dim(object)
+  dimnames(result) <- dimnames(object)
+
+  class(result) <- "convergence.mhmkin"
+  return(result)
+}
+
+#' @export
+print.convergence.mhmkin <- function(x, ...) {
+  class(x) <- NULL
+  print(x, quote = FALSE)
+  cat("\n")
+  if (any(x == "OK")) cat("OK: Fit terminated successfully\n")
+  if (any(x == "E")) cat("E: Error\n")
+}
+
+#' @export
+AIC.mhmkin <- function(object, ..., k = 2) {
+  res <- sapply(object, function(x) {
+    if (inherits(x, "try-error")) return(NA)
+    else return(AIC(x$so, k = k))
+  })
+  dim(res) <- dim(object)
+  dimnames(res) <- dimnames(object)
+  return(res)
+}
+
+#' @export
+BIC.mhmkin <- function(object, ...) {
+  res <- sapply(object, function(x) BIC(x$so))
+  dim(res) <- dim(object)
+  dimnames(res) <- dimnames(object)
+  return(res)
+}
diff --git a/R/saem.R b/R/saem.R
index e2b4d21c..370de3d8 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -58,6 +58,9 @@ utils::globalVariables(c("predicted", "std"))
 #' f_saem_sfo <- saem(f_mmkin_parent["SFO", ])
 #' f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
 #' f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
+#' illparms(f_saem_dfop)
+#' update(f_saem_dfop, covariance.model = diag(c(1, 1, 1, 0)))
+#' AIC(f_saem_dfop)
 #'
 #' # The returned saem.mmkin object contains an SaemixObject, therefore we can use
 #' # functions from saemix
@@ -125,6 +128,7 @@ saem.mmkin <- function(object,
   fail_with_errors = TRUE,
   verbose = FALSE, quiet = FALSE, ...)
 {
+  call <- match.call()
   transformations <- match.arg(transformations)
   m_saemix <- saemix_model(object, verbose = verbose,
     degparms_start = degparms_start,
@@ -184,6 +188,7 @@ saem.mmkin <- function(object,
     transform_rates = object[[1]]$transform_rates,
     transform_fractions = object[[1]]$transform_fractions,
     so = f_saemix,
+    call = call,
     time = fit_time,
     mean_dp_start = attr(m_saemix, "mean_dp_start"),
     bparms.optim = bparms_optim,
@@ -579,3 +584,32 @@ saemix_data <- function(object, verbose = FALSE, ...) {
     ...)
   return(res)
 }
+
+#' @export
+logLik.saem.mmkin <- function(object, ...) return(logLik(object$so))
+
+#' @export
+update.saem.mmkin <- function(object, ..., evaluate = TRUE) {
+  call <- object$call
+  # For some reason we get saem.mmkin in the call when using mhmkin
+  # so we need to fix this in order to avoid exporting saem.mmkin
+  # in addition to the S3 method
+  call[[1]] <- saem
+  update_arguments <- match.call(expand.dots = FALSE)$...
+
+  if (length(update_arguments) > 0) {
+    update_arguments_in_call <- !is.na(match(names(update_arguments), names(call)))
+  }
+
+  for (a in names(update_arguments)[update_arguments_in_call]) {
+    call[[a]] <- update_arguments[[a]]
+  }
+
+  update_arguments_not_in_call <- !update_arguments_in_call
+  if(any(update_arguments_not_in_call)) {
+    call <- c(as.list(call), update_arguments[update_arguments_not_in_call])
+    call <- as.call(call)
+  }
+  if(evaluate) eval(call, parent.frame())
+  else call
+}