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

14 files changed, 460 insertions, 24 deletions

diff --git a/NAMESPACE b/NAMESPACE
index 96cd92e4..50cf1486 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,21 +1,29 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method("[",mhmkin)
 S3method("[",mmkin)
+S3method(AIC,mhmkin)
 S3method(AIC,mmkin)
+S3method(BIC,mhmkin)
 S3method(BIC,mmkin)
 S3method(aw,mkinfit)
 S3method(aw,mmkin)
 S3method(confint,mkinfit)
+S3method(convergence,mhmkin)
 S3method(convergence,mmkin)
 S3method(f_time_norm_focus,mkindsg)
 S3method(f_time_norm_focus,numeric)
+S3method(illparms,mhmkin)
 S3method(illparms,mkinfit)
 S3method(illparms,mmkin)
+S3method(illparms,saem.mmkin)
 S3method(intervals,saem.mmkin)
 S3method(loftest,mkinfit)
 S3method(logLik,mkinfit)
+S3method(logLik,saem.mmkin)
 S3method(lrtest,mkinfit)
 S3method(lrtest,mmkin)
+S3method(mhmkin,list)
 S3method(mixed,mmkin)
 S3method(mkinpredict,mkinfit)
 S3method(mkinpredict,mkinmod)
@@ -27,8 +35,11 @@ S3method(plot,mixed.mmkin)
 S3method(plot,mkinfit)
 S3method(plot,mmkin)
 S3method(plot,nafta)
+S3method(print,convergence.mhmkin)
 S3method(print,convergence.mmkin)
+S3method(print,illparms.mhmkin)
 S3method(print,illparms.mmkin)
+S3method(print,mhmkin)
 S3method(print,mixed.mmkin)
 S3method(print,mkinds)
 S3method(print,mkindsg)
@@ -50,6 +61,7 @@ S3method(summary,saem.mmkin)
 S3method(update,mkinfit)
 S3method(update,mmkin)
 S3method(update,nlme.mmkin)
+S3method(update,saem.mmkin)
 export(CAKE_export)
 export(DFOP.solution)
 export(FOMC.solution)
@@ -78,6 +90,7 @@ export(max_twa_hs)
 export(max_twa_parent)
 export(max_twa_sfo)
 export(mean_degparms)
+export(mhmkin)
 export(mixed)
 export(mkin_long_to_wide)
 export(mkin_wide_to_long)
diff --git a/NEWS.md b/NEWS.md
index 93f688ff..4cffcb81 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,10 +1,14 @@
 # mkin 1.1.2
 
+- 'R/mhmkin.R': New method for performing multiple hierarchical mkin fits in one function call, optionally in parallel.
+
 - 'R/saem.R': Implement and test saemix transformations for FOMC and HS. Also, error out if saemix transformations are requested but not supported.
 
-- 'R/convergence.R': New generic to show convergence information with a method for 'mmin' objects.
+- 'R/saem.R': 'logLik' and 'update' methods for 'saem.mmkin' objects.
+
+- 'R/convergence.R': New generic to show convergence information with methods for 'mmkin' and 'mhmkin' objects.
 
-- 'R/illparms.R': New generic to show ill-defined parameters with methods for 'mkinfit' and 'mmkin' objects.
+- 'R/illparms.R': New generic to show ill-defined parameters with methods for 'mkinfit', 'mmkin', 'saem.mmkin' and 'mhmkin' objects.
 
 - 'R/summary.mmkin.R': Summary method for mmkin objects.
 
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
+}
diff --git a/log/test.log b/log/test.log
index 67ba8841..8be4a512 100644
--- a/log/test.log
+++ b/log/test.log
@@ -1,11 +1,11 @@
 ℹ Testing mkin
 ✔ | F W S  OK | Context
 ✔ |         5 | AIC calculation
-✔ |         5 | Analytical solutions for coupled models [3.3s]
+✔ |         5 | Analytical solutions for coupled models [3.4s]
 ✔ |         5 | Calculation of Akaike weights
 ✔ |         3 | Export dataset for reading into CAKE
 ✔ |        12 | Confidence intervals and p-values [1.0s]
-✔ |     1  12 | Dimethenamid data from 2018 [32.9s]
+✔ |     1  12 | Dimethenamid data from 2018 [32.4s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip (test_dmta.R:98:3): Different backends get consistent results for SFO-SFO3+, dimethenamid data
 Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
@@ -16,6 +16,7 @@ Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
 ✔ |        14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.8s]
 ✔ |         4 | Test fitting the decline of metabolites from their maximum [0.4s]
 ✔ |         1 | Fitting the logistic model [0.2s]
