aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--NAMESPACE13
-rw-r--r--NEWS.md8
-rw-r--r--R/illparms.R70
-rw-r--r--R/mhmkin.R184
-rw-r--r--R/saem.R34
-rw-r--r--log/test.log21
-rw-r--r--man/illparms.Rd35
-rw-r--r--man/mhmkin.Rd75
-rw-r--r--man/saem.Rd9
-rw-r--r--tests/testthat/convergence_fits_synth_const.txt0
-rw-r--r--tests/testthat/convergence_hfits_synth.txt0
-rw-r--r--tests/testthat/illparms_hfits_synth.txt0
-rw-r--r--tests/testthat/setup_script.R4
-rw-r--r--tests/testthat/test_mhmkin.R31
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))
+})

Contact - Imprint