+✔ |         5 | Batch fitting and diagnosing hierarchical kinetic models [14.5s]
 ✔ |     1  12 | Nonlinear mixed-effects models [0.2s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip (test_mixed.R:68:3): saemix results are reproducible for biphasic fits
@@ -27,23 +28,23 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve
 ✔ |         8 | mkinmod model generation and printing [0.2s]
 ✔ |         3 | Model predictions with mkinpredict [0.4s]
 ✔ |        16 | Evaluations according to 2015 NAFTA guidance [1.8s]
-✔ |         9 | Nonlinear mixed-effects models with nlme [8.3s]
-✔ |        16 | Plotting [10.6s]
+✔ |         9 | Nonlinear mixed-effects models with nlme [8.7s]
+✔ |        16 | Plotting [10.1s]
 ✔ |         4 | Residuals extracted from mkinfit models
-✔ |        28 | saemix parent models [182.2s]
+✔ |        28 | saemix parent models [181.2s]
 ✔ |         2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s]
-✔ |         7 | Fitting the SFORB model [3.9s]
+✔ |         7 | Fitting the SFORB model [3.7s]
 ✔ |         1 | Summaries of old mkinfit objects
 ✔ |         5 | Summary [0.2s]
 ✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.1s]
-✔ |         9 | Hypothesis tests [7.7s]
-✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
+✔ |         9 | Hypothesis tests [7.8s]
+✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [2.1s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 266.8 s
+Duration: 279.6 s
 
 ── Skipped tests  ──────────────────────────────────────────────────────────────
 • Fitting this ODE model with saemix takes about 15 minutes on my system (1)
 • Fitting with saemix takes around 10 minutes when using deSolve (1)
 
-[ FAIL 0 | WARN 0 | SKIP 2 | PASS 228 ]
+[ FAIL 0 | WARN 0 | SKIP 2 | PASS 233 ]
diff --git a/man/illparms.Rd b/man/illparms.Rd
index 7f70229d..90adf2bb 100644
--- a/man/illparms.Rd
+++ b/man/illparms.Rd
@@ -5,6 +5,9 @@
 \alias{illparms.mkinfit}
 \alias{illparms.mmkin}
 \alias{print.illparms.mmkin}
+\alias{illparms.saem.mmkin}
+\alias{illparms.mhmkin}
+\alias{print.illparms.mhmkin}
 \title{Method to get the names of ill-defined parameters}
 \usage{
 illparms(object, ...)
@@ -14,6 +17,12 @@ illparms(object, ...)
 \method{illparms}{mmkin}(object, conf.level = 0.95, ...)
 
 \method{print}{illparms.mmkin}(x, ...)
+
+\method{illparms}{saem.mmkin}(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...)
+
+\method{illparms}{mhmkin}(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...)
+
+\method{print}{illparms.mhmkin}(x, ...)
 }
 \arguments{
 \item{object}{The object to investigate}
@@ -23,14 +32,32 @@ illparms(object, ...)
 \item{conf.level}{The confidence level for checking p values}
 
 \item{x}{The object to be printed}
+
+\item{random}{For hierarchical fits, should random effects be tested?}
+
+\item{errmod}{For hierarchical fits, should error model parameters be
+tested?}
 }
 \value{
-For \link{mkinfit} objects, a character vector of parameter names
-For \link{mmkin} objects, an object of class 'illparms.mmkin' with a
-suitable printing method.
+For \link{mkinfit} or \link{saem} objects, a character vector of parameter
+names. For \link{mmkin} or \link{mhmkin} objects, a matrix like object of class
+'illparms.mmkin' or 'illparms.mhmkin'. The latter objects have a suitable
+printing method.
 }
 \description{
-Method to get the names of ill-defined parameters
+The method for generalised nonlinear regression fits as obtained
+with \link{mkinfit} and \link{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.
+}
+\details{
+The method for hierarchical model fits, also known as nonlinear
+mixed-effects model fits as obtained with \link{saem} and \link{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.
 }
 \examples{
 fit <- mkinfit("FOMC", FOCUS_2006_A, quiet = TRUE)
diff --git a/man/mhmkin.Rd b/man/mhmkin.Rd
new file mode 100644
index 00000000..0ef1599e
--- /dev/null
+++ b/man/mhmkin.Rd
@@ -0,0 +1,75 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mhmkin.R
+\name{mhmkin}
+\alias{mhmkin}
+\alias{mhmkin.list}
+\alias{[.mhmkin}
+\alias{print.mhmkin}
+\title{Fit nonlinear mixed-effects models built from one or more kinetic
+degradation models and one or more error models}
+\usage{
+mhmkin(objects, backend = "saemix", algorithm = "saem", ...)
+
+\method{mhmkin}{list}(
+  objects,
+  backend = "saemix",
+  ...,
+  cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(),
+  cluster = NULL
+)
+
+\method{[}{mhmkin}(x, i, j, ..., drop = FALSE)
+
+\method{print}{mhmkin}(x, ...)
+}
+\arguments{
+\item{objects}{A list of \link{mmkin} objects containing fits of the same
+degradation models to the same data, but using different error models.}
+
+\item{backend}{The backend to be used for fitting. Currently, only saemix is
+supported}
+
+\item{algorithm}{The algorithm to be used for fitting (currently not used)}
+
+\item{\dots}{Further arguments that will be passed to the nonlinear mixed-effects
+model fitting function.}
+
+\item{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 \code{\link[parallel:detectCores]{parallel::detectCores()}} are used, except on Windows where the default
+is 1.}
+
+\item{cluster}{A cluster as returned by \link{makeCluster} to be used for
+parallel execution.}
+
+\item{x}{An \link{mhmkin} object.}
+
+\item{i}{Row index selecting the fits for specific models}
+
+\item{j}{Column index selecting the fits to specific datasets}
+
+\item{drop}{If FALSE, the method always returns an mhmkin object, otherwise
+either a list of fit objects or a single fit object.}
+}
+\value{
+A two-dimensional \link{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'.
+
+An object of class \code{\link{mhmkin}}.
+}
+\description{
+The name of the methods expresses that (\strong{m}ultiple) \strong{h}ierarchichal
+(also known as multilevel) \strong{m}ulticompartment \strong{kin}etic models are
+fitted. Our kinetic models are nonlinear, so we can use various nonlinear
+mixed-effects model fitting functions.
+}
+\seealso{
+\code{\link{[.mhmkin}} for subsetting \link{mhmkin} objects
+}
+\author{
+Johannes Ranke
+}
diff --git a/man/saem.Rd b/man/saem.Rd
index 0c066dd2..a202f52f 100644
--- a/man/saem.Rd
+++ b/man/saem.Rd
@@ -49,9 +49,9 @@ saemix_data(object, verbose = FALSE, ...)
 
 \item{transformations}{Per default, all parameter transformations are done
 in mkin. If this argument is set to 'saemix', parameter transformations
-are done in 'saemix' for the supported cases. Currently this is only
-supported in cases where the initial concentration of the parent is not fixed,
-SFO or DFOP is used for the parent and there is either no metabolite or one.}
+are done in 'saemix' for the supported cases, i.e. (as of version 1.1.2)
+SFO, FOMC, DFOP and HS without fixing \code{parent_0}, and SFO or DFOP with
+one SFO metabolite.}
 
 \item{degparms_start}{Parameter values given as a named numeric vector will
 be used to override the starting values obtained from the 'mmkin' object.}
@@ -120,6 +120,9 @@ f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE)
 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
diff --git a/tests/testthat/convergence_fits_synth_const.txt b/tests/testthat/convergence_fits_synth_const.txt
new file mode 100644
index 00000000..e69de29b
--- /dev/null
+++ b/tests/testthat/convergence_fits_synth_const.txt
diff --git a/tests/testthat/convergence_hfits_synth.txt b/tests/testthat/convergence_hfits_synth.txt
new file mode 100644
index 00000000..e69de29b
--- /dev/null
+++ b/tests/testthat/convergence_hfits_synth.txt
diff --git a/tests/testthat/illparms_hfits_synth.txt b/tests/testthat/illparms_hfits_synth.txt
new file mode 100644
index 00000000..e69de29b
--- /dev/null
+++ b/tests/testthat/illparms_hfits_synth.txt
diff --git a/tests/testthat/setup_script.R b/tests/testthat/setup_script.R
index 62cdf78e..af059668 100644
--- a/tests/testthat/setup_script.R
+++ b/tests/testthat/setup_script.R
@@ -8,8 +8,8 @@ if (identical(Sys.getenv("NOT_CRAN"), "true")) {
   n_cores <- 1
 }
 
-# We are only allowed one core on travis, but they also set NOT_CRAN=true
-if (Sys.getenv("TRAVIS") != "") n_cores = 1
+# Use the two available cores on travis
+if (Sys.getenv("TRAVIS") != "") n_cores = 2
 
 # On Windows we would need to make a cluster first
 if (Sys.info()["sysname"] == "Windows") n_cores = 1
diff --git a/tests/testthat/test_mhmkin.R b/tests/testthat/test_mhmkin.R
new file mode 100644
index 00000000..5243f971
--- /dev/null
+++ b/tests/testthat/test_mhmkin.R
@@ -0,0 +1,31 @@
+context("Batch fitting and diagnosing hierarchical kinetic models")
+
+test_that("Multiple hierarchical kinetic models can be fitted and diagnosed", {
+
+  skip_on_cran()
+  fits_synth_const <- suppressWarnings(
+    mmkin(c("SFO", "FOMC"), ds_sfo[1:6], cores = n_cores, quiet = TRUE))
+
+  fits_synth_tc <- suppressWarnings(
+    update(fits_synth_const, error_model = "tc"))
+
+  expect_known_output(
+    convergence(fits_synth_const),
+    "convergence_fits_synth_const.txt")
+
+  hfits <- mhmkin(list(fits_synth_const, fits_synth_tc))
+
+  expect_known_output(
+    convergence(hfits),
+    "convergence_hfits_synth.txt")
+
+  expect_known_output(
+    illparms(hfits),
+    "illparms_hfits_synth.txt")
+
+  expect_equal(which.min(AIC(hfits)), 3)
+
+  hfit_sfo_tc <- update(hfits[["SFO", "tc"]],
+    covariance.model = diag(c(0, 1)))
+  expect_equal(illparms(hfit_sfo_tc), character(0))
+